aa r X i v : . [ phy s i c s . pop - ph ] S e p Derivation of Solar Position Formulae
Ross Ure Anderson ∗ Abstract.
Derivation of the following formulae for solar position as seen from orbitingplanet based on a simplified model: sunrise direction formula, solar declination for-mula, sunrise equation, daylight duration formula, solar altitude formula, solar azimuthformula. Use of notion of effective axial tilt, and Rodrigues Rotation Formula, and re-flections of the orbital quadrants to reduce the general case to the simpler case of the dayof the winter solstice. Derivation of equations for solar time to clock time conversion,and sunrise, sunset, and solar noon times. Implementation of an analemma calculator.Comparison of the sunrise direction formula with 304 point dataset of actual sunrisedata from Earth, obtaining average accuracy of 1 . ◦ , and estimate of Earth’s axial tiltof 23 . ◦ — within 0 . ◦ of the currently accepted value 23 . ◦ . This article first derives a formula for the sunrise direction for a planet orbiting a central sun, in termsof the day of the year and the latitude, under a simplified model. Then, building on the method of proofused a number of further formulae relating to the solar position [1], [2] are derived. The approach adoptedreduces the cases of the four orbital quadrants to a single quadrant and then reduces the latter case to asingle point (the winter solstice). The simplifying assumptions made throughout are :1. The planet is a sphere2. The orbit of the planet is an ellipse with the sun at one focus
3. The planet radius is negligible compared with its minimal distance from the sun4. Whilst orbiting, the planet rotates about a fixed direction axis through its center5. The sun can be approximated as a point source of light6. A day of the year is approximated as a single stationary point in the orbit at which the planet completesa full 360 ◦ revolutionThe sunrise direction formula gives the sunrise direction θ as an angle north of east : θ = − arcsin (cid:18) sin α cos ψ cos δ (cid:19) (1)where : θ ∈ [ − ◦ , ◦ ], α = planet axial tilt ∈ [0 ◦ , ◦ ], δ = latitude ∈ ( − ◦ , ◦ ), ψ = day of year angle = orbital angle swept out from winter solsticeIn a circular orbit/constant speed model the stationary days are evenly spaced around the circle andthe orbital angle ψ can be defined as dN (360 ◦ ) for day d ∈ [0 , N − d is the day offset from theday of the winter solstice, and N is the number of days in the year. A more accurate relation between theday of the year and ψ , in the case of Earth, is described in [3] with a formula expressing the angle sweptout from the spring equinox as a function of time — we can then take this angle at say 12 noon on a day and ∗ E-mail: [email protected]; web: archive.org/details/@ross ure anderson The proofs below do not specifically make use of the elliptical property — the only orbital parameter they require is theorbital angle ψ , thus they would apply to any planar orbit shape so long as the other simplifying assumptions applied. In thecase of Earth the gravity of the moon and other planets perturbs its orbit away from a perfect ellipse. ◦ to obtain ψ , the time in the formula being measured in units of mean solar days of 24 hours duration from UTCmidnight on 1 st January 2013. In the comparison with real data from Earth in Appendices A.1 and B the former modelis used. To use the latter model instead, the program code which computes the formula’s predictions can be changed asdescribed in Appendix C. The former model is a reasonable approximation for Earth whose orbit is very close to circular,but for a planet with more eccentric orbit a procedure such as in [3] would be required.In the above simplified model the sunset position is symmetrically opposite to the sunrise position, on the west sideof the horizon, at angle θ north of west. A negative θ means south of east/west. An alternative expression is : θ = − arcsin (cid:18) sin α cos δ (cid:19) (2)where sin α = sin α cos ψ (3)defines an ‘effective axial tilt’ α at orbital angle ψ (Figure 5), which makes the geometric situation identical to thewinter solstice day of ψ = 0 ◦ , when viewed from a different angle. The effective axial tilt equals the actual axial tilt α onthe winter solstice, and falls to zero at the spring equinox, so α ∈ [0 , α ]. Any day of the first quadrant Q α replacing α . Days in Q Q Q Q
1. Thuswith the substitution of α for α the simpler case of the winter solstice in Figure 1 applies to the entire orbit, with thesun being placed on the opposite side of that in Figure 1 in the cases of Q Q
3. The effective axial tilt is definedonly for Q , α ] but if the definition (3) is extended throughout the whole year to a value α in therange [ − α, α ] then it equals minus the solar declination λ ( § θ cos δ + sin α cos ψ = 0 (4)or, to estimate the axial tilt from the sunrise direction (Appendix B) :sin α = − sin θ cos δ cos ψ (5)provided cos ψ = 0, ie we are not at an equinox, where the formula becomes 0 /
0. Similarly (4) can be rearranged toexpress day of the year in terms of sunrise position at a known latitude (or Pole Star elevation, Appendix E) if the axialtilt is known.The term ‘winter solstice’ refers to the standard notion of ‘winter’ of the northern hemisphere, even though towardsthe equator the terms winter and summer lose their normal meanings — for example on the equator maximal solarradiation is received at the equinoxes and is minimal at both summer and winter solstices [4].Though astronomical terminology is used throughout, the problems of the solar position under the simplified modelare purely geometric involving a sphere on an elliptical path intersected by parallel rays from the direction of one of thepath’s foci. For example the point of sunrise/sunset on a planet is simply when these rays become parallel to the tangentplane at the point. The desired sunrise/sunset direction is then the direction of the rays wrt true north [5] in the localtangent (or horizon) plane.In Appendix A.1, the sunrise direction formula with a circular orbit/constant speed model is compared with actualdata for the case of the Earth, using 304 data points at 8 latitudes for the year 2018–2019 — the actual sunrise directionsbeing taken from [6] and the latitudes from
Google Maps . Using an Earth axial tilt of 23 . ◦ [7],the formula then matches the actual sunrise directions to within an average error of 1 . ◦ , which is a quite good ap-proximation for the simplified model. Using the formula’s own predicted axial tilt of 23 . ◦ (Appendix B) the error is 1 . ◦ .A chief interest of the sunrise direction formula is how in relating the above quantities measurable by naked eye fromthe Earth using only primitive instruments, the hypothesis of a spherical Earth in a circular orbit and rotating about afixed axis is arrived at using only a geometric argument — since the close correlation with the actual data would be highlyunlikely to occur unless the above simplified model was a good approximation in practice. Thus it is one method fromwhich the ancients could have determined the Earth shape and motions from only naked eye observations and geometry,without the benefit of modern astronomy.The naked eye measurements corresponding to the above parameters are :2A) n , the day offset from the winter solstice, where the winter solstice is determined by the day on which the Sun isat its lowest southerly elevation at solar noon (with the convention of a northerly Sun elevation being > ◦ ), and n is counted from that day starting at zero. For southerly latitudes such lowest southerly elevation would actuallybe in the summer as there the Sun reaches into the northern half of the sky and is low in the northern sky in thewinter — the day offset would thus be counted from the summer solstice when the Sun is at its highest northerly elevation at solar noon, which is at the same time as the northern winter solstice.(B) δ , the latitude, is obtained from the fixed angle of elevation of the Pole Star at the given location. Where the PoleStar is no longer visible a corresponding fixed point in the southern sky would be used to determine latitude.(C) θ , the position of the sunrise, is measured wrt due east, where due east is determined wrt true north, and true northis determined from the position of the Pole Star projected onto the horizon orthogonally (the ‘azimuthal’ position).The latitude δ being measured as elevation of the Pole Star is based on the hypothesis of a spherical Earth and verydistant Pole Star (see Appendix E). Measurement (A) may be recorded by some kind of calendrical system markingmid-winter and mid-summer.Originally ancient astronomers would only have these observed quantities to work with, gathering data over longperiods of time, and hypothesizing what may be the larger scale structure that caused them [8], [9], [10], [11].The formula would then relate the above quantities with the right choice of axial tilt, which could be hypothesizedby some other method or be estimated from actual sunrise data using equation (5) (Appendix B). The constancy of therhs of equation (5) for actual data would give evidence for the hypotheses 1 – 6 above.In short the above mathematical relationship between the observable quantities (A) – (C) gives a strong evidence forthe hypotheses 1 – 6 with a circular orbit describing, to a good approximation, the Earth shape and motion throughspace — since it would be very difficult for any other configuration to produce the same correlation.The derivation below of the sunrise direction formula envisages a prograde (wrt orbital direction) rotation axis oftilt up to 90 ◦ . For a retrograde rotation ( α > ◦ ) the axial tilt angle (180 − θ ) ◦ could be used, and would be clock-wise, so the sunrise and sunset positions would simply be swapped, and the formula would remain applicable, sincesin α = sin(180 − α ). The ‘north’ side of the orbital plane is defined as the side from which the orbit appears anti-clockwise. In using (5) to estimate axial tilt, the result chosen is α ≤ ◦ or α > ◦ according as the rotation is progradeor retrograde, the latter being when the sun rises in the west and sets in the east. From § α ∈ [0 ◦ , ◦ ).The sunrise direction formula applies to all latitudes δ except for the poles, at 90 ◦ and − ◦ , where the sunrise/sunsetis caused by the orbital motion of the planet rather than the daily rotation — eg on Earth at the North Pole the Sunrises once per year on the spring equinox and sets once per year on the autumnal equinox, the opposite being true forthe South Pole, and this type of sunrise takes much longer than a diurnal sunrise, around 30 hours versus a few minutes,for the sun’s disc to fully cross the horizon, and the type of discrete sunrise direction we consider here is not applicable[12], [13], [14], [15], [16].For non-pole latitudes beyond one of the arctic circles such non-diurnal sunrises and sunsets also occur at certaintimes of the year — when periods of perpetual day or night are entered. Outside those periods the sunrise and sunsetare diurnal. These changes occur as the effective axial tilt changes so the actual latitude travels above and below theeffective arctic circle (the latter becoming a point at the pole on the spring equinox). In the sunrise direction formulathe periods of perpetual day and night correspond with the arcsin argument going out of range, for there is no sunsetnor sunrise at these times — and this can be used to determine the range of dates for these periods for a given polarlatitude, as shown in some examples in Appendix A.2. 3 Sunrise Direction Formula
The formula is first proved for the simplest case of the winter solstice (ie ψ = 0), from which the summer solstice casethen readily follows. The equinoxes are readily checked because formula (1) then gives θ = 0 at all latitudes which isclear from the symmetry — the effective axial tilt is zero, ie. vertical, and each latitude has equal lengths of day andnight, the sunset and sunrise being along the east/west line everywhere (except the poles which are seeing the non-diurnalsunset/sunrise). The general case for the first quadrant Q Q Q
3, and Q Q We require to show : θ = − arcsin (cid:18) sin α cos δ (cid:19) (6) The winter solstice case ( ψ = 0) is shown in Figures 1 and 2, with latitude δ ∈ [0 , − α ] in the northern hemisphere.The case of δ = 90 − α means points A and B coincide and formula (6) reduces to − arcsin 1 = −
90, ie southerly, as alsoseen from Figure 1. The case of δ = 0 means B coincides with O in Figure 1 and (6) reduces to the required angle of θ = − α .The east-west line at any point is the intersection of the latitude plane and the horizon plane at that point. Thus toobtain the east-west line we can take the cross product of the normals n L , n H to these two planes.From Figure 1, we can choose n L = (sin α, cos α, n H = (0 , BC, OC ), so : n L × n H = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i j k sin α cos α BC OC (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = ( OC cos α, − OC sin α, BC sin α ) , ∴ parallel to (cos α, − sin α, BCOC sin α )From Figure 1 the easterly direction at point B must have a positive x -coordinate and a negative y -coordinate andthus we can choose the following as east and west (non-unit) vectors at B : e = (cos α, − sin α, BCOC sin α ) , w = ( − cos α, sin α, − BCOC sin α )Then taking the mirror image of these in the xy -plane, which reverses the z -component, with east mapping onto westand vice-versa, east/west vectors at B ′ are : e ′ = ( − cos α, sin α, BCOC sin α ) , w ′ = (cos α, − sin α, − BCOC sin α )From the assumptions 3 and 5 the rays from the sun striking every point of the planet are parallel and come from the s = − i direction, and it is clear geometrically that the direction of s in the horizon plane at B ′ is south of the easterly e ′ at B ′ . Thus the desired magnitude of angle θ of sunrise satisfies : e ′ · s = | e ′ | | s | cos θ, ∴ cos θ = cos α | e ′ | But | e ′ | = 1 + BC OC sin α and so we need to express BC/OC in terms of α and δ . From Figure 2, OC = √ R − BC , and BC equals thedistance OB from Figure 1. From the triangle △ AOB of Figure 1, we have :sin(90 − α − β ) OB = sin(90 + α ) R , ie. cos( α + β ) OB = cos αR OR DAYHEMISPHERE NIGHTHEMISPHERE x TERMINATOR y z = ⊙ E Q UA T O R A R CT I C C I R C L E AN T A R CT I C C I R C L E P O L A R AX I S A L A T I T U D E αβα BS N Latitude δ = 90 ◦ − ( α + β ),0 < β < ◦ − α , β < ◦ ⇒ perpetual night, β > ◦ − α ⇒ perpetual dayThus BC of Figure 2 satisfies : BC = R cos( α + β )cos α = R sin δ cos α , and so : OC = R s(cid:18) − sin δ cos α (cid:19) = R cos α p cos α − sin δ ⇒ BCOC = sin δ p cos α − sin δ ⇒ | e ′ | = 1 + BC OC · sin α = (cos α − sin δ ) + sin δ sin α cos α − sin δ = cos α − sin δ cos α cos α − sin δ = cos α cos δ cos α − sin δ ⇒ cos θ = p cos α − sin δ cos δ ⇒ sin θ = 1 − (cid:18) cos α − sin δ cos δ (cid:19) = 1 − cos α cos δ = sin α cos δ ⇒ sin θ = sin α cos δ . The special case of α = 90 ◦ can be checked separately using a plan view. θ above was an angle south of east the required angle north of east is given by equation (6).Figure 2: Winter Solstice - View Towards Sun R T E R M I N A T O R z y x = ⊙ AB B ′ L A T I T U D E H O R I Z O N P L AN E C O BAB ′ = latitude circle, A = solar noon, B = sunset, B ′ = sunrise, x > ⇒ night, x < ⇒ day The case δ ∈ [ α − ,
0) in the southern hemisphere is shown in Figures 3 and 4.From the triangle △ AOB of Figure 3, we have :sin( α + β − OB = sin(90 − α ) R ie. − cos( α + β ) OB = cos αR Thus the y -coordinate BC of Figure 4, which is negative of OB of Figure 3, satisfies : BC = R cos( α + β )cos α = R sin δ cos α , as before — thus we obtain the same east-west line as above and since the easterly direction at point B still musthave a positive x -coordinate and a negative y -coordinate the same expressions for e ′ and w ′ , and hence θ , are obtainedas in § α + β > ◦ O xy z = ⊙ P O L A R AX I S A L A T I T U D E αβ BS N Figure 4: Winter Solstice - View Towards Sun - Case α + β > ◦ z y x = ⊙ AB B ′ L A T I T U D E H O R I Z O N P L AN E C O .2 Summer Solstice Here ψ = 180 ◦ so we require to show θ is minus what it is in the winter solstice case. The situation is as Figure 1 withthe sun at the opposite side so s = i , and then it is clear a sunset θ south of west becomes a sunrise θ north of east, anda sunrise θ south of east becomes a sunset θ north of west, as required. Consider a day in the first quadrant Q xyz -axes are fixed to the planet, but have constantdirection as the planet orbits (note this is a different xyz -axis system than was used in Figures 1–4 of § N S axis by some angle ρ about vector b into plane Γ the effective rotation axis N ′ S ′ isproduced, making effective axial tilt angle α with the vertical. This is equivalent to simply viewing the sun and planetfrom a different angle, by rotating ourselves backwards by ρ from plane Γ. After that rotation the situation is thenexactly the same as the winter solstice of Figure 1, except with the axial tilt angle α instead of α . Thus the requiredangle θ is given by equation (6) with α replaced by α . Thus to obtain formula (1) we need to show sin α = sin α cos ψ .Let a = sin α i + cos α k be the N S axis unit vector, and let a ′ be the rotated N ′ S ′ axis unit vector. Then the requiredangle α satisfies : cos α = a ′ · k Wrt right hand rule the rotation about b = cos ψ i + sin ψ j is by − ρ , thus from the Rodrigues Rotation Formula(Appendix D) : a ′ = (cos ρ ) a + (sin ρ ) ( a × b ) + (1 − cos ρ ) ( b · a ) b Viewing along the b axis the rotation angle ρ satisfies : a ⊥ · k = | a ⊥ | cos ρ where a ⊥ is the component of a perpendicular to b . Since a · b = sin α cos ψ , a ⊥ = a − ( a · b ) b = a − (sin α cos ψ ) b = sin α sin ψ i − sin α sin ψ cos ψ j + cos α k Thus | a ⊥ | = sin α sin ψ + sin α sin ψ cos ψ + cos α = sin α sin ψ + cos α = 1 − sin α cos ψ Since a ⊥ · k = cos α we then have (noting the case of the denominator zero ( ⇒ ψ = 0) has already been covered) :cos ρ = cos α p − sin α cos ψ and, sin ρ = 1 − cos α − sin α cos ψ = 1 − sin α cos ψ − cos α − sin α cos ψ = sin α sin ψ − sin α cos ψ Using a × b = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i j k sin α α cos ψ sin ψ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = − cos α sin ψ i + cos α cos ψ j + sin α sin ψ k and the k -components of the above Rodrigues formula for a ′ we obtain a ′ · k = cos ρ cos α + sin ρ sin α sin ψ Q Q Q Q ψ xyz P l a n e Γ b S NαS ′ N ′ α PLANET O ψ = day of year angle swept out, Γ = vertical plane through sun and planet, N ′ S ′ = N S axis rotated aboutvector b into plane Γ , α = true axial tilt (in xz -plane), α = effective axial tilt (in plane Γ )Then, using the above expressions for cos ρ and sin ρ ,cos α = ( a ′ · k ) = (cos α + sin α sin ψ ) − sin α cos ψ = (1 − sin α cos ψ ) − sin α cos ψ = 1 − sin α cos ψ ∴ sin α = sin α cos ψ ∴ sin α = sin α cos ψ as required.For quadrant Q
4, the situation at day of year angle (360 − ψ ) ◦ is exactly the same as for Q ψ ,with the same effective axial tilt α , just on the opposite side of the solstice axis. Thus θ must be the same, and then ascos(360 ◦ − ψ ) = cos ψ , the formula (1) gives the correct result for Q Q
3, the situation at day of year angle (180 + ψ ) ◦ is exactly the same as for Q ψ except with the sun in the exact opposite direction, and thus as with the case of deriving the summer solstice case fromthe winter solstice case above the sunrise direction θ just changes sign. Since cos(180 ◦ + ψ ) = − cos ψ , formula (1) thengives the correct sunrise angle for Q Q
2, the situation at day of year angle (180 − ψ ) ◦ is exactly the same as for Q ψ ) ◦ ,just on the opposite side of the solstice axis. Thus θ must be the same, and then as cos(180 ◦ − ψ ) = cos(180 ◦ + ψ ),formula (1) gives the correct sunrise angle for Q For a given α and δ the equation (1) in radians has form f ( x ) = − sin − ( k cos x ), for x ∈ [0 , π ], where k = sin α/ cos δ ,and k ∈ (0 ,
1) for non-polar latitudes (ie | δ | < π − α ). Then : f ′ ( x ) = k sin x √ (1 − k cos x )and f ′′ ( x ) = k (1 − k ) cos x (1 − k cos x ) / Thus the curve is ‘smooth’ at non-polar latitudes, with the following ‘bell shape’ gradient : x π/ π π/ πf ′ ( x ) 0 + ↑ k + ↓ ↓ − k – ↑ δ = ± ( π − α )), k = 1 and so in Q x ∈ [0 , π ]), f ( x ) = − sin − (cos x ) = x − π , ie a linearfunction. Then by reflecting into Q Q
3, and Q . ◦ .Beyond the polar circles, k > f ( x ) is undefined at certain times of the year when there is perpetual day or night,namely where | cos x | > /k . Figure 6: Sunrise Direction Curves0 90 180 270 360 − − Orbital Position ψ Sun r i s e D i r ec t i o n θ SpringEquinox AutumnalEquinoxSummer Solstice x ◦ f ( x ) ◦ Equatorial sunriseSunrise at latitude 64 ◦ Polar sunrise10
Solar Declination Formula
The solar declination angle [19] λ is the angle made by the direction of the sun from the planet’s center with the planet’sequatorial plane, one of the two coordinates of a celestial body within the equatorial coordinate system [20]. FromFigure 1 it is clear that λ = − α on the winter solstice. And throughout Q α replaced with the effective axial tilt α , we must have λ = − α . Thus in Q λ ∈ [ − α,
0] and :sin λ = − sin α cos ψ (7)In Q
4, at angle (360 − ψ ) ◦ , by symmetry λ is the same as at angle ψ in Q
1, and as cos(360 − ψ ) = cos ψ , formula (7)thus gives the correct result for Q ψ is substituted by (360 − ψ ) ◦ .In Q ψ ) ◦ the sun is in the exact opposite direction in space wrt planet’s center as it was in Q ψ , andthe solar declination is − λ ∈ [0 , α ]. But from (7) :sin( − λ ) = − sin α cos(180 + ψ )so that (7) gives the correct solar declination for Q Q
2, at (180 − ψ ) ◦ , by symmetry solar declination is as Q ψ ) ◦ , ie. − λ ∈ [0 , α ]. But from (7) :sin( − λ ) = − sin α cos(180 − ψ ) , thus (7) gives the correct solar declination for Q λ = − arcsin(sin α cos ψ ) , λ ∈ [ − α, α ] , (8)and the sunrise direction formula can be written as : θ = arcsin (cid:18) sin λ cos δ (cid:19) , θ ∈ [ − ◦ , ◦ ] (9)The formulae in the sections below are expressed in terms of λ and equations (7) and (8) relate λ to the day of the yearvia ψ . The relationship of ψ to day of the year depends on the orbital model used, for example the two orbital modelsconsidered in § Consider a day in Q α replaced with the effective axial tilt α . The latitude circle isshown face on in Figure 7. A is the point of solar noon on the latitude circle, B is the point of sunset, and B ′ is the pointof sunrise. The angle σ is the sunrise (and sunset) solar hour angle [21] which is a measure of the time between the sunrise(or sunset) and solar noon [22]. The time of solar noon is uniquely defined for any day of the year and for any arbitrarylocation L on the planet except for the poles which are stationary throughout the day (in the simplified model). It isdefined as the time when the unique half plane containing the location L and hinged on the N S axis intersects the sun —solar midnight is the time when this half plane is in the exact opposite direction from that. The general angular positionof the half plane through L is a measure of time at L throughout the day as the planet rotates at a constant angularvelocity about the N S axis. The solar hour angle τ of L at a given time is the angle of L ’s half plane wrt the solar noonposition, and is an angle in the range [0 ◦ , ◦ ), with the +ve angle direction defined by the RH rule wrt axis N S . Thesolar hour angle is 0 ◦ at solar noon and 180 ◦ at solar midnight, and is negative at sunrise and positive at sunset. If thesun were placed in the exact opposite direction then the solar hour angle wrt the new sun position would be the previoussolar hour angle plus 180 ◦ (this is used in § − δ in the southern hemisphere is obtained from the latitude + δ by rotating Figure 1 about O by 180 ◦ so in Figure 7 the day/night arc lengths in (ii) are just the reverse of those in (i). Considering case (i) of northernhemisphere the angle σ is in the range [0 , ◦ ] depending on α and δ , eg it is 90 ◦ for all δ at the spring equinox (12 hoursof day and 12 hours of night) and is 0 ◦ at the artic circle on the winter solstice where solar noon, sunrise, and sunset allcoincide. If R L is radius of latitude circle then from Figure 1 : R L R = sin( α + β ) = cos δ ⇒ R L = R cos δ To define solar noon at a pole a meridian longitude would have to be selected for the pole (see § B R L B ′ A d dσσ N I G H T D A Y B R L B ′ A d dσσ N I G H T D A Y (i) Latitude δ ≥
0, Northern Hemisphere (ii) Latitude − δ <
0, Southern Hemisphere A = solar noon, B = sunset, B ′ = sunriseand from Figure 7 we have sin σ = d/R L . But from Figure 2, d = OC , and from § OC = R cos α q cos α − sin δ ∴ sin σ = OCR cos δ = p cos α − sin δ cos α cos δ ∴ cos σ = cos α cos δ − cos α + sin δ cos α cos δ = sin δ sin α cos α cos δ = tan δ tan α ∴ cos σ = tan α tan δ = − tan λ tan δ ∴ σ = arccos( − tan λ tan δ ) (10)The identity (10) also applies in case (ii) of the southern hemisphere, at latitude − δ , since from Figure 7(ii) therequired sunrise solar hour angle is 180 ◦ − σ , and arccos as a function from [ − ,
1] to [0 , ◦ ] satisfies the conditionarccos( − x ) = 180 ◦ − arccos( x ). The equation (10) is called the ‘Sunrise Equation’ [23], [24]. The argument to the arccoswill go out of range for latitudes beyond the effective polar circles (ie the polar circles corresponding to the effective axialtilt), for at these latitudes and times there is no sunrise nor sunset.We have shown (10) holds for any day in Q
1, using Figures 1 and 2 with an effective axial tilt α ∈ [0 , α ]. For a dayin Q − ψ ), by symmetry the situation of the planet is the same as for a day in Q ψ , thus σ isthe same. Then since cos(360 − ψ ) ◦ = cos ψ and from (8) the solar declination λ is the same, equation (10) produces thecorrect sunrise solar hour angle σ for Q Q ψ ) ◦ and with solar declination − λ , the geometric situation of the planet is identicalto a day in Q ψ with solar declination λ ( < Q
1, equation (10) gives the value of σ with the sun at its normal side of Figure 1, and then the required sunrise solar hour angle for the sun at the opposite side is (180 − σ ) ◦ , since the dayand night arcs on the latitude circles will be reversed. Thus the required sunrise solar hour angle for the day in Q − σ ) ◦ . But from (10), for any latitude δ , we have :180 ◦ − σ = arccos(tan( − λ ) tan δ )12nd so (10) gives the correct sunrise solar hour angle for the day in Q λ .For a day in Q − ψ ) ◦ , by symmetry the situation of the planet is the same as for a day in Q ψ ) ◦ , and so σ is the same. Then since cos(180 − ψ ) ◦ = cos(180 + ψ ) ◦ and the solar declination λ from (8) is thesame, the equation (10) produces the correct sunrise solar hour angle σ for Q σ ∈ [0 , ◦ ]. Writing tan λ in terms ofthe orbital angle ψ , (8) and (10) give :tan λ = ± sin λ √ (1 − sin λ ) = − sin α cos ψ √ (1 − sin α cos ψ ) (noting declination is -ve in Q Q σ = arccos (cid:18) tan δ · sin α cos ψ √ (1 − sin α cos ψ ) (cid:19) (11)To convert the total solar hour angle 2 σ of daylight to the daylight duration D we simply multiply the rotationalperiod T of the planet by the fraction 2 σ/ ◦ = σ/ ◦ : D = T σ ◦ (12)Although we are using the model of a stationary day, we could use as T the mean solar day [25], [26], of 24 hoursduration in the case of Earth, rather than the sidereal day [27] length of 23.93447 hours, to compensate for the progrademotion of the Earth in its orbit during the day, which means it takes slightly longer than a sidereal day for the sun toreturn to the local meridian position at solar hour angle 0 ◦ . A planet with a retrograde rotation would have a mean solarday shorter than the sidereal day.In the case of Earth two sources of inaccuracy in equation (12) cause underestimation of the daylight duration D : • The non-zero angular diameter [14] of the solar disc, approximately 0 . ◦ , and the convention of measuring sunriseand sunset wrt the upper edge of the solar disc crossing the horizon means official sunrise time comes a few minutesbefore and official sunset time a few minutes after the times predicted by the simple geometric model. • The effect of atmospheric refraction [28] of light from an object in space arriving at an observer on Earth causesthe altitude of the object appear slightly higher than it really is. This is because entering the slower medium ofthe atmosphere from the vacuum causes the light ray to bend towards the normal between the two media, ie theangle of incidence reduces, and this implies a downwards bending, thus the eye sees the light source as higher up,similarly to a stick under water appearing higher up. This means the Sun can be seen shortly before sunrise andshortly after sunset, whilst its upper edge is still geometrically below the horizon, again lengthening the day.[3] provides details on adjustments which can be applied to (12) to compensate for these two factors.Note: the two sunrise angles σ and θ should not be confused — σ is the angle of the sunrise wrt solar noon, whilst θ is the angle of the sunrise wrt due east on the horizon — the former is an angle in the latitude plane and the latteris an angle in the horizon plane. Both are functions of time of year and latitude and from the equations (9) and (10)cos σ/ sin θ = − sin δ/ cos λ , which varies throughout the year. In this section we relate the clock time CT (or civil time) to the solar time for the Earth . The solar time is based on thedaily passage of the Sun, with a day defined as the time between successive solar noons ( § − MST varies within a range of about -14/+16 minutes according to the ‘equation of time’ [30].Clock time in the various time zones is based on the UTC time standard, which is within 0.9 seconds of anotherstandard called UT1 [31] which measures MST at longitude 0 ◦ (UT1 replaces the older GMT standard). Thus unlessgreater than this level of accuracy is required UTC can be taken as the MST at longitude 0 ◦ , [32]. Each time zone is then We shall assume a 24 hour clock. , and has a defining meridian [32] located at longitude thecorresponding integral multiple of 15 ◦ . The time zone time, ie the clock time CT, or civil time, shared by each locationin the time zone, is defined as the MST at the longitude of the defining meridian for the zone. To determine the MSTfor any other location in the time zone we need to adjust by 4 minutes for every degree of longitude difference ∆ L ◦ fromthe longitude of the defining meridian for the zone, with the easterly direction being positive [32], ie.MST = CT + 4 minutes × ∆ L ◦ ∴ MST = CT + ∆ L ◦
15 (in decimal hours) (13)Daylight Saving Time (DST) is an adjustment to civil clock time of +1 hour for the summer months (eg from Marchto October in UK), but here we will define the civil clock time CT to be exclusive of any DST changes.When requiring the time of day in formulae (eg the solar hour angle parameters σ in the sunrise equation in § τ in the altitude and azimuth formulae in § §
8) the AST is convenient, since geometrically it correspondsstraightforwardly with the solar hour angle parameter, with 24 hours of AST being equivalent to a solar hour angle 360 ◦ .AST of 12 noon is solar noon, and AST of 12 midnight is solar midnight. An ‘hour’ of AST would be a 24 th part of thesolar day and would vary from day to day, within a range of approximately 3599 to 3601 atomic seconds. To obtain thetime of day solar hour angle τ from AST :(AST −
12) = τ ·
24 = τ ∴ τ = 15 (AST − , AST in decimal hours, τ in degrees (14)To obtain τ from the clock time CT we can use (13) and (14) together with the equation of time (EOT) which expressesAST in terms of MST, [30], [32], [33], and which is a continuous function of time throughout the year, repeating itselfover a yearly cycle. We can approximate EOT as a constant for each day of the year d so that EOT( d ) equals a constantdifference AST − MST throughout the day d . EOT( d ) can thus produce AST from a known MST, and vice-versa. Thedifference EOT( d ) has two components : EC( d ) the difference due to the eccentricity of the Earth’s orbit, and OB( d )the difference due to the obliquity of the ecliptic plane. EC( d ) is an approximate sine wave of period one year with zerosat perihelion and aphelion, which are respectively about 2 weeks after the winter and summer solstices (the exact dayvarying with the year), [34]. OB( d ) is an approximate sine wave of period a half year with zeros at the equinoxes andsolstices. Using respective amplitudes of 7.66 mins and 9.87 mins for these from [30], and taking a perihelion/aphelionoffset from the winter/summer solstice of 14 days (the 2020 figure — other years will require a different value [34], [35]),and refering to the graphs in [30], [33], approximate formulae for the discretized EC( d ) and OB( d ) in decimal hours forthe day which is at offset d from the winter solstice are :EC( d ) = − . (cid:18) d − · ◦ (cid:19) OB( d ) = − . (cid:18) d · ◦ (cid:19) and EOT( d ) = EC( d ) + OB( d ) (in decimal hours) (15)Using the relation EOT( d ) = AST − MST and equations (13) and (14) we then have : τ = 15 (MST + EOT( d ) − L ◦
15 + EOT( d ) − ⇒ τ = 15 CT − ◦ + ∆ L ◦ + 15 EOT( d ) (CT, EOT in decimal hours, τ in degrees) (16)This is an approximate relation from clock time CT (ex-DST) to solar hour angle τ for the day d , and is applicable toall time zones. Values for the approximate EOT( d ) can be calculated using the above sine wave formulae or obtained fromonline calculators such as [36] and [37]. Note the approximation is based on a discretized form of the EOT — however itis sufficiently accurate to demonstrate the effect of the analemma and to generate approximate analemma graphs ( § In a few time zones, eg in India, a non-integral number is used. Sunrise, Sunset, and Solar Noon Times
The solar noon local clock time CT for Earth is derived from equation (16) by putting τ = 0 :CT = 12 − ∆ L ◦ − EOT( d ) (17)If σ ∈ [0 , ◦ ] is the sunrise/sunset solar hour angle of § τ = − σ at sunrise clock time CT R ,and τ = σ at sunset clock time CT S . Then from (16) :CT R = CT − σ
15 (18)CT S = CT + σ
15 (19)with σ given by equation (10) or (11). Example
Consider the location Madrid on the date 15 th May 2019. Longitude is − . ◦ and latitude δ = 40 . ◦ . Thetime zone is UTC+1 (ex-DST) with defining meridian longitude at +15 ◦ . Thus ∆ L ◦ = − . ◦ . Day offset from wintersolstice 21 st Dec 2018 is d = 10 + 31 + 28 + 31 + 30 + 15 = 145. Assuming a circular orbit/constant speed model tocalculate an approximation to ψ , we have ψ = ( d/ · ◦ = 143 . ◦ . Using EOT formula (15), EOT( d ) = 0 . = 12 + 18 . − . = 14:11:00 (DST) = 13:11:00 (ex-DST)and from equation (11) : σ = arccos (cid:18) tan 40 . · sin 23 . · cos 143 . √ (1 − sin . · cos . (cid:19) = 106 . ◦ ⇒ CT R = 13 . − . /
15 = 6.083 hours = 06:04:59 (ex-DST)and CT S = 13 .
19 + 106 . /
15 = 20.30 hours = 20:18:00 (ex-DST)Actual CT R = 06:58:00 (DST) = 05:58:00 (ex-DST)Actual CT S = 21:23:00 (DST) = 20:23:00 (ex-DST)where the actual sunrise, sunset, and solar noon times are obtained from [38]. The actual sunrise comes about 7 minutesbefore the time predicted by the model, and the actual sunset comes about 5 minutes after the time predicted by themodel — this correlates with the two sources of discrepancy described in §
4. The sunrise/sunset times given by [6] takeinto account the effect of refraction [39]. The solar noon prediction by contrast was more accurate and is not impactedby the above two factors.
The solar altitude angle µ [40] is the angle the line from planet to sun makes with an observer’s local horizon plane. Itis well-defined at all points on the planet at all times and is in the range [ − ◦ , ◦ ], with the sun visible for µ ≥ µ <
0, and with µ = 0 at sunrise and sunset. The formula derived here gives µ as a function of latitude δ ,the day of the year (via λ , as described in § τ the time of day expressed as the solar hour angle, § τ is negativebefore solar noon and positive after.Assume we are at a day in Q α = − λ , where λ is the solar declination ( §
3) — ie. inFigure 1, α is replaced by α — and consider a local horizontal coordinate system [41] at the point A in Figure 1 atlatitude δ with orthogonal right-handed basis vectors u (east), v (north), w (upwardly vertical) given (in terms of the xyz -coordinate system of Figure 1) by u = (0 , , v = (cos β, sin β, w = ( − sin β, cos β, A could be any pointon the half circumference of the planet including the poles N , S and β ∈ [ − α, ◦ − α ]. The other half circumference canbe obtained by rotating about N S by 180 ◦ — this ensures v always points north (see § u ′ (east), v ′ (north), w ′ (upwardly vertical) for thepoint on latitude δ of Figure 1 which is at solar hour angle τ wrt solar noon can be obtained by rotating u , v , w by angle τ about the polar axis N S according to the RH rule. Then at this angle τ , if s = − i is the direction from the planet to15he Sun, making altitude angle µ ∈ [ − ◦ , ◦ ] with the u ′ v ′ -plane, then : µ ≥ ⇒ s · w ′ = cos(90 ◦ − µ ) = sin µ and µ < ⇒ s · w ′ = cos(90 + | µ | ) = − sin | µ | = sin µ ∴ ∀ µ ∈ [ − ◦ , ◦ ] , sin µ = s · w ′ = − i · w ′ (20)Since in Figure 1 the N S rotation axis a = (sin α , cos α , w ′ = (cos τ ) w + (sin τ )( a × w ) + (1 − cos τ )( a · w ) aa × w = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i j k sin α cos α − sin β cos β (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (0 , , sin α cos β + cos α sin β )= (0 , , sin( α + β )) = (0 , , cos δ ) , and a · w = cos( α + β ) = sin δ ∴ w ′ = cos τ ( − sin β, cos β,
0) + sin τ (0 , , cos δ ) + (1 − cos τ ) sin δ (sin α , cos α , ∴ sin µ = − w ′ · i = cos τ sin β − (1 − cos τ ) sin δ sin α Since δ = 90 − ( α + β ), sin β = cos( α + δ ) ∴ sin µ = cos τ (cos α cos δ − sin α sin δ ) − (1 − cos τ ) sin δ sin α ∴ sin µ = cos τ cos α cos δ − sin δ sin α ∴ at any point in Q , sin µ = cos τ cos δ cos λ + sin δ sin λ, µ ∈ [ − ◦ , ◦ ] (21)Note the above derivation is valid for any latitude on the planet, including the poles which are stationary throughoutthe day with δ = ± ◦ and the formula (21) reducing to µ = λ (for N) and µ = − λ (for S) which is clear intuitively asthe horizon plane at a pole is parallel to the equatorial plane.The above derivation was based on Figure 1 which covers case Q α in place of α in thediagram. In Q − ψ ) ◦ , the situation of the planet wrt the sun is as Q ψ so µ must be the same. Since λ given by (8) is also the same, the angle µ ∈ [ − ◦ , ◦ ] given by (21) is then correct for Q
4, for any latitude δ and time τ .In Q ψ ) ◦ is as Q ψ , except with the sun at the exact oppositeside. Considering this symmetric Q usual side , as depicted in Figure 5, the solar hour angleis τ + 180, where τ is the solar hour angle for the Q N S axis via the RH rule — thus to go from a sunat one side to the complete opposite side add 180 ◦ to the solar hour angle), and the solar declination is − λ where λ isthe solar declination for the Q Q µ = cos( τ + 180) cos δ cos( − λ ) + sin δ sin( − λ ) , (22)and hence sin( − µ ) = cos τ cos δ cos λ + sin δ sin λ (23)where µ is the solar altitude for the Q − µ is the required solar altitude for the point in Q
3. Thus equation (21)gives the correct solar altitude for the Q Q − ψ ) ◦ is as Q ψ ) ◦ so µ must be the same. Sincecos(180 − ψ ) ◦ = cos(180 + ψ ) ◦ and λ given by (8) is the same, the angle µ ∈ [ − ◦ , ◦ ] given by (21) is then correct for Q
2, for any latitude δ and time τ .Thus equation (21) gives the solar altitude µ ∈ [ − ◦ , ◦ ] at all latitudes δ (including the poles), at all times of theyear, and all times τ of the solar day. The solar declination parameter λ depends on the day of the year as described in §
3. Equation (21) also gives an alternate route to the Sunrise Equation (10) by setting altitude µ = 0 in (21) and thensolving for solar hour angle τ . From equation (21), since the rhs is maximal at τ = 0 the sun can only be directly overhead when τ = 0, ie at solar noon.The condition for that to happen for a day is then : sin 90 ◦ = cos δ cos λ + sin δ sin λ ⇒ δ − λ ) ∴ as δ − λ ∈ ( − ◦ , ◦ ) , δ − λ = 0or sin δ = sin λ (as δ, λ ∈ [ − ◦ , ◦ ]) ∴ from (7) sin δ = − sin α cos ψ (24)16utside the tropics, where | δ | > α , (24) has no solution in ψ . On the Tropic of Capricorn at δ = − α , there is a singlesolution ψ = 0, ie on the winter solstice. On the Tropic of Cancer at δ = α , there is a single solution ψ = 180 ◦ , ie thesummer solstice. Inside the tropical region : | δ | < α ⇒ ≤ (cid:12)(cid:12)(cid:12)(cid:12) sin δ sin α (cid:12)(cid:12)(cid:12)(cid:12) < ∴ there are two distinct solutions to (24) in ψ ∈ [0 , ◦ ), ie on two days of the year the sun is directly overhead atsolar noon. In the case of the equator, these solutions are ψ = 90 ◦ and ψ = 270 ◦ , ie the two equinoxes. Let the azimuthal angle in the horizon plane be φ ∈ [0 ◦ , ◦ ) using the convention of measuring φ clockwise from duenorth [42]. Continuing from the main Q § p be the component of s parallel to the horizon plane,ie parallel to the plane of u ′ and v ′ . Assume p = , ie | µ | < ◦ , otherwise φ is undefined. Then as ( s · w ′ ) w ′ is theperpendicular component of s : p = s − ( s · w ′ ) w ′ = − i + ( i · w ′ ) w ′ = − i − w ′ sin µ , (from (20)) ∴ | p | = p · p = 1 + sin µ + 2(sin µ ) w ′ · i ∴ | p | = cos µ (since from (20) sin µ = − i · w ′ ) (25)Now as φ is defined clockwise from north (ie from v ′ ) in the u ′ v ′ -plane, φ satisfies : | p | sin φ = p · u ′ , and | p | cos φ = p · v ′ , ∴ | p | sin φ = ( − i − w ′ sin µ ) · u ′ = − u ′ · i (since basis vectors u ′ , v ′ , w ′ are orthogonal) , (26)and | p | cos φ = ( − i − w ′ sin µ ) · v ′ = − v ′ · i (27)But from the Rodrigues Rotation Formula : u ′ = (cos τ ) u + (sin τ )( a × u ) + (1 − cos τ )( a · u ) aa × u = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i j k sin α cos α
00 0 1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cos α , − sin α , , and a · u = 0 ∴ u ′ = cos τ (0 , ,
1) + sin τ (cos α , − sin α , ∴ u ′ · i = sin τ cos α = sin τ cos λ (28)and v ′ = (cos τ ) v + (sin τ )( a × v ) + (1 − cos τ )( a · v ) aa × v = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i j k sin α cos α β sin β (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (0 , , sin α sin β − cos α cos β )= (0 , , − cos( α + β )) = (0 , , − sin δ ) , and a · v = sin( α + β ) = cos δ ∴ v ′ = cos τ (cos β, sin β,
0) + sin τ (0 , , − sin δ ) + (1 − cos τ ) cos δ (sin α , cos α , ∴ v ′ · i = cos τ cos β + (1 − cos τ ) cos δ sin α Since δ = 90 − ( α + β ), cos β = sin( α + δ ) ∴ v ′ · i = cos τ sin( α + δ ) + (1 − cos τ ) cos δ sin α = cos τ (sin α cos δ + cos α sin δ ) + (cos τ −
1) cos δ sin α = cos τ sin δ cos α + cos δ sin α = cos τ sin δ cos λ − cos δ sin λ (29)Thus from (25), (26), (27), (28), and (29)sin φ = − sin τ cos λ cos µ (30)and cos φ = cos δ sin λ − cos τ sin δ cos λ cos µ (31)17quation (31) can also be written without a τ term :cos φ = (1 − sin δ ) sin λ − cos τ sin δ cos λ cos δ cos µ cos δ = sin λ − sin δ (cos τ cos δ cos λ + sin δ sin λ )cos µ cos δ = sin λ − sin δ sin µ cos µ cos δ , (from (21)) (32)Then equations (30) and (31) (or equations (30) and (32)) uniquely define φ ∈ [0 ◦ , ◦ ), for the case of Q
1. For thecase of Q
4, at an angle (360 − ψ ) ◦ the situation of the planet is as at angle ψ in Q φ applies. Since λ given by (8) is the same at these two angles it follows that the rhs’s of (30) and (31) produce the same results, andhence give the correct φ for Q Q ψ ) ◦ , the same process of reflection as in § Q ψ can be performed, with theresult that (30) and (31) produce negative of the Q φ . Thus the required angle for Q φ + 180 ◦ ) mod 360 ◦ , is produced by (30) and (31). (Note the cos µ term does not change sign in the reflection).In Q − ψ ) ◦ is as Q ψ ) ◦ so φ must be the same. Sincecos(180 − ψ ) ◦ = cos(180 + ψ ) ◦ and λ given by (8) is the same, the angle φ ∈ [0 ◦ , ◦ ) given by (30) and (31) is thencorrect for Q
2, for any latitude δ and time τ .Thus equations (30) and (31) (or equations (30) and (32)) uniquely define the solar azimuth angle φ ∈ [0 ◦ , ◦ ) atall latitudes δ at all times of the year, and all times τ of the solar day, provided an azimuth is defined ie the sun is not atthe zenith or nadir position. Since at the poles cardinal directions of north, south, east, and west are not well-defined ameridian longitude (a half great circle) has to be selected, wrt which the cardinal directions can be defined. Then solarnoon time is well-defined as the time when the sun crosses the selected meridian, even though on a stationary day inthe orbit (in the simplified model) the sun’s altitude is constant (possibly negative) all throughout the day. North at apole would then be defined as the continuation of northwards on the selected longitude line, which would be towards thenorth pole and away from the south pole. With this convention equations (30) and (31) produce the correct azimuthangle using the north clockwise convention, with τ measured wrt to the selected meridian longitude :North pole ⇒ µ = λ ⇒ sin φ = − sin τ and cos φ = − cos τ ∴ φ = τ + 180 ◦ (mod 360 ◦ )which is the correct angle wrt due north as τ is wrt due south, andSouth pole ⇒ µ = − λ ⇒ sin φ = − sin τ and cos φ = cos τ ∴ φ = 360 ◦ − τ (mod 360 ◦ )ie φ = − τ (mod 360 ◦ )which is the correct angle wrt due north since at the south pole north points into the selected meridian half plane, and apositive rotation τ from solar noon will produce a negative solar azimuth φ according to the north clockwise convention.The above process will work for any choice of meridian longitude for a pole.The azimuth cosine formula (32) provides an alternate route to the sunrise direction formula (9) since putting µ = 0in (32) and remembering θ is measured anti-clockwise from due east we obtain : θ + φ = 90 ◦ ⇒ sin θ = cos φ = sin λ cos δ This completely characterizes θ since from § θ ∈ [ − ◦ , ◦ ]. At the north pole every direction points to the south and at the south pole every direction points to the north. xample Continuing with the example of Madrid on 15 th May 2019 in §
6, compare the above altitude and azimuthformulae with the actual solar altitude of 50 ◦ and azimuth of 249 ◦ at DST 16:47 from [38]. We have :CT = 15:47 = 15.78 hours τ = 15 · . − − .
72 + 15 · . . ◦ , ∴ cos τ = 0 . , sin τ = 0 . λ = − sin 23 . · cos 143 . . , cos λ = 0 . ∴ sin µ = 0 . · (cos 40 . · . . · . . ∴ µ = 50 . ◦ , cos µ = 0 . . And sin φ = − . · . / . − . φ = (0 . − (sin 40 . · . / (0 . · cos 40 .
42) = − . , from (32) ⇒ φ = 248 . ◦ Thus the formulae predictions are accurate to within 0 . ◦ when compared with [38]. The solar noon altitude [24] µ for a location is by definition the solar altitude when the solar hour angle τ for the locationis 0 ◦ . From (21) this is the maximum value of µ throughout the day, ie the highest position of the sun. On a day ofperpetual darkness this will be negative, eg as in Figure 1 at a point above the arctic circle on the winter solstice (wherecase 1 below applies). From Figure 1 with the effective axial tilt α , covering cases Q Q
4, it is clear the azimuth φ of the sun at solar noon for any day and latitude is either :1. southerly, ie φ = 180 ◦ , for latitudes above the effective Tropic of Capricorn (at latitude − α )2. northerly, ie φ = 0 ◦ , for latitudes below the effective Tropic of Capricorn3. undefined on the effective Tropic of Capricorn, as the sun is then directly overhead at solar noon (ie µ = 90 ◦ )The effective Tropic of Capricorn is the southern tropic corresponding with the effective axial tilt α . The effectivetropics converge to the equator at the spring equinox where α becomes zero. Note cases (1) and (2) do not correspondwith the northern and southern hemispheres, but with the two sides of the effective tropic which changes throughout theyear. For example in Figure 1 a southern latitude between − α and 0 ◦ will see the noon sun in the southern half of thesky, as a northern latitude will. The dividing line is the effective tropic, which at the spring equinox becomes the equator.For Q Q α . Forpractical purposes which of the above three cases applies would be determinable whenever the sun is visible at solar noon.Considering case 1, entering τ = 0 , φ = 180 , µ = µ in equations (21) and (31) we obtain :sin µ = cos δ cos λ + sin δ sin λ cos µ = sin δ cos λ − cos δ sin λ ie. cos(90 − µ ) = cos( δ − λ )and sin(90 − µ ) = sin( δ − λ ) ∴ (90 − µ ) − ( δ − λ ) = integral multiple of 360But then, defining the solar zenith angle [40] ν at solar noon to be the complement of angle µ : µ ∈ [ − , , δ ∈ ( − , , λ ∈ [ − α, α ] , α ∈ [0 ,
90) (33) ⇒ (90 − µ ) − ( δ − λ ) ∈ ( − , ⇒ (90 − µ ) − ( δ − λ ) = 0 ⇒ ν = δ − λ (34)Considering case 2, entering τ = 0 , φ = 0 , µ = µ in equations (21) and (31) we obtain :sin µ = cos δ cos λ + sin δ sin λ cos µ = cos δ sin λ − sin δ cos λ ie. cos(90 − µ ) = cos( δ − λ )and sin(90 − µ ) = − sin( δ − λ )Thus in the xy -plane on the unit circle the angles 90 − µ and δ − λ are reflections of one another in the x -axis.19 − µ = 360 − ( δ − λ ) (mod 360) ∴
270 + µ − δ + λ = integral multiple of 360but from (33), µ − δ + λ ∈ ( − , ∴
270 + µ − δ + λ ∈ (0 , ∴
270 + µ − δ + λ = 360 ∴ ν = λ − δ (35)In case 3, entering τ = 0 , µ = 90 in equation (21) we obtain :sin 90 = cos δ cos λ + sin δ sin λ ie. 1 = cos( δ − λ ) ∴ δ − λ = integral multiple of 360 ⇒ from (33), δ − λ = 0 ⇒ λ = δ (36)The three solar noon equations (34), (35), and (36) give a relationship between latitude λ and solar declination δ depending on which of the above cases 1–3 applies. These equations will still apply even if µ <
0, ie a day ofperpetual darkness, however in case 3 the one possibility not considered, ie µ = −
90, with the sun at the nadirposition at solar noon, is not possible, as we would expect intuitively, but formally by (21) such a position would imply :sin( −
90) = − δ − λ ) ⇒ δ − λ −
180 = multiple of 360. But since from (33), δ − λ − ∈ ( − ,
10 Analemma
One form of analemma [43] is a graph of the Sun’s position in the sky as seen from a fixed location on Earth, as altitudeplotted against azimuth, at the same local clock time CT (ex-DST) each day for a period of a year. With equal horizontaland vertical scales this will have the same shape as a trace of the Sun’s actual positions over these times, for example ascaptured in analemma photographs. The irregular figure 8 shape of the analemma (Figure 8) is a result of the deviationthroughout the year of the AST from the MST, as given by the equation of time ( § τ for the day d , which is then entered in equations (21), (30), and (31)to produce the altitude µ and azimuth φ . The parameter λ also varies daily, according to equation (8). If the solar daylength were constant throughout the year so that MST = AST, and EOT ( d ) = 0 ∀ d , then τ would be constant for clocktime CT each day, and the analemma would be a simple curve without any loops because the solar declination parameter λ retraces its values in [0 , ◦ ] in the return journey from 180 ◦ to 360 ◦ due to the cos ψ term in equation (8).The Perl script calc-analemma.pl in Appendix C.3 calculates the analemma for a given location and time usingthe above formulae, producing a table of altitude and azimuth values over a year that can be directly input into theTIKZ/pgfplots package to produce a graph (eg by compiling the file analemma-graph.tex ). The script uses the circularorbit/constant speed model to calculate an approximate orbital angle ψ for each day of the year. To use the script specifythe desired local clock time CT (ex-DST), the latitude δ , and the difference ∆ L ◦ in degrees between the longitude of thelocation and the longitude of the defining meridian for the time zone .Uncommenting a line in the script shows the effect of a constant solar day length, and other lines can be uncommentedto show the separate effects of the eccentricity and obliquity components of the EOT.Note that an analemma in the tropics will become ill-defined towards solar noon because at such locations at certaintimes of the year the Sun is directly overhead at solar noon ( § /
0. Note because of equation (21),which implies µ is maximal at τ = 0, the Sun directly overhead can only occur at solar noon (and likewise the nadirposition only at solar midnight). Note the longitude itself is not sufficient information as a single longitude can span more than one time zone, eg the UK and Spain spantime zones UTC and UTC+1. § . ◦ .In Figure 8 below two analemmas generated by the script are shown :(i) Athens, Greece : latitude 37 . ◦ N, longitude = 23 . ◦ E, UTC+2, ∆ L ◦ = − . ◦ , CT = 16:00Script command line : perl calc-analemma.pl 37.98 16 -6.27 > analemma-athens.dat Photo : (ii) Kumagaya, Japan : latitude 36 . ◦ N, longitude = 139 . ◦ E, UTC+9, ∆ L ◦ = 4 . ◦ , CT = 7:00Script command line : perl calc-analemma.pl 36.15 7 4.38 > analemma-kumagaya.dat Photo : https://earthsky.org/todays-image/todays-image-analemma-2013
230 240 250 260 27010203040 (i) Athens, Greece, 16:00 hoursAzimuth φ ( ◦ ) A l t i t ud e µ ( ◦ )
80 90 100 110 1200102030 (ii) Kumagaya, Japan, 07:00 hoursAzimuth φ ( ◦ ) A l t i t ud e µ ( ◦ ) Figure 8: Analemmas Generated by Formulae (16), (21), (30), and (31)For a planet other than the Earth, with an axial tilt α ∈ [0 , ◦ ), provided the original assumptions 1–6 in § d to orbital angle ψ would be required.With the required changes the script should then be able to generate analemmas for the planet.
11 Acknowledgements
Dedicated to my nephew Rory, a great sportsman, who passed away all too soon in 2015, who loved astronomy and whoalways encouraged my studies. God bless Rory, your light shines each day. And with much appreciation also to myparents, my teachers, and my brother and sister and their families. b Whatsoever thy hand findeth to do, do it with thy might c ppendix A Comparison With Actual Sunrise DataComparison With Actual Sunrise Data A.1 Sunrise Predictions
The actual sunrise directions θ for locations at 8 different latitudes (4 northern hemisphere, 3 southern hemisphere, 1 equa-torial) at 10 day intervals over the year 2018-2019 are shown in Table A.1 below (source : ,[6]). gives the sunrise and sunset directions to the nearest degree. Where in the figures the sunrise angle north of east differs from the sunset angle north of west on a day the average of these two isquoted in the table below. The day offset n is from the 2018 winter solstice date of 21/12/2018.Day Offset Abu Dhabi Edinburgh Melbourne Milan Quito Reykjavik Rio Stanley n . ◦ . ◦ − . ◦ . ◦ − . ◦ . ◦ − . ◦ − . ◦ θ for 2018-2019The charts below compare the actual sunrise directions shown in red with the curve of sunrise directions predictedby the sunrise direction formula (with the currently accepted Earth axial tilt α = 23 . ◦ ) shown in blue. The circularorbit/constant speed model is used with the orbital angle ψ in the formula calculated as dN (360 ◦ ) for day offset d ∈ [0 , N − N = 365 is the number of days in the year. The average error over the 8 chartsis 1 . ◦ . When the formula’s own axial tilt estimate of 23 . ◦ (Appendix B) is used the average error is 1 . ◦ . − − . ◦ Latitude δ = 55 . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inEdinburgh, 2018-19
Actual sunriseComputed sunrise0 50 100 150 200 250 300 350 400 − − − . ◦ Latitude δ = 45 . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inMilan, 2018-19
Actual sunriseComputed sunrise
50 100 150 200 250 300 350 400 − − − . ◦ Latitude δ = 64 . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inReykjavik, 2018-19
Actual sunriseComputed sunrise0 50 100 150 200 250 300 350 400 − − . ◦ Latitude δ = 24 . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inAbu Dhabi, 2018-19
Actual sunriseComputed sunrise
50 100 150 200 250 300 350 400 − − − . ◦ Latitude δ = − . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inMelbourne, 2018-19
Actual sunriseComputed sunrise0 50 100 150 200 250 300 350 400 − − − − . ◦ Latitude δ = − . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inStanley, Falkland Islands, 2018-19
Actual sunriseComputed sunrise
50 100 150 200 250 300 350 400 − − . ◦ Latitude δ = − . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inRio De Janeiro, 2018-19
Actual sunriseComputed sunrise0 50 100 150 200 250 300 350 400 − − . ◦ Latitude δ = − . ◦ Day Offset ( n ) from Winter Solstice D e g r ee s ( θ ) N o rt h o f E a s t Actual/Computed Sunrise inQuito, Ecuador, 2018-19
Actual sunriseComputed sunrise .2 Perpetual Day/Night Predictions Using the same circular orbit/constant speed model as in Appendix A.1, consider the location of Jan Mayen, Norway atlatitude 71 ◦ north. The argument to the arcsin goes out of range when either:cos ψ > cos δ sin α = cos 71 ◦ sin 23 . ◦ ≃ . ≃ cos 35 ◦ , or cos ψ < − cos δ sin α ≃ − . . These correspond with ψ ranges of [0 , , [145 , , [325 , ◦ south, the estimated dates for 2018-2019 are as for Jan Mayen and theactual dates are:Actual date 31/1/2019 19/5/2019 26/7/2019 14/11/2019Estimated date 25/1/2019 17/5/2019 27/7/2019 15/11/2019Error (days) 6 2 1 1Average error 2.5 daysFor Longyearbyen, Svalbard, Norway at latitude 78 ◦ north, the arcsin goes out of range when cos ψ > . ≃ cos 58 . ◦ or cos ψ < − . ψ ranges [0 , . , [121 . , . , [301 . , ◦ south, the estimated dates for 2018-2019 are as for Longyearbyenand the actual dates are :Actual date 20/2/2019 25/4/2019 19/8/2019 24/10/2019Estimated date 18/2/2019 23/4/2019 20/8/2019 23/10/2019Error (days) 2 2 1 1Average error 1.5 daysThe average errors overall is 3.1 days which though not as good as the sunrise predictions of § A.1 or the axial tiltprediction of Appendix B are still reasonable approximations for the simplified model.27 ppendix B Axial Tilt EstimatesAxial Tilt Estimates
The axial tilt estimates of the sunrise direction formula obtained by applying equation (5) to the sunrise data of Table A.1,using the same circular orbit/constant speed model as in Appendix A.1, are shown in Table B.1 below. The estimatescan become abnormally inaccurate when observed values of sin θ and cos ψ are close to zero near the equinoxes ( n ≃ n ≃
270 for autumn), for then small absolute changes in these two terms can result in large % changes,causing large % change in sin α on lhs of (5). The rhs of (5) is thus unstable near the equinoxes (and at the equinoxes ithas form 0 / δ term does not cause instability so long as we are not too close to the poles. A good directionapproximation can have a large % error when we are near zero - eg approximating direction of 0 . ◦ by 0 . ◦ has 100%error but is still a good direction approximation. Thus the particularly bad 5 estimates near the spring equinox in thetable which were negative are excluded in calculating the averages — though any other bad estimates are left in. Overthe remaining 299 data points the overall average estimate for the Earth’s axial tilt is 23 . ◦ , which is within 0 . ◦ of thecurrently accepted value of 23 . ◦ .Day Offset Abu Dhabi Edinburgh Melbourne Milan Quito Reykjavik Rio Stanley Overall n . ◦ . ◦ − . ◦ . ◦ − . ◦ . ◦ − . ◦ − . ◦ Average0 22.63 22.89 24.01 23.09 23.00 22.64 23.82 24.0010 22.99 22.80 24.39 22.81 23.37 22.54 24.20 24.3820 23.16 22.48 24.01 22.56 23.45 22.58 23.46 23.9230 23.09 22.24 24.36 22.25 23.16 22.24 23.38 24.1440 22.57 21.93 24.08 22.57 23.59 22.20 23.47 24.0550 22.64 21.70 23.86 21.59 23.39 21.66 22.93 23.7060 21.68 20.84 23.52 20.75 21.86 21.08 22.90 23.8270 20.74 19.80 23.73 19.90 21.40 19.63 22.38 22.9580 19.27 16.19 23.17 16.61 21.25 16.03 19.51 23.1190 -21.67 -42.94 -18.69 -34.67 0.00 -62.12 0.00 14.56100 31.92 33.47 24.40 34.72 27.70 35.40 25.36 25.58110 26.68 28.15 22.93 27.37 26.03 28.78 23.84 23.98120 25.54 26.21 23.73 26.36 25.96 26.88 23.79 23.79130 25.48 25.58 23.24 25.13 24.73 25.72 24.23 23.05140 24.73 24.88 23.44 24.87 24.54 25.15 23.77 23.02150 23.73 24.57 23.20 24.44 23.80 24.76 22.93 22.81160 23.57 24.36 22.79 24.36 23.86 24.45 22.88 22.88170 24.11 24.35 23.08 24.31 23.58 24.35 23.49 23.00180 23.54 24.20 22.54 24.37 23.02 24.21 22.94 22.98190 23.73 24.18 22.72 23.93 23.20 24.23 23.12 22.64200 23.76 24.26 22.85 24.25 23.10 24.34 23.10 22.71210 24.07 24.38 22.90 24.31 23.74 24.40 22.82 22.92220 24.11 24.64 22.73 24.34 24.05 24.72 23.23 22.83230 24.29 24.94 23.27 24.65 23.77 24.99 23.19 22.52240 25.42 25.54 23.38 25.91 24.20 26.01 23.96 22.73250 25.91 27.29 24.41 27.17 25.90 27.05 23.73 22.98260 28.24 29.77 24.24 29.58 26.47 29.80 24.25 24.43270 47.61 49.16 32.29 49.32 32.75 55.46 29.90 24.78280 8.51 10.48 18.72 9.84 18.97 4.06 17.42 23.75290 18.42 17.44 21.93 18.03 20.31 17.54 20.41 22.95300 21.22 20.04 22.10 19.51 20.99 19.97 22.61 23.04310 22.14 21.04 23.29 21.16 22.64 20.86 22.42 23.42320 21.86 21.58 23.34 21.57 22.69 21.54 23.48 23.61330 22.20 22.03 23.43 21.91 23.27 21.97 23.62 23.84340 23.04 22.29 23.67 22.69 23.22 22.21 23.34 23.94350 22.52 22.33 24.11 22.60 23.32 22.49 23.75 24.07360 22.72 22.54 24.11 22.54 23.09 22.51 23.92 24.09365 22.63 22.89 24.01 23.09 23.00 22.64 23.82 24.00Average ◦ ) excluding any negative estimates of axial tilt. ppendix C Perl ScriptsPerl Scripts The scripts C.1 and C.2 below use the circular orbit/constant speed model for Earth with stationary days evenly spacedaround the circle and the orbital angle ψ defined as dN (360 ◦ ) for day d ∈ [0 , N − d is the day offset from theday of the winter solstice, and N = 365 is the number of days in the year. To use the more accurate mapping from dayof the year to ψ described in [3] make the following changes to these scripts (the time of noon is used for each day, andeach day duration is taken to be a mean solar day of 24 hours). A similar change can be made to script C.3. • replace the line my $winter_solstice = DateTime->new(day => 21, month => 12, year => 2018); with my $base_date = DateTime->new(day => 1, month => 1, year => 2013); • replace the line my $day_offset = $date->delta_days($winter_solstice)->in_units(’days’); with my $t = $date->delta_days($base_date)->in_units(’days’) + 0.5; • replace the line $psi = ($day_offset / 365) * pi2;) with the formula of [3] as a function of $t , plus $pip2 (ie90 ◦ ) C.1 Calculate Sunrise Direction rintf("@fields\n");next;} Sample output :
C:\>perl calc-sunrise.pl iedinburgh.txtActual Sunrise Computed Sunrise Abs Error-------------- ---------------- ---------....-43.0 -44.42 -1.42-40.0 -41.97 -1.97-36.0 -38.16 -2.16....EdinburghAverage magnitude abs error = 2.14 degrees38 data points .2 Calculate Axial Tilt ay => $day,month => $month,year => $year); Sample output :
C:\>perl calc-axial-tilt.pl iedinburgh.txtPSI cos(PSI) Actual Sunrise Axial Tilt------ ---------- -------------- ----------....147.95 -0.85 39.0 24.57157.81 -0.93 43.0 24.36167.67 -0.98 46.0 24.35....Edinburgh37 data pointsAverage axial tilt = 23.88 .3 Calculate Analemma perihelion_offset = 14;}if ($adjust_azimuth_angle_range eq ’.’) {$adjust_azimuth_angle_range = 0;}if ($num_days eq ’.’) {$num_days = 365;}if ($day_step eq ’.’) {$day_step = 1;}if ($longitude_diff eq ’.’) {$longitude_diff = 0;}my $clock_time_ex_dst = pop @ARGV;my $latitude_degrees = pop @ARGV; EOT component amplitudes in decimal hoursmy ($EC_amplitude, $OB_amplitude) = (0.1277, 0.1645);my $EC_component = $EC_amplitude * -sin( ($day_offset - $perihelion_offset) * pi2 / 365 );my $OB_component = $OB_amplitude * -sin( $day_offset * 2 * pi2 / 365 );my $EOT_hours_adjustment = $EC_component + $OB_component;
Sample output :
C:\>perl calc-analemma.pl 37.98 16 -6.27 | perl -S tee.pl analemma.datDay Altitude Azimuth analemma-graph.tex : \documentclass[a4paper]{article}\usepackage[margin=0.3in]{geometry}\usepackage{tikz}\usepackage{pgfplots}\pgfplotsset{width=12cm,compat=1.14}\newcommand{\dg}{^{\circ}}\begin{document}\begin{tikzpicture}[node font=\Large]\definecolor{Grid Color}{HTML}{aaaaaa};\definecolor{Analemma Color}{HTML}{ff0000}\begin{axis}[clip=false,axis equal=true,xlabel={Azimuth ($\dg$)},xlabel style={anchor=north, below=8pt},ylabel={Altitude ($\dg$)},ylabel style={anchor=south, above=8pt},grid=both,grid style={Grid Color, line width=0.3pt, opacity=0.3},]\addplot [Analemma Color,only marks,mark size=1.1pt,] table[x index=2, y index=1, skip first n=0] {analemma.dat} -- cycle;\path (axis description cs:0.5, 1.1) node[] {$<$Title$>$};\end{axis}\end{tikzpicture}\end{document} ppendix D Rodrigues Rotation FormulaRodrigues Rotation Formula In the following EV denotes the set of all Euclidean vectors. Rotations about an axis follow the right hand rule convention.
Theorem.
Given a unit vector u ∈ EV then the rotation r ′ of vector r ∈ EV about the axis u is given by : r ′ = (cos θ ) r + (sin θ ) ( u × r ) + (1 − cos θ ) ( u · r ) u (37) Proof.
As shown in Figure D.1 express r as a sum of components parallel and perpendicular to u : r = r + r ⊥ , where r = ( u · r ) u , and r ⊥ = r − r = r − ( u · r ) u . Since we are using the right hand rule convention the image of r ⊥ under the rotation is given by : r ′⊥ = (cos θ ) r ⊥ + (sin θ ) ( u × r ⊥ )= (cos θ ) r ⊥ + (sin θ ) ( u × r ) , since u × r = , and clearly r ′ = r , thus : r ′ = r ′ + r ′⊥ = ( u · r ) u + (cos θ ) ( r − ( u · r ) u ) + (sin θ ) ( u × r )which readily rearranges to (37). QED PP ′ rr r ⊥ u × r ⊥ r ′⊥ θ u r ′ Figure D.1: Rotating Vector r by θ About Axis u ppendix E Latitude as Elevation of Pole StarLatitude as Elevation of Pole Star In Figure E.1, an observer at latitude δ above the equator sees an angle of elevation γ of the Pole Star above their horizon.Because the distance OP is vastly much greater than the planet radius R the lines P A and
P N are close to parallel,thus the angles of incidence δ and γ of the horizon line onto these two lines are close to equal. (And as we approach theequator the Pole Star falls to the horizon with γ < δ , δ → sin − ( ROP ) + , γ → + , and R ≪ OP ). Thus γ always wellapproximates δ as a direction. OSN P (Pole Star) A (Observer)EQUATOR R δ H O R I Z O N P L AN E δ γ Figure E.1: Latitude as Elevation of Pole Star38 eferences
All web links retrieved on 31 st August 2020.[1] Position of the Sun, https://en.wikipedia.org/wiki/Position_of_the_Sun [2] NOAA Solar Calculator, [3] Alejandro Jenkins (2013), The Sun’s position in the sky, https://arxiv.org/abs/1208.1043 [4] Equatorial seasons and climate, https://en.wikipedia.org/wiki/Equator [5] True north, https://en.wikipedia.org/wiki/True_north [6] Sun Calculator, [7] Earth, Axial Tilt and Seasons, https://en.wikipedia.org/wiki/Earth [8] A. Pannekoek (1990),
A History of Astronomy , Dover Books.[9] Kristen Lippincott et al. (1999),
The Story Of Time , Merrell Hoberton Publishers (in association with NationalMaritime Museum).[10] The Calendar, BBC Radio 4 In Our Time podcast, [11] The Measurement of Time, BBC Radio 4 In Our Time podcast, [12] How long are sunset and sunrise at the North and South Pole? [13] Sunrise and Sunset at the Poles, [14] Angular diameter, https://en.wikipedia.org/wiki/Angular_diameter [15] Twilight Definitions, https://en.wikipedia.org/wiki/Twilight [16] How long does a sunrise or sunset take? https://astronomy.stackexchange.com/questions/12824/how-long-does-a-sunrise-or-sunset-take [17] Rodrigues Rotation Formula, https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula [18] Ecliptic longitude of soltices and equinoxes, https://en.wikipedia.org/wiki/Ecliptic [19] Declination, https://en.wikipedia.org/wiki/Declination [20] Equatorial coordinate system, https://en.wikipedia.org/wiki/Equatorial_coordinate_system [21] Solar hour angle, https://en.wikipedia.org/wiki/Hour_angle [22] Solar noon https://en.wikipedia.org/wiki/Noon [23] Sunrise equation, https://en.wikipedia.org/wiki/Sunrise_equation [24] Solar Geometry, http://mypages.iit.edu/~maslanka/SolarGeo.pdf [25] Solar time, https://en.wikipedia.org/wiki/Solar_time [26] Apparent and mean solar day, https://en.wikipedia.org/wiki/Day [27] Sidereal time, https://en.wikipedia.org/wiki/Sidereal_time [28] Atmospheric refraction, https://en.wikipedia.org/wiki/Atmospheric_refraction [29] Apparent solar time, https://en.wikipedia.org/wiki/Solar_time [30] Equation of time, https://en.wikipedia.org/wiki/Equation_of_time [31] Universal Time, https://web.archive.org/web/20190821105541/https://aa.usno.navy.mil/faq/docs/UT.php [32] The Equation of Time, https://web.archive.org/web/20190820154547/http://aa.usno.navy.mil/faq/docs/eqtime.php [33] Hughes, Yallop, Hohenkerk (1989), The Equation of Time,
Monthly Notices of the Royal Astronomical Society ,vol 238, p1529-1535, http://adsabs.harvard.edu/full/1989MNRAS.238.1529H [34] Perihelion, Aphelion and the Solstices, [35] Solstices & Equinoxes, [36] Equation of Time Calculations, [37] Sundial Time Correction - Equation of Time, http://mb-soft.com/public3/equatime.html [38] Sun Calculator - Madrid, May 2019, [39] What Is Refraction of Light? [40] Solar zenith angle, https://en.wikipedia.org/wiki/Solar_zenith_angle [41] Horizontal coordinate system, https://en.wikipedia.org/wiki/Horizontal_coordinate_system [42] Solar azimuth angle, https://en.wikipedia.org/wiki/Solar_azimuth_angle [43] Analemma, https://en.wikipedia.org/wiki/Analemmahttps://en.wikipedia.org/wiki/Analemma