In preparation for some work I will be doing I've been working on calculation of insolation as a function of time and latitude. In this post I will explain the method and compare results to a readily obtainable plot of insolation at various latitudes.
Note, all of the following calculations and the comparison with another source are for top of atmosphere and do not take into account absorption/scattering by the atmosphere, or reflection from clouds.
Insolation is a term for the solar radiation received by the earth, most of the insolation is in the visible spectrum, indeed that is the reason our eyes see in that spectrum (hence visible spectrum). The calculations have been carried out in an Excel 2007 spreadsheet, this should work in earlier and later versions of Excel, whether it will work in other programs (such as Libre Office) I do not know. The spreadsheet is available here. The maths and theory behind the calculation has been obtained from these lecture notes. Thanks are due once again to an academic making their course material available to the public!
Before going into detail, here is a screenshot of the upper left corner of the spreadsheet. Note that the statement 'assume circular orbit' is something I forgot to remove from the initial version of the calculation, an elliptical orbit is assumed, as will be explained.
I've colour coded it so as to try to make the various elements more clear. The largest single region, ranging from greens to reds is the actual results of calculated insolation. To the left of that is a pink region with latitude in degrees and in radians.
The constants used are at the top of the sheet, tilt, orbital period, eccentricity of the earth, also solar constant (1361W/m^2) and vernal equinox (20 March).
The year day is given in five day increments, it's easy to change that if needed because of the way the formula has been implemented. Above that is the ratio of the average earth to sun distance (ro) and the actual distance factoring in the eccentricity of the orbit of the earth. The vernal day is the day of a year centred on the spring equinox (20 March ~ Day 80). And the declination gives the variation of the angle of the sun throughout the year (in radians) due to the tilt of the earth.
Now for the nuts and bolts.
The length of day is H/Pi (I don't have the easy facility to use greek notation here, the English version of such notation is italicised from here on). H itself is the hour angle at sunset. H is calculated using latitude (phi) and the declination of the sun (delta), presented in terms of Excel formulae the formula for H is:
H = ACOS(-TAN(phi)*TAN(delta))
And here we hit a problem, in the Arctic the sun never rises for periods in the winter, and never sets for periods in the summer, the result is that the -TAN(phi)*TAN(delta) produces unphysical results which are greater than +/-1 so the ACOS (inverse cosine) generates a numeric error at high latitudes in summer and winter. A cosine always returns a result between +1 and -1 so the inverse cosine cannot be fed with a number outside that range.
H = ACOS(
Already that looks a bit daunting, but it isn't too bad really. ABS() is a function that produces the absolute numeric value, without the sign of the number. SIGN() is a function that produces a +1 if the sign of the number is positive, and a -1 if the sign of the number is negative. So the above is described by:
If ABS(-TAN(phi)*TAN(delta)) is greater than 1 then use the sign of the "-TAN(phi)*TAN(delta)" part, otherwise use the actual value "-TAN(phi)*TAN(delta)".
This means the result is always equal to or less than 1 and stops the numeric errors at the polar regions.
From the lecture notes linked to above the following equation is provided, it describes the effect of the elliptical orbit of the earth on the earth to sun distance, which obviously impacts received insolation.
This gives the basics of a term (ro/r)^2 in the main equation, described below. Using the latter of the above two equations I use the following for the (ro/r)^2 term:
(ro/r)^2 = 1/(1-Eccentricity*COS(2*PI()*((YearDay-4)/OrbitalPeriod)))
Note that the average distance between earth and sun (ro) and the distance adjusted for the ellipse of the earth's orbit (r) are not actually used, because both 'r' terms cancel out. The result is the ratio of the two which is all we need.
So from the lecture notes linked to at the top of this post we have the main equation.
We have the term 'So' (Solar Constant) as the insolation from the sun received in space at the disantance from earth to sun perpendicular to the sun. 1361W/m^2 is used for this, sourced from the Wikipedia page on Solar Constant. We have the (ro/r)^2 term, which is implemented as a row in the spreadsheet (row 7), term H is as described above and is implemented within the equation in each spreadsheet cell where insolation is calculated. Latitude (phi) is calculated in radians and provided in column E of the spreadsheet. Declination (delta) is calculated and provided in row 10 of the spreadsheet.
The variable setting rows and columns then provide the drivers for a rectangular matrix of calculation of insolation as a function of latitude (vertical) and time (horizontal). The equation has been set up so that copying and pasting, also fill down, work even if the size of the range is changed to increase or decrease time or latitudinal resolution (using the $ operator in cell notation).
The main equation as used in the sheet is:
The 'H' term being highlighted in blue to aid reading.
So how does it all work out? A commonly used presentation of insolation as a function of time and latitude is available from PhysicalGeography.net. This plot uses a set of latitude bands.
Using that as a reference of the correct answer worked out by a proper scientist (I presume), I can compare it with my own amateur calculation.
The results are pretty close, but not exactly (is anything ever so?). Note that for 60 and 30degN in my calculation the insolation is the same in mid summer, whereas in the official plot, 60degN is slightly lower. Also the mid summer peak for 90degN is too high by a small margin in my plot. However overall the agreement is good, and for my purposes what looks like a few percent of error at high latitudes is not likely to be a problem.