Sun-Earth Kinematics, the Equation of Time, Insolation and the Solar Analemma

1 Sun-Earth Kinematics, the Equation of Ti me, Insolation and the Solar Analemma Phil Lucht Rimrock Digital Technology, Salt Lake City, Utah 84103 last update: January 8, 2013 Maple code is available upon request. Comments and errata are welcome. The material in this document is copyrighted by the author. The graphics look ratty in Windows Adobe PDF viewers wh en not scaled up, but look just fine in this excellent freeware viewer: http://www.tracker-software.com/pdf-xchange-products-co mparison-chart . Overview .................................................................................................................................................. 4 Summary.................................................................................................................................................. 5 1. Geometry of th e Earth's Orbit......................................................................................................... 10 (a) The Heliocentric Picture................................................................................................... ............. 10 (b) Development of the Geocentric Picture ........................................................................................ 13 (c) Relation between the th ree angular velocities .............................................................................. .15 (d) Equation of the ellipse.................................................................................................... ............... 16 (e) The Anomalies.............................................................................................................. ................. 17 2. The Sun's Apparent Path in the Sky .......................................................................................... .....18 3. The Celestial Sphere and the Sun's Daily and Annual Paths Thereon ........................................ 20 (a) The Celestial Sphere and the Tropics of Cancer and Capricorn.................................................... 20 (b) Why is R >> R Earth for the celestial sphere? ................................................................................ 22 (c) Annual apparent motion of the sun around the earth: the ecliptic................................................. 23 (d) Precession effects .......................................................................................................................... 27 (e) Daily sun paths for various observation points............................................................................. .29 (f) The Arctic and An tarctic Circles ........................................................................................... ........ 33 4. How the Earth's Orbital Azimuthal Position ψ Varies with Time due to Eccentricity .............. 35 (a) A simple solution for ψ(t).............................................................................................................. 35 (b) A systematic solution for ψ(t) and estimate of error in using (4.17)............................................. 40 5. How the Earth's Tilt Angle θt Varies with Time ........................................................................... 44 6. Coordinates of the Sun on the Celestial Sphere ............................................................................. 47 (a) Cartesian coordinates x,y,z and velocity of the sun in Frame S .................................................... 47 (b) Spherical coordinates θ, φ of the sun in Frame S.......................................................................... 48 (c) The Equation of Time....................................................................................................... ............. 51 (d) Plots of φ (right ascension of the sun) and Δφ (equation of time) ................................................ 52 (e) Plots of θt (declination of the sun) ................................................................................................ 59 7. Coordinates of the Sun in an Arbi trary Frame S' Glued to the Earth......................................... 61 (a) Setup: θ1, φ1, φb, φb0, τ, tMe........................................................................................................... 61 (b) Obtaining a value for φb0 in terms of φG0...................................................................................... 63 (c) Cartesian coordinates x',y',z' a nd velocity of the sun in Frame S' ................................................. 64 (d) Spherical coordinates θ',φ' of the sun in Frame S'......................................................................... 67 2 8. The Equation of Time Revisited ............................................................................................. ........ 70 (a) Definition of new azimuths Φ ' and Φ's.......................................................................................... 70 (b) Local Time, Sundial Time and the Equation of Time ................................................................... 72 9. Noon and Midnight, Sunrise and Sunset....................................................................................... ..75 (a) Conditions for Noon and Midnight ........................................................................................... ....75 (b) Local time of Noon and Midnight : approximations..................................................................... 77 (c) Another approach to noon ................................................................................................... .......... 79 (d) Sunrise and Sunset, Length of Day, Noon Elevation of the Sun................................................... 80 10. Insolation ................................................................................................................ ........................ 86 (a) The t ss function and hours of daylight .......................................................................................... 86 (b) Daily Insolation versus latitude and day of year ........................................................................... 90 (c) Annual Insolation pl otted versus latitude ..................................................................................... 95 11. The Sola r Analemma ......................................................................................................... ............. 97 (a) The Sun Sphere and the Analemma ............................................................................................ ..97 (b) Qualitative appearance of the inscribed analemma in photos ....................................................... 99 (c) Calculation of the an alemma in angle space................................................................................ 101 1. Location of the actual sun and the mean fictitious sun ............................................................. 102 2. Relationship between local time LMT and master time τ : shooting times τN......................... 103 3. Location of the actual and fictitious strobed suns in terms of LMT ......................................... 104 4. View of the continuous θ',Φ' data before strobing .................................................................... 105 (d) Calculation of the analemma in film-space................................................................................. 106 1. Overview.................................................................................................................... ............... 106 2. The transformation from Frame S to Frame S'.......................................................................... 107 3. The transformation from Frame S' to Frame S c of the camera (eye transformation)................ 107 4. The transformation from Frame S c to the camera film (viewport transformation)................... 109 5. Putting the pieces together ................................................................................................. ....... 110 (e) A tour of Maple code which comput es angle-space and screen-space analemmas..................... 112 (f) An Analemma Ga llery ................................................................................................................. 118 1. Analemmas for an observer at the North Pole .......................................................................... 118 2. Analemmas for observers at other latitudes in the northern hemisphere. ................................. 121 (g) Analemma case study from the Crimea (Ukraine) ...................................................................... 128 Appendix A: A Note on Various Arctangent Functions ............................................................... 132 (a) arctan2Pi(y,x) and comparison with arctan(y/x) and arctan(y,x) ................................................ 132 (b) arctan2Pis(y,x)............................................................................................................ ................. 137 Appendix B: The Active and Passive Views of Vector Rotation................................................... 141 (a) Active view.................................................................................................................................. 141 (b) Passive view ............................................................................................................... ................. 141 (c) Summary.................................................................................................................... .................. 142 Appendix C: Shape of the north pole analemma when e = 0 ......................................................... 144 (a) Parametric equation for the analemma ....................................................................................... .144 (b) Locus of the analemma in angle space: symmetry..................................................................... 146 (c) How β varies along the analemma............................................................................................... 148 (d) Limit of analemma shape for small θtmax : Gerono!................................................................... 149 3 Appendix D: Calculation of how much the analemma is tipped on film....................................... 151 (a) Calculation of the tip angle from vertical................................................................................. ...151 (b) Physical explanation of the tip angle behavior............................................................................ 155 References............................................................................................................................................ 157 Overview 4 Overview This paper gives an elementary presentation of the kinematics of the earth-sun system with a few basic applications. An interested reader would do well to r ead the Summary first, to acquire a familiarity with the general thread of the paper and with th e symbols used for the variables involved. The term kinematics refers to the set of variables used to describe the mechanics of the sun-earth system. Much of this paper is concerned with findi ng expressions for and relations among these variables. One might refer to the way the variables vary in time ("motion") as the dynamics of the earth-sun system characterized by the kinematics. The paper is almost completely self-contained and derives from scratch precise mathematical expressions for all quantities of inter est. The only math required is trigonometry, some geometry, the use of 3x3 matrices, and a basic knowledge of spherical coordinates. A small handful of derivatives and simple one-dimensional integrals are encountered al ong the way. One first-order differential equation appears in Section 4. Due to the high historical cost of creating and printing graphics, authors of technical books often short-changed their readers with regard to "figures", having to replace proper detailed drawings with complicated text descriptions of what those drawings would look like. The author of this paper, though hardly an artist, benefits greatly from reasonably pr ecise and well-labeled drawings, and for that reason has produced many of them for this document, even though that increased th e file size to several megabytes. Most were quickly created in Visio, though some were cribbed fro m the web (with reference always given). Our drawings are "line drawings", rather than elaborate renderings of intersecting curved surfaces. A particular graphic trick used several ti mes is the superposition of two screen-clip images in Visio where the front one is made partially transparen t, allowing us to compare one of our graphs with another graph found elsewhere. Finally we come to Maple. The author uses a 1995 ve rsion of the Maple computer algebra system that works beautifully and is unencumber ed with GUI enhancements. The internal Maple system is deep and powerful and one gets a lot done with a few lines of code (see author's Maple User's Guide). Maple code is readable by a person familiar with any programming language. The main use of Maple in this paper is to graphically display how interesting quantities vary with respect to one or maybe two parameters. The analemmas in the last section were created by ge nerating time-strobed data arrays for analemma coordinates and having Maple plot the data with its poi ntplot routine. Our Maple code is freely available by email to any interested party (it is presented only in screen-clip form in this document). It will run in the worksheet of the current Maple release, or can be translated into Mathematica, MATLAB, Octave (freeware), etc. It is hard to imagine how "those w ho went before" investigated things with only hand calculations. It does generate a sense of respect. • Equation numbers for equations which are repeat s of the original equation are italicized. • Unit vectors are usually marked with a "hat", as in x^ . Exceptions are the e n and e'n unit vectors. Summary 5 Summary This is an abnorm ally detailed summary which the author feels is useful to introduce the names of the many kinematic variables, and to make visi ble the meandering thread of the paper. Section 1 presents the basic facts of the earth's orbit around the sun, then shows how this picture is transformed into one in which the sun apparently orbits the earth. These are the heliocentric and geocentric views of the earth-sun relationship. Th ree different angular velocities are involved, called Ω0, ω and ωs. The earth's rotation axis is tilted by θtmax ≈ 23o relative to the plane of its slightly elliptical orbit which has eccentricity e ≈ .017 and whose major and minor axes are rotated by an angle φ per ≈ 13o relative to another pair of axes associated with seasonal moments of the year known as solstices and equinoxes. The equation of the ellipse is stated and Kepler's "anomalies" are described. The basic heliocentric picture is this (Fig 1.3), which shows two azimuthal angles φ' and ψ for the earth's position, and another angle α. Angle φ' is measured from March equinox, where time τ = 0, while ψ is measured from the point of closest approach (perihelion) where a different time t = 0. The figure shows that α = π/2 - φ per (1.9) φ' = ψ - α . (1.10) t = τ + tMe (7.2) where t Me is the time it takes for the earth to move from perihelion to March equinox. Section 2 describes the sun's apparent daily path in the sky at the two solstice and two equinox days of the year, as viewed by an observer at the North Po le. The two equinox paths are the same, so these paths form three circles in the sky which are parallel to th e earth's equatorial plane. On any day of the year, there is one circle on which the sun rides, and it li es somewhere between the extremal solstice circles. Section 3 introduces the celestial sphere as a visualization tool for describing the angular location (declination = latitude, right ascension = longitude) of distant astronomical objects as well as local objects Summary 6 like the sun and planets whose projections move on this sphere. The three special circular daily paths of the sun (solstices and equinox) are projected onto this celestial sphere and also onto the earth, where they become the Tropics of Cancer and Ca pricorn and the Equator. These thre e projected celestial sphere paths of the sun are later plotted from the point of view of inhabitants of the earth at different latitudes, using a few lines of Maple code. These abstract circles, when clipped to the local tangent plane of the earth, become the familiar sun paths of everyday life. Afte r a brief discussion of the radius of the celestial sphere, the annual path of the sun on the celestial sphere is described, a path known as the ecliptic. The earth's rotation axis precesses slowly and this h as implications which are duly noted. Finally, the significance of the Arctic and Antarctic Circles is explained. Section 4 shows how the angular velocity Ω = ψ• of the sun along its elliptical path (and along its circular projected ecliptic path on the celestial sphere) varies slightly with sun position due to orbital eccentricity (Kepler's second law). This causes the angular position ψ of the sun to differ slightly from the position of a circular-orbit sun. The azimuthal position ψ (t) is given approximately by (t=0 is perihelion Jan 3), ψ(t) = Ω 0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t ) (4.31) which is one of two principle ingredient s in the so-called "equation of time". Section 5 discusses the way the earth's tilt angle toward the sun, called θt, varies during the year between θtmax = 23o at June solstice and -23o at December solstice. At either equinox, θt = 0. A physical interpretation of θt is developed in which θt really is the angle of "earth's tilt toward the sun". This same angle θt is in fact the declination of the sun on the celestial sphere ecliptic and is shown to be θ t(t) = sin-1 { sinθtmax sin(φ')} = sin-1 { sinθtmax sin(ψ-α) } . (5.5) Section 6 obtains expressions for the sun's position on th e ecliptic first in Cartesian coordinates r = x,y,z and then in spherical coordinates R, θ,φ where R is the radius of the celestial sphere and θ ,φ are the sun's polar and azimuthal angles on the celestial sphere. Since θ t is the sun's declination (latitude), θ = π/2 - θt. This angle φ (the sun's right ascension) is measured from a certain point on the celestial sphere (March equinox) where τ = 0. Since the ecliptic is tilted on the celestial sphere, angle φ and the angle φ' of Section 1 are not the same and we find that, φ(t) = tan -1[cosθtmax tanφ'(t)] = tan-1[cosθtmaxtan(ψ-α)] . (6.13) A "mean fictitious sun" is imagined to exist in the sky which operates with θt = 0 (no tilt) and e = 0 (circular orbit). The azimuth of this fictitious sun is shown to be φf(t) = Ω0t - α . (6.14) The difference Δφ = φ(t) - φ f( t ) (6.15) Summary 7 is called "the equation of time", which then includes the azimuthal variation effects of both the tilt θtmax and the eccentricity e. Since φf(t) is a simple linear function of time t, it can be regarded as time on a clock, and then Δφ describes the variation of "sun time" from "clock time", which is the real meaning of the phrase "the equation of time" -- it is a temporal correction. Sundials follow φ(t), whereas mechanical or electronic clocks follow φ f(t). Section 7 considers an observer's reference Frame S' wh ich is located on the surface of the earth at colatitude (polar angle) θ1 and longitude φ1 east of Greenwich. In the celestial sphere frame of reference, the location of the earth's Frame S' origin is described by a vector b(t) which moves (on a cone) because the earth rotates. The polar angle of b is θ1 and its azimuth is φb(t). This azimuth φb(t) can be projected out onto the celestial sphere where it becomes in eff ect the right ascension of the observer site. Since the earth rotates at sidereal rate ωs, one has φb(τ) = ωsτ + φb 0 (7.3) where then φb0 is the right ascension the observer site had at the last March equinox τ = 0. The value of φb0 is given by φb0 = φ1 + φG0 where φG0 is the right ascension that Greenwich had at τ = 0, and this in turn can be obtained from the publis hed GMT time of March equinox. We mention all this detail just to get all the symbols laid out on the table. The sun's position in Frame S' is computed first in Cartesian coordinates r' = x',y',z' and then in spherical coordinates R, θ',φ' : x' = R cos θt sin(φ –φb) y' = R [sin θt sinθ1 – cosθt cosθ1 cos(φ –φb) ] z' = R [sin θt cosθ1 + cosθt sinθ1 cos(φ-φb) ] (7.14) and θ' = cos-1[sinθt cosθ1 + cosθt sinθ1 cos(φ-φb)] φ' = tan-1 ( [sinθt sinθ1 – cosθt cosθ1 cos(φ –φb) ] / [cosθt sin(φ –φb)] ) . (7.29) In these equations, θ t and φ are given by (5.5) and (6.13) quoted above. Section 8 replaces the azimuth φ' in (7.29) with Φ' = π/2 - φ' so the second of equations (7.29) describing the sun's azimuthal position in Frame S' becomes Φ' = π + tan-1 ([cosθt sin(φb-φ)] / [– sinθt sinθ1 + cosθt cosθ1 cos(φb-φ) ] ) Φ'f = π + tan-1 ( [sin(φb-φf)] / [cosθ1 cos(φb-φf) ] ) . (8.3) where the second line is for the fictitious sun. For northern hemisphere observers at latitude above 23 o, Φ' is an azimuth that "runs with the sun" at all times of the year, and solar noon always occurs at Φ' = π which is "due south" in Frame S'. At the North Pole ( θ1= 0), we find ΔΦ' NP ≡ Φ'NP – Φf'NP = – (φ – φf) = - Δφ . (8.6) Summary 8 At the north pole, Φ' "runs with the sun" while Φf "runs with the clock" and their difference is again the equation of time, albeit negated. The notions of Lo cal Mean Time (LMT) and Local Sundial Time (LST) are defined, and we arrive then at our ve rsion of the "official" equation of time, EOT(t) ≡ [LST – LMT]minutes = – (720/ π) { tan-1[cosθtmax tan(ψ-α)] – (Ω0t-α) }. (8.16) LST varies from LMT by a maximum of about 15 minutes during the year. Section 9 uses the equations developed above to compute the Local Mean Time (LMT) of (actual sun) noon and midnight in Frame S' on any day of th e year. Similarly, th e times of sunrise (t' sr) and sunset (t'ss) are obtained relative to noon, so the nominal length of day equals t' ss – t'sr, with time of day light being up to 8 minutes longer (ball park). The Frame S' azimuth for sunrise (φ 'sr) is computed and plotted for different latitudes and days of the year. Finally , an expression is found for the noontime elevation of the sun. These results are: LMT(noon) = [12:00 + (720/ π)Δφ] m o d 2 4 (9.5) LMT(midnight) = [00:00+ (720/ π)Δφ)] mod 24 . (9.6) t'sr = – (1/ω ) cos-1[– tanθt cotθ1] < 0 // sunrise t'ss = + (1/ω) cos-1[– tanθt cotθ1] > 0 // sunset (9.15) tday = 2 t'ss = (24/π ) cos-1[- tanθt cotθ1] hours (9.16) tanφ'sr = (tanθ t/sinθ1) / 1 - tan2θt cot2θ1 (9.17) θ1+θt range: 0 π/2 π θel noon = : no noon θ1+θt π - (θ1+θt) no noon (9.20) These expressions are approximations which assume (correctly) that θt(t) and Δφ(t) vary slowly during a day so one can approximate these functions at precision time t by their computed values at the start of a day of interest. Section 10 is concerned with the amount of solar energy per second loaded into a square meter of black ground on the earth (insolation) as a function of latitude and time. If the solar power flux is S watts/m 2, then the amount captured by that black square meter of ground is Scos θ', a fairly trivial fact known as Lambert's Law (θ ' is the polar angle of the sun given by th e first of equations (7.29) above). Daily insolation is Scos θ' averaged over a day, while annual insola tion is this same quantity averaged over a year. Various plots are generated by Maple code and are then compared with other references. The same approximation mentioned in the previous para graph is used in this insolation work. Section 11 introduces an imaginary sphere called the sun sphe re which is a version of the celestial sphere which rotates once a year such that the mean fictitious sun remains at a fixed point on the sun sphere all year long. The actual sun then moves in a pattern in the neighborhood of this fixed point on a path that is called the solar analemma. A qualitative description of the analemma viewed from various locations is presented. The analemma path can be exposed by the use of stroboscopic photography of the real sun wherein a multiple exposure image is accumulated over the period of a year by taking a heavily filtered Summary 9 exposure at the same time each day (or every nthday) in the general direction of the fictitious sun's fixed point. The path that is the analemma is computed in a few dozen lines of Maple code in both angle-space (θ',Φ') and in the screen-space (x s,ys) of the film in the camera. An analemma "gallery" then presents Maple-computed images of the analemma as it woul d appear shot from various positions in the northern hemisphere. The final section treats a certain Crimean-s hot analemma as a case study, and the computed analemma is compared with the one captured on film. Appendix A describes a special arctangent function impl ementation which returns a well-defined result in the range 0 to 2 π and introduces the special notation tan -1( [y] / [x] ). This notation is used in the second line of (7.29) above and the corresponding func tion is implemented in Maple for evaluation of φ'. Several other arctangent functions are also consider ed and are then compared on a branch diagram. Appendix B discusses the active and passive views of the rotation of a vector. Appendix C studies the shape of the analemma for a circula r-orbit earth. It is shown to be symmetric in both directions, and in the limit of small θ tmax it becomes a curve known as the lemniscate of Gerono. Appendix D develops a chart to show how much a photograp hed analemma will be tipped left or right on the film, assuming a camera platform whose pointing direction is set by certain tilt and pan angles. The final section contains a few References. Section 1: Geometry of the Earth's Orbit 10 1. Geometry of the Earth's Orbit (a) The Heliocentric Picture The following "heliocentric" (sun-centered) camera view of the earth 's orbit is taken from above the earth's orbital plane, in the sense th at the earth's north pole is up, a nd the camera is elevated perhaps 30 degrees above that plane. Note the basic fact that the ear th rotates about its axis in the same sense that the earth rotates around the sun: both rotations are coun terclockwise viewed from above. This is why the noon-to-noon day is a little longer than the sidereal (rela tive to stars) day. Starting at noon, after rotating one sidereal day the earth has to rotate an extr a 4 minutes to reach noon again because the earth has moved ahead a bit in its orbit. The average noon-to-noon day (mean solar day) is exactly 24 hours (by definition), whereas the sidereal day is about 23 hours and 56 minutes ( ≈ 23.93447). A sidereal year is 365.25 mean solar days long (365.256363004 in 2000), a nd is the time for the earth to completely traverse its orbit, returning to its starting point relative to th e distant stars. Since time is now based on an atomic standard, and since the earth has tiny variations in its rotational speed, the mean solar day might vary from 24 hours by 1 second per year (called a leap second), but we ignore this 10 -8 fine detail. The sun is located at one of the two focal points of the earth's elliptical orb it. The distances for closest approach (perihelion) and farthest (aphelion) are stated in millions of km. F i g 1 . 1 The slanted arrows indicate the angular momentum vector of the rotating earth. The times of the year when this tilt arrow points toward or away from th e sun are called solstices, and the times at which the arrow points neither toward nor away from the sun ar e called equinoxes. When the arrow tilts toward the sun, earth inhabitants in the northern hemisphere ex perience "summer". At the summer solstice, the mean daily insolation reaches a maximum (but not close to the equator), as we show in some detail in Section 10. This is basically due to the fact that the sun is close to overhead much of the time during the day. At the same time, the insolation reaches a minimum in the southern hemisphere, so Argentineans experience winter. In hemisphere-neutral language, it is best to refer to "December so lstice" and "June solstice", and as well "March equinox" and "September equinox". The term "vernal equinox" is heavily embedded into astronomy literature meaning the March equinox, though it is not verdant spring in Argentina at that time. Section 1: Geometry of the Earth's Orbit 11 The plane of the earth's orbit around the sun is calle d "the ecliptic plane", but one more normally uses this term when thinking of the plane of the sun's apparent orbit around the earth. It is the same plane relative to the stars in either case, ju st viewed by different observers. The plane of the earth's equator is tilted by an angle of about 23o relative to the ecliptic plane of the earth's orbit. For reasons seen later, we refer to this angle as θtmax, and its official name is "the obliquity of the ecliptic". Here both planes are seen edge-on, Fig 1.2 Viewed from directly above, the elliptical orb it would appear very close to circular, having eccentricity e = c/a = .0167. The ratio of the extremal earth-sun distances is about 152/147 ≈ 1.034, resulting in a solar energy flux ratio of (1.034)2 = 1.069 between the two solstices (ignoring the 13o in Fig 1.1). Solar flux is therefore 6.9% greater at summer sols tice in the southern hemisphere (Dec) than it is in the northern hemisphere (June). To see more clearly the position of the earth at the equinoxes, we now swing the camera up and view the ecliptic plane from directly above (ellipse eccentricity exaggerated). The earth's angular position can be measured by angle φ' from March equinox ( τ=0), or by angle ψ from perihelion (t=0) where τ and t are two different time variables with different zero points: Fig 1.3 Section 1: Geometry of the Earth's Orbit 12 The official way of measuring the 13o angle is shown below as the "longitude of perihelion" ≈ 283o. Here is some data from http://aom.giss.nasa.gov/srorbpar.html showing how this angle, the tilt angle and even eccentricity change over a small number of years: [ 283o = 270o + 13o ] Table 1.1 For 2012 the angle is therefore 13.101o. Here is data at two larger time scales, Table 1.2 Table 1.3 The 13 o angle indicates a relationship between the earth's elliptical orbit and the direction of its tilt vector. This angle φ per ≈ 13o along with the tilt angle θtmax ≈ 23o and the eccentricity e ≈ .017 all vary in time as part of the complicated solution of a gravitational N-body problem involving the earth, sun, moons, planets, and asteroids. In this solution, the total change in a parameter like φ per can be Section 1: Geometry of the Earth's Orbit 13 approximately viewed as being a sum of changes each associated with a certain physics "mechanism", some of which are long term and some short term relative to human time scales. One such mechanism causing a change in φ per (called astronomical precession, see Section 3 (d)) is that the earth's rotation axis precesses in a complete conical circle every 26,000 years. Although this precession does not affect θtmax, other mechanisms do, and θtmax cycles between 22.0o and 24.5o with a period of 41,000 years. Another mechanism is that the ellipse of the earth's orbit is itself rotating in its plane relative to the stars, completing a rotation every 112,000 years. Each of the two extremal points of an ellipse is called an apsis (perihelion = periapsis, aphelion = apapsis, plural apsides), and the line connecting these points is the line of apsides, so it is this line that rotates relative to the stars (called apsidal or perihelion precession). Yet another mechanism is that the ecliptic plane of the earth's orbit wobbles back and forth with a similarly long period. Fu rthermore, the eccentricity of the orbit varies as well (range .005 to .058) with cycl es of 413,000 years and 100,000 years. These "mechanisms" are all due to forces of the moon and other objects acting on the earth. They were studied by the Serbian geologist/climatologist Milutin Milankovic (1879-1958) who was interested in the effect of these changing orbital parameters on climate. His theory called orbital forcing shows how changing orbital parameters could cause ice ages on the earth. Short term mechanisms also affect the paramete rs. As the earth's rotation axis precesses, it also nutates ±20 arc seconds with a 19 year cycle due the perihelion precession of the moon's elliptical orbit around the earth. Another form of nutation is associated with the "wander" of the earth's poles, such as the 433 day Chandler wobble which moves the poles ±15 feet. These effects are caused by the non-rigidity and non-roundness of the earth. Ignoring all these effects due either to their long time scale or their small size, we may regard the earth's rotation axis as pointing in a direction which is fixed relative to the stars. (b) Development of the Geocentric Picture First, consider a scenario in which some object A moves counterclockwise in a closed path about object B which is at rest. Neither object is rotating about its own axis. Each object has a glued-on set of axes. This scenario appears on the left below (black) which is the view taken by an observer on object B in Frame S: Fig 1.4 If we call the path p(t), t being time, then on the left we have in Frame S, rB(t) = 0 // B is at rest rA(t) = p(t) . // A moves around the path p(t) Section 1: Geometry of the Earth's Orbit 14 Now we go to a new frame of reference S' defined by a time-dependent translation by - p(t) relative to Frame S, r'(t) = r(t) – p( t) . In this new frame of reference we find that the positions of A and B are given by, r A'(t) = rA(t) – p(t) = p(t) – p( t) = 0 // A is at rest at the origin of Frame S' rB'(t) = rB(t) – p(t) = 0 – p(t) = – p(t) . // B moves around the path -p (t) In Frame S' the object A is now at rest, and object B is traversing a negated path – p(t) which is an inversion of the original path through the Frame S orig in (blue dotted path). This situation is shown on the right above. If we now assume that p(t) is a planar path, as suggested by the drawing, then inverting the original path is the same as rotating it 180 degrees about the Frame S origin in the plane of paper. This is so because then R z(π) p(t) = - p(t) since p z(t) = 0. The two points of this exercise are the following: (1) To obtain the picture on the right, rotate the path p(t) 180 degrees about the Frame S origin on the left (to get the dotted blue path), then translate points A and B together by - p(t) which causes B to lie on the inverted path. (2) If A is going counterclockwise around B on th e left, then B is going counterclockwise around A on the right. The sense of the path motion is not reversed. The new path is -p (t), not p(-t). Example. Suppose path p(t) is a tilted ellipse and suppose A lies at a focal point. Application of the above procedure takes the left picture below (A is at rest) into the right picture (B is at rest) : F i g 1 . 5 If A is the non-rotating sun and B is the non-rotating earth, then the left side picture is the heliocentric one and the right is the geocentric one. Of course for the earth-sun application the eccentricity is much less than shown here. Since the sun is so massive, it can be regarded as a fixed point A in the left picture. The Section 1: Geometry of the Earth's Orbit 15 conclusions are that in the geocentric picture on the right, the earth lies at a focus of the elliptical sun orbit, and the sun moves counterclockwise around the earth. Notice where C and D end up on the right. Without further ado, we now transform Fig 1.3 to the geocentric view : Heliocentric Picture Geocentric Picture F i g 1 . 6 (c) Relation between the three angular velocities Let diy = 365.256363004 be the number of days in a year. In one (sidereal) year's time, the earth experiences diy number of 24 hour periods, normally re ferred to as "mean solar days". But during this same time, the earth has in fact rotated diy+1 tim es relative to the stars. If a person faces the noon sun each day and is the subject of a stroboscope movie f ilmed from above the ecliptic, with a frame shot each day at that person's noon (shots are 24 hours apart) , that movie will show that person making one complete revolution during the year. Had that movie shown a static person, one would say that the earth rotated exactly diy times during the year. But in fact the earth rotates a little more than 360o each 24-hour day to keep that person facing the sun and then the total rotations is diy + 1. So we have T year = diy*T T ≡ 24 hours, length of solar day Tyear = (diy+1)*T s T s ≡ length of sidereal day (1.1) Therefore, Ts = [ d i y / ( d i y + 1 ) ] T ( 1 . 2 ) = [365.256363004/366.256363004]*24 = 23.93447 hours which agrees with the number quoted above. The (mean) rotation rate of the earth around the sun is Ω 0 ≡ 2π/Tyear ( 1 . 3 ) Section 1: Geometry of the Earth's Orbit 16 The solar-day angular frequency ω is, using (1.3) then (1.1), ω ≡ 2π /T = Ω0 (Tyear/T) = Ω 0 * diy => Ω0 = ω/ d i y . ( 1 . 4 ) Inverting equation (1.2) then multiplying by 2π gives 2π/T s = [1 + 1/diy] 2 π/ T or ω s = [1 + 1/diy] ω . ( 1 . 5 ) Then using (1.4), ωs = ω + Ω0 (1.6) which is the main result of this little section ( ω s is a little larger than ω). As for units, most often we shall use ω = π/12 radians/hr and Ω0 = 2π radians/year. (d) Equation of the ellipse We ref er to the 13o angle in Fig 1.3 as φper (perihelion). Notice in this same figure another angle ψ which is the usual angle one uses when describing an ellip se in polar coordinates (a = semimajor axis, e = eccentricity) r = a(1-e2) / (1 + e cos ψ ) ( 1 . 7 ) so that ψ = 0 gives the shortest (perihelion) value of r. The figure shows that ψ - φ' = π /2 - φper => φ' = ψ + φper - π/ 2 ( 1 . 8 ) To reduce symbol count later on, we define α ≡ π/2 - φper // α ≈ 77o (2012) (1.9) and then φ' = ψ - α . ( 1 . 1 0 ) Section 1: Geometry of the Earth's Orbit 17 (e) The Anomalies When Kepler wrote up his laws in 1609 and 1619, he used the term "anomaly" for an angle that exhibited some unusual behavior. As we shall find in Section 4, the angle ψ (and therefore φ') varies with time in a non-linear fashion; it is linear only when the elliptical orbit is circular. The angle ψ (Fig 1.3 or Fig 1.6) measured from perihelion is called "t he true anomaly". Kepler's third law says that the period of an elliptical orbit (for two given objects) is a function only of the semimajor axis a (T2 ~ a3). One can circumscribe a circle of radius a around an ellipse and if the orbiting object had that circular orbit, the orbital period would be the same as that of the inscrib ed elliptical orbit. The angular position along this alternate circular orbit is linear in time and is some times called M(t), the "mean anomaly", so one has then M(t) = 2 π(t/T) where T is the orbital period. If the elliptically orbiting object's position is projected "vertically" to that circumscribed circle, a third angle appears, E(t), called "the eccentric anomaly". After time t elapses from perihelion, the three anomalies ψ, M and E all have different values. Fig 1.7 Kepler showed that M = E - e sin(E) and tan( ψ/2) = ( 1+e / 1-e ) tan(E/2), so this provided a numerical method of finding ψ(t) by first solving the transcendental equation for E. In Section 4 we shall calculate ψ(t) using an approximation method good for small e, a nd we shall make use only of this true anomaly ψ. Ref [1] Chap 3 (3-18) and (3-20) provides a ge neral formal solution to the Kepler problem. Section 2: The Sun's Apparent Path in the Sky 18 2. The Sun's Apparent Path in the Sky As noted above, the earth tips toward the sun at June solstice, and away at December solstice. This figure shows a composite of three different earth positions relative to the sun, F i g 2 . 1 For the moment, let's call the maximum earth tilt 23 o (Table 1.1 shows 23.4382 for 2012). By the term "tilt" we mean "amount of tilt toward the sun" which is 0 in the middle image above. This tilt will later get the name θ t, and the maximum value for θt is 23o ≡ θtmax. When the earth is in the orientation shown at the left, a June solstice observer on the north pole sees the sun elevated 23o above the horizon. As the earth turns, that sun appears to sweep around in a circle, holding its 23o elevation angle above the horizon. An equinox north pole observer instead sees the sun sweep right around at the horizon, with half the sun above the ho rizon and half below (the sun is very far away). A December solstice north pole observer never sees the sun, but if this observer had X-ray vision the sun would appear to be going in a circular path 23o below the horizon. In Fig 2.1 Frame S is a coordinate system: (1) w hose origin is at the center of the earth; (2) whose axes are aligned with the earth but don't rotate with th e earth; (3) whose axis directions are fixed relative to the stars; (4) whose z axis is along the red line and whose x and y axes lie in the blue equatorial plane. Since the earth moves slowly around the sun, it has some small acceleration relative to the stars, so Frame S is not quite an inertial frame. [The red lines in Fig 2.1 point in different directions only because the "camera" is attached to the earth-sun line and is at different orientations for the three earth images. See Fig 5.1 and discussion there. ] We now construct a rotating Frame S' glued to th e north pole that rotates with the earth. Later we shall move Frame S' to other locations on the earth. If one hops aboard Frame S' wh ich is rotating at rate ω relative to Frame S, then objects which were at rest in Frame S appear to be rotating at rate - ω relative to Frame S'. In each of the three di fferent situations shown above, the ω vector lies along the red line. What one sees from Frame S' in these three situations is shown below, where the earth's size is highly exaggerated relative to its distance from th e sun (earth should really be a dot), Section 2: The Sun's Apparent Path in the Sky 19 F i g 2 . 2 During the cycle of the year, the green daily circular sun path moves. If we start at a June solstice, it is the upper green circle (but the circle is viewed edge- on so appears as a green line segment). Moving toward September, that circle drops down and at Septem ber equinox becomes the middle green circle above. Moving toward December, the sun continues to move down until the December solstice when it follows the lower green circle. During the next 6 months the path of the sun reverses and ne xt June it is at the top again. So all possible sun paths are bounded by the upper and lower green circles. The daily sun rotates clockwise as viewed from above the north pole (sun rises in the east). For the positions of the sun marked on the right, the sun is moving toward the observer. In reality, the green circular daily sun paths shown above are turns of a dense continuous green slinky spiral having 365.25 nearly-horizontal turns (half of the turns spiraling down, the other half spiraling back up), and we have drawn three of the spiral turns. In the next Sections, we shall refine this brief qualitative description of the sun's apparent motion in the sky and shall provide appropriate equations describing this motion. Section 3: The Celestial Sphere 20 3. The Celestial Sphere and the Sun' s Daily and A nnual Paths Thereon (a) The Celestial Sphere and the Tropics of Cancer and Capricorn The celestial sphere is an imaginary spherical surface surrounding the earth at some unspecified radius R upon which one imagines that the sun and other astronomical objects are replaced by their projections onto this sphere. An infinite ray is drawn from the center of the earth through an astronomical object, and the projection point is the intersection of this ray with the celestial sphere. Astronomical objects could be inside or outside the celestial sphere. The "marked up" (with star projections) celestial sphere is at rest with respect to the distant stars (Frame S), and the ear th rotates within it. Of course in our Frame S' glued to the north pole of the earth, the earth is at rest and the celestial sphere rotates clockwise around the earth causing stars to "rise in the east" and "set in the west". We shall now draw the projections of the three green circular sun paths of Fig 2.2 onto both the earth and onto a celestial sphere. The projections of sun pa ths onto the earth are shown here in blue, while projections onto the celestial sphere are shown in red, and we slightly displace the equatorial projections just so they are visible : F i g 3 . 1 The two extremal blue circles on the earth are at latitude ± 23 o. These latitude circles have the special historical names Tropic of Cancer in the north, and Tropic of Capricorn in the south. The history of these names is obscure, but the word "tropic" indicates a "turning point", and indeed these extremal circles mark the turning points of the sun path (project ed onto the earth) as the yearly cycle progresses, as described in the previous Section. The "t ropics" lie between these bounding latitudes. . // OED2 The two extremal red circles on the celestial sphere are of course also at latitude ± 23o on the celestial sphere. The radius of the celestial sphere (discussed be low) is supposed to be much larger than the radius of the earth, so to make the above picture more accura te, we redraw it with the earth shrunk to a small dot. Section 3: The Celestial Sphere 21 F i g 3 . 2 A very convenient model for thinking about the sun paths is to think of them as being those red circles on a celestial sphere (usually called the celestial sphere). In the above drawing, the radius R of the celestial sphere is perhaps 25 million miles, sinc e the sun is 93 million miles away. Below we are going to ask the question: What do the three red circles look like from various observation points on the earth? We want an exact answer to this question. For an observer on the north pole, the answer is clear and is represented by the left pi cture below (the winter circle is not visible to this observer), while for an observer on the equator the answ er is also clear and is represented by the right picture below. In these two pictures, the vertical direct ion represents "up" for the observer, south is to the left and east is toward the viewer. Observer at North Pole Observer at the equator Fig 3.3 In these drawings the celestial sphere is suggested by th e dim gray circles. These gray circles form a half sphere, and the flat bottom of the half sphere coincid es with a plane tangent to the earth at the observer's location. This plane contains the observer's horizon, so portions of the red circles below this plane are not visible to the observer. The ending points of partial red circles (as on the right) are identified therefore with sunrises and sunsets (the near ones are sunrises since east is toward the viewer). Section 3: The Celestial Sphere 22 (b) Why is R >> R Earth for the celestial sphere? The celestial sphere construct simply does not "work" if it is very close to the earth. For example, suppose the celestial sphere were 500 feet above the surface of the earth. One could imagine for some observer the stars projected onto this sphere and their positions ma rked with paint dots on the otherwise transparent sphere. As the observer moves to a different location, those paint dots no longer align with the stars they represent. The projections of the stars for this cel estial sphere are dependent on the location of the observer on the earth. On the other hand, if the mathematical celestial sphere is placed 100 earth radii out into space, then as the observer wanders from A to B, the paint dots made at A will be pretty close to the stars seen from location B. Maybe the maximum angular error will be 1/100th of a radian or about 1/2 degree (parallax error). As the celestial sphere is mentally moved further away, this error becomes smaller, and in the limit R→∞ the paint dots are in the same place for all observers on the earth (the stars are imagined to be infinitely far away). So we wa nt R large to prevent parallax error for viewers on the earth, then a single celestial sphere map of the stars can be used by everyone. A separate but related issue concerns the distan ce of observed objects from the earth. A balloon at 500 feet could be projected onto the celestial sphere, but its projection must be different for different observers, and the celestial sphere is pretty useless in this case. The moon is about 61 earth radii away, so if it were projected onto the celestial sphere, its posit ion would be the same value for all earth observers with perhaps 1 degree of error or less. The sun is 24,000 earth radii from the earth, so the parallax error for any earth observer is then very small, perhaps .002 degrees. Eventually we will set R = 1 unit and have mind perhaps R Earth = .0000001 unit. Section 3: The Celestial Sphere 23 (c) Annual apparent motion of the sun around the earth: the ecliptic Viewed from directly above the ecliptic plane, the pa th of the sun around the earth was shown on the right side of Fig 1.6, which we reproduce here as its own sep arate figure. The elliptical orbit of the sun lies in the plane of paper and the earth's tilt arrow points 23o out of the plane of paper and to the right. F i g 3 . 4 We have added in red the intersection of the sun's ecliptic plane with the celestial sphere, which sphere we now show being outside the sun's orbit. Since the earth lies at the center of the celestial sphere, this intersection is a great circle on the celestial sphere. This red circle lies in the plane of paper above and is commonly referred to as the ecliptic . The March equinox is the azimuthal reference point on both the ecliptic and on the earth's elliptical orbit. The angle φ' from the ecliptic reference point to the sun's projection on the ecliptic is the same as the angle φ ' going from the elliptical orbit reference point to the actual sun (similarly for ψ). The elliptical orbit is thus projected onto the circular ecliptic. If we now think of the red disk (with earth) as a fl at damper disk hinged on the vertical axis, we can rotate that disk 23o so that the June solstice rises above the plane of paper, the December solstice drops behind the plane of paper, and the two equinoxes (being on the hinge) stay in the plane of paper. Such a rotation causes the angular momentum vector of the earth to point directly at the viewer. This situation would be precisely the view of the earth-sun system from the north pole of the celestial sphere. Here then is a nice 3D representation of the celestial sphere and its red ecliptic, where June solstice is the high point on the ecliptic on the right. One now sees the abovementioned damper disk in its rotated position. Section 3: The Celestial Sphere 24 Fig 3.5 http://astro.wsu.edu/worthey/ast ro/htm l/lec-celestial-sph.html The ecliptic circle intersects the celestial sphere equator at the two (hinge) equinox points on opposite sides of the celestial sphere. The March equinox (nor thern hemisphere's vernal equinox) is used as the azimuthal reference on the celestial sphere. This azimuth may be rega rded as longitude on the celestial sphere, measured in the same sense that longitude is measured on the earth east of Greenwich. To avoid confusion with various other azimuths, the longitu de or azimuth on the celestial sphere is called right ascension . It could be measured in degrees or hours, and in the latter case one refers to the "hour angle" of the azimuth, as in the picture above. Simila rly, latitude on the celestial sphere is called declination and runs from +90 o at the north pole to -90o on the south pole with 0o at the celestial equator. http://hyperphysics.phy-astr.gsu.edu/hbase/eclip.html Fig 3.6 A star on the celestial sphere located at a certain right ascension will "ascend" from the eastern horizon into the night sky (as viewed from the earth) at a time associated with that star's right ascension. The term "right" has a historical meaning, but it is best to think of it meaning the angle is measured to the "right" of the vernal equinox as shown in the above figure. Section 3: The Celestial Sphere 25 If Fig 3.5 were "live", one would see the celestial sphere fixed, the earth rotating counter clockwise once per day (viewed from the north celestial pole), and the sun slowly traversing the red ecliptic circle counterclockwise once per year. The sun traverses the ecliptic circle over a peri od of 12 months, and each month's segment of the ecliptic lies "in" one of the 12 constellations of the zodiac (circle of animals, like Leo the lion, but only half of the 12 constellations are animals). Roughly in February, for example, the sun is "in the house of Aquarius". The houses don't exactly align with official months. One might wonder how the ancients knew where the sun was on the ecliptic since the bright sun makes the stars on the celestial sphere invisible. They might have done it this way: (1) observe the elevation of the noon sun (sun on meridian) and deduce from one's latitude the corresponding celestial sphere declination of the sun (Eq. 9.20); (2) at mi dnight, observe what stars are located at this same declination on the midnight meridian ; (3) the sun must then be on that same declination on the opposite meridian; (4) consult a celestial sphere map to see where the sun must therefore be. The mass-averaged plane of the planetary orbits is called the invariable plane . The orbits of the 8 major planets all lie fairly close to this invariable plane (< ~ 6 o). For that reason, the orbits of the planets all lie close to that of the ecliptic ( < ~ 7o). Therefore, one looking for the planets in the nighttime sky should look near the path in the sky the sun took that day. Since the moon's orbit is inclined only 5o relative to the ecliptic, when visible the moon will be located in that same region of the sky. So the planets and the moon all wander through the houses of the zodiac. Notice that the annual path of the sun shown in red in Fig 3.5 is very different from the daily paths of the sun which are red circles which align with lines of latitude on the celestial sphere as shown in Fig 3.2. When the sun on the above tilted red ecliptic circle is viewed from our north-pole Frame S' which is glued to the rotating earth, its path over a year is a red spir al (with 365.25 turns) which runs back and forth once between the two extremal daily red circles shown in Fig 3.2. Here is an illustration of that spiral where we draw 12 turns instead of 365.25 turns: Fig 3.7 On any given day, that day's spiral turn is very clo se to lying on a fixed declination circle of the celestial sphere, and that spiral turn is th en the sun path for that day. Here is an $80 commercial celestia l sphere with various stars, galaxies and the milky way projected onto its surface: Section 3: The Celestial Sphere 26 Fig 3.8 http://starlab.com/shop/Cel estial-Globe-SKU1094.htm l The yellow Milky Way shown in Fig 3.8 reflects the pos ition and orientation of the earth relative to the galaxy as indicated in this Britannica picture ( the earth is in the Orion Arm of the Milky Way) : http://www.britannica.com/EBchecked/media/10140 7/Three-views-of-the-Milky -Way-Galaxy Fig 3.9 Section 3: The Celestial Sphere 27 The picture shows that the ring of the Milky Way runs close to the north and south poles of the earth and one can then imagine its location on the celestial sphere (yellow in Fig 3.8). Commercial celestial globe makers face a conundrum: if the sphere is accurate, then all constellations appear inverted when the sphere is viewed from the outside, which is where it is normally viewed from. Sometimes manufacturers transform the true celestial sphere so what one sees from the outside is what one would see on a true celestial sphere from the insi de, which is immensely confusing. Transparent spheres at least allow the possibility of doing it righ t so one can view through the sphere. One could imagine a fancy celestial sphere with a tiny movable camera on the earth which one can hook to one's "PC". A simpler alternative is to just download the freeware Stellarium http://www.stellarium.org/ . (d) Precession effects The reason a tilted spinning top precesses is that gravity applies a torque N = r x m g about the center of mass of the top, relative to the point of c ontact. This tangential torque causes a change d ω = Ndt which is perpendicular to the rotation vector ω of the top and causes the top to precess. Since the earth rotates and is not rigid, it has a larg er radius at the equator than at the poles (by about 21 km). The plane of the earth's fat equator makes some significant angle relative to the earth-moon line, as in this drawing which assumes a set of positions and orientations at some instant in time (the moon- earth line can be at most 5 o out of the ecliptic ). Fig 3.10 As a result, the moon exerts a net torque on the earth due to the gradient of its gravitational field at the earth (the same gradient which causes the bulging of the tides). This torque causes the earth to precess. If the earth were perfectly spherical and uniform, it w ould act as a point mass at its center upon which the moon could exert no torque. The sun creates a torque in this same manner, but the torque is smaller because the sun's distance outweighs its larger mass. The net resu lt is that the sun's effect is about half that of the moon, as is the case for earth tides. The precession rate of the earth's rotation axis is about one degree every 72 years (25,800 year cycle). The earth's angular momentum vector ω points from the center of the earth to its north pole. If one traces out the path of this vector, it rotates around a cone of half angle 23 o once every 25,800 years, as shown in Fig 5.1 below. One consequence of the earth's precession is that the time between March equinoxes is slightly less than a year. After the sun does exactly one revolution about the earth starting at a March equinox, it expects to be at the next March equinox where the earth is again tilted neither toward nor away from the sun. But the tilt axis has precessed 1/72 of one degree (50 arc sec), so the equinox is early (it turns out) by this amount of time : Section 3: The Celestial Sphere 28 Δt = (1/72) deg *[1 yr/ 360 deg] * [ 365.25 day/yr] * [24*60 min/day] = 20.3 minutes (3.1) This effect is called "the Precession of the Equinoxes", and of course the Se ptember equinox and the two solstices experience this same anomaly. The seasona l "year" between solstices is called "the tropical year". The sidereal year is "too long" and if it c ould be shortened a bit, it would match the seasonal year. But as it is, the seasons are slowly sliding backwards relative to a sidereal calendar and in fact slide entirely through the calendar every 25,800 years. The solstices are associated with the high and low points of the ecliptic circle on the celestial sphere, and those points in turn are associated with certain of those zodiac zoo animals, so the 12 members of the zodiac slide along with the seasons relative to the sidereal calendar. But most people don't use a sidereal calendar. Inst ead, they use one that is on average 365.25 days long, since every 4th year (leap year) adds a day on Fe b 29 (the Julian calendar). Now the sidereal year is 365.25 days + 9.17 min long, so the tropical year is then 365.25 days - 11.13 min long (taking off the 20.3 minutes in (3.1)). So the Julian calendar is too long by 11.13 minutes, which means the seasons are sliding backwards through the Julian calendar by 11.13 minutes every year. This adds up to one day after 129.4 years. A scheme was concocted to knock out 3 days every 400 years which is one day every 133.3 years, close enough for government work. The scheme is to not have a leap year every XX00 year, except when XX is a multiple of 4. This plan is called the Gregorian calendar and was set up in 1582 by Pope Gregory. Since Easter was associated with a certa in seasonal date (first Sunday after full moon after vernal equinox) and with a certain calendar date (vernal equinox around March 20-21), the Pope didn't like the idea of these two associated dates drifting apart. It is all a very long story. Since the north pole of the celestial sphere always a ligns with that of the earth, one can think of this precession effect in terms of the earth's rotation axis staying fixed vertically and the celestial sphere (and all the stars attached to it) doing the precessing. This is a very difficult concept to illustrate with line drawings, but here is a crude attempt: F i g 3 . 1 1 Think of the present celestial sphere on the left as a globe with all the stars and the sun's ecliptic painted on it. The usual axis of this globe is the vertical black arrow, but consider the blue arrow shown on the Section 3: The Celestial Sphere 29 left as another axis glued to the celestial sphere. As the celestial sphere precesses, this blue arrow traverses the large vertical cone (half angle 23o), so that 25,800/2 years from now the blue arrow will be as shown on the right. The positions of all stars on our precessed globe are then transferred to the celestial sphere of that future time. Polaris which was the north star moves to the new location shown, 46o down from the celestial north pole. The small black ro tation axis of the earth meanwhile has moved on its cone, which is the pink cone of Fig 5.1. The red ecliptic circle wobbles around with the celestial sphere and after half a wobble ends up as on the right. Every poi nt on the ecliptic, and in fact every point on the original celestial sphere, traverses a circle on the sphere of half angle 23o, Polaris being an example. As the ecliptic wobbles around, its intersections with th e celestial equator move, and these intersections are the equinoxes which are then moving 1/72nd of a degr ee per year on the celestial sphere. To understand this picture, one has to hold a physical globe in tw o hands, imagine a blue arrow sticking out, and then wobble the globe in such a way that the blue arrow traverses the cone. During this precession of the celestial sphere, star charts have to be updated. Since any two orientations of a sphere are related by a specific ro tation R, the coordinates (ephemerides, see below) of all stars at t = t can be obtained from those at t = 0 by applying some uniform rotation R(t) to their coordinates at t = 0. Ephemeris. A table or book of computed future positions of astronomical objects on the celestial sphere as a function of time is called an ephemeris (< Gk. journal, diary, calendar). The positions themselves are called ephemerides, plural of ephe meris. The phrase "ephemeris time" (ET) refers to a certain precision time scale (of which there are now many) which is inde pendent of the vagaries of the earth's (rotational) motion and which agrees exactly with the dynamical time variable "t" in Newton's Laws governing the motions of objects in the solar system, as based on observations of these objects. Ephemeris time was used in ephemeris tables circa 1950-1980, but is now replaced by more elaborate precision times based on atomic clocks, some of which take into account effect s of special (time dilation) and general (gravitational curvature) relativity. (e) Daily sun paths for various observation points In Section 7 below we shall derive the equations used in this section to plot three red daily sun-path circles (two solstice paths, one equinox path) as vi ewed from various observation sites in the northern hemisphere. Here we just display the results for the r eader's possible interest. Recall that "up" is up, south is to the left, and east is toward the viewer. For ea ch picture both the polar angle and latitude of the observer site are given, where polar angle + latitude = 90o. Representative cities near the indicated latitude are noted. Basically one sees here how the two pictures shown earlier are "interpolated" as we move from equator to north pole. The plane of any of the red curves shown makes an angle with the horizontal site tangent plane which is equal to the polar angle shown. We are just viewing the three red circles of Fig 3.1 or Fig 3.2 from various places on the earth. To obtain the figures for the southern hemisphere, the images below should be mirror-reflected left to right, that is, in the east-west-up plane. Then the red curves all tip to the right. Section 3: The Celestial Sphere 30 polar = 90 lat = 0 Singapore, Quito (equator) polar = 80 lat = 10 Caracas, Costa Rica polar = 70 lat = 20 Honolulu, Mexico City, San Juan polar = 60 lat = 30 Cairo, Houston polar = 50 lat = 40 New York, Salt Lake City polar = 40 lat = 50 Brussels Section 3: The Celestial Sphere 31 polar = 30 lat = 60 Helsinki polar = 25 lat = 65 Reykjavik polar = 20 lat = 70 Barrow polar = 10 lat = 80 Eureka Fig 3.12 polar = 0 lat = 90 North Pole Section 3: The Celestial Sphere 32 Maple Code Notes The code is based on equations (7.15) which ar e derived below in Section 7 (see notes there), x' = – R cos θt sin(ωt) y' = R [sinθt sinθ1 – cosθt cosθ1 cos(ω t) ] z' = R [sinθt cosθ1 + cosθt sinθ1 cos(ωt ) ] . (7.15) The functions are first entered with R = 1 and we then verify that x' 2+y'2+z'2 = 1, Next, the plot library is activated and constants are set in, The Maple code generates one plot at a time, and for this plot θ1 = 50o. The three red partial circles are generated in a display list p 1 : Three values of θ t are used here to make three red circles. The z part of the view option causes these circles to be plotted only for z' ≥ 0 which means only for daytime. The grey half sphere latitude guide circles are generated using equations (6.7) with R = 1, and then finally both the p1 and p2 plots are superpos ed, displayed and oriented by the last command, Section 3: The Celestial Sphere 33 (f) The Arctic and Antarctic Circles In going from the polar = 25 to the polar = 20 drawing in Fig 3.12, one notices that the winter solstice sun path has completely disappeared. In the following picture, we assume θt = -20o which is not quite at December solstice (which would be -23.4o). F i g 3 . 1 3 This picture is hard to draw because we really want the sun extremely far away. This is why the arrows which say "to sun" do not point to the close-in sun, but do point to the correct far-away sun. The bottom green line is the edge-on path of the close-in sun. The point here is that someone on the north pole cannot ever see the sun in this picture, no matter where it is on its green circular path. This is also true for any observer above the upper blue circle which is at θ 1 = 20o polar angle. In general, the condition for 24 hour night near the north pole is this θt < 0 AND θ1 < |θt| // north pole region always night The largest angle that θ 1 can be, with 24 hour da rkness above it, is 23.4o and this occurs at Dec solstice. The blue circle at this limiting point is called the Arc tic Circle (not shown). Below that circle it can never be 24 hour night. At the same time, an observer south of the mirror-image lower blue circle always sees the sun, as the lower arrow shows. The condition for 24 hour daylight here is θ t < 0 AND π- θ1 < |θt| south pole region always day The upper limit of the lower blue circle is called the Antarctic Circle. Above this circle it can never be 24 hour day. Section 3: The Celestial Sphere 34 Now consider the situation near June solstice F i g 3 . 1 4 The conditions now for all night and all day are θ t > 0 AND θ1 < |θt| north pole region always day θt > 0 AND π- θ1 < |θt| south pole region always night We can then summarize these four conditions: θ t < 0 AND θ1 < |θt| north pole region always night θt < 0 AND π- θ1 < |θt| south pole region always day θt > 0 AND θ1 < |θt| north pole region always day θt > 0 AND π- θ1 < |θt| south pole region always night (3.2) Section 4: Eccentricity effect on ψ(t) 35 4. How the Earth's Orbital Azimuthal Position ψ Varies with Time due to Eccentricity (a) A simple solution for ψ(t) Recall fro m Section 1 the equation of the e lliptical orbit of the earth around the sun, r(ψ ) = a(1-e2) / [1 + e cos ψ ] . (1.7) Using 2012 data, a = 149.598261 million km and eccentricity e = .0167. Here a is the semi-major axis of the ellipse, and the semi-minor axis will be called b below. Ignoring small external torques on the sun-earth system, the angular momentum L of the system is a conserved quantity. If m is the reduced mass of the system (m ≈ M Earth ) then Kepler's second law of orbital motion states that (see for example Ref [1] p 60 Eq (3-8) ) ψ• r2 = L / m = c o n s t a n t . ( 4 . 1 ) The area of a small triangular wedge of angle d ψ is dA = (1/2)r(rd ψ) = (1/2)r2ψ• dt = (L/2m) dt (4.2) and here we see (4.1) stated as the "equal areas swep t out in equal times" rule. Integrating both sides around the entire ellipse then gives A = ( L / 2 m ) T , ( 4 . 3 ) where T is the period of the orbit. But we know that A = πab for an ellipse, and also that b = a 1-e2 , so the area is A = π a21-e2 . Therefore, (L/m) = 2A/T = 2 πa 21-e2 / T . ( 4 . 4 ) We define Ω0 to be the mean angular velocity of the orbit, Ω0 ≡ (2π/T) // radians/time (4.5) and we then have (L/m) = Ω 0 a21-e2 . ( 4 . 6 ) Then we can then rewrite (4.1) as ψ• r 2 = Ω0a21-e2 . ( 4 . 7 ) Combining this with (1.7), r = a(1-e2) / [1 + e cos ψ ], we get Section 4: Eccentricity effect on ψ(t) 36 ψ• {a(1-e2) / [1 + e cos ψ ]}2 = Ω0a21-e2 or ψ• = Ω0 (1-e2)-3/2 [1 + e cos ψ ]2 or dψ/dt = C [1 + e cos ψ ]2 C = Ω0 (1-e2)-3/2 . (4.8) This equation can be integrated from t = 0 to t and ψ = 0 to ψ , dψ/[1 + e cos ψ ] 2 = C dt I ≡ ∫0 ψ dψ'/[1+ e cosψ ' ]2 = C t . // ψ(t=0) = 0 (4.9) Using the following integrals from p172 of Ref [2] with a = 1, b = e, A = 1, B = 0, and n = 2 : the integral I may be evaluated, I = (1-e 2)-1 { (2/ 1-e2 ) tan-1[1-e 1+e tan(ψ/2) ] – esinψ 1+ecosψ } = (1-e2)-3/2 { 2 tan-1[1-e 1+e tan(ψ/2) ] – e1-e2 sinψ 1+ecosψ } Since C = Ω0 (1-e2)-3/2 we then find that Ω0t = 2 tan-1[1-e 1+e tan(ψ/2)] – e1-e2 sinψ 1+ecosψ // radians (4.10) Later we shall use years as the unit of time, so then Ω0 = 2π radians/year and the time span of an orbit is from t = 0 at perihelion until t = 1 year at the next perihelion. Equation (4.10) appears in Ref [3] as (8). Looking at Fig 1.3, we can identify th ese angles with the solstices and equinoxes, Section 4: Eccentricity effect on ψ(t) 37 ψ = π/2-φper ≡ α March equinox // α was defined this way in (1.9) ψ = π-φper June solstice ψ = 3π/2 - φper Sept solstice ψ = 2π - φper Dec solstice Then using (4.10), we may compute the time (0-1) and day (0-365.25) at which these events occur : ( 4 . 1 1 ) In the case of a circular orbit, we expect to find that ψ = Ω 0t and that is what (4.10) gives with e → 0. In this case, keeping the "perihelion" at the same location, the time of the March equinox is determined by Ω 0tMe = ψ = π/2-φper ≡ α so that t Me = α/Ω0 = (π /2-φper)/2π / / e = 0 ( 4 . 1 2 ) and then the four times or interest are close to, but different from, those shown above, and they are separated from each other by 1/4 orbit : Section 4: Eccentricity effect on ψ(t) 38 (4.11) e=0 Equation (4.10) gives t(ψ ) which can be numerically inverted to obtain ψ(t). Since e is so small, it is much simpler to back up and use an approximation for ψ(t). Consider again (4.8), dψ/dt = C [1 + e cos ψ ] 2 (4.8) ≈ C [ 1 + 2 e cos ψ ] . ( 4 . 1 3 ) The fractional error in this approximation is ~e 2 = (.0167)2 = 3 x 10-4. Since the orbit is nearly circular, we approximate ψ = Ω0t inside the sine to get (another order e2 error contribution) dψ/dt ≈ C [ 1 + 2 e cos(Ω 0t ) ] . ( 4 . 1 4 ) Taking the mean of the left hand side over a year gives <d ψ/dt> = C, but in our approximation this must be C = Ω o, so defining Ω ≡ dψ/dt as the instantaneous angular veloc ity of the orbit at time t we get Ω(t) = Ωo [ 1 + 2 e cos(Ω 0t ) ] . ( 4 . 1 5 ) Here is a plot of Ω(t)/Ω 0 as in (4.15) which can be correlated with Fig 1.3. Starting at t = 0 (Jan 3 perihelion), Ω(t) first slows down, reaching a minimum va lue at aphelion, then speeds back up: Section 4: Eccentricity effect on ψ(t) 39 Fig 4.1 Alternatively, we can integrate (4.14) with C = Ωo to get an expression for ψ(t), dψ/dt = Ω 0 [ 1 + 2ecos( Ω0t ) ] (4.14) ψ(t) - ψ (0) = Ω0t + 2e Ω 0 ∫0 t dt cos(Ω 0t) = Ω0t + 2e Ω 0(1/Ω0) [sin(Ω0t)]|t 0 = Ω 0t + 2e sin( Ω0t ) . ( 4 . 1 6 ) If we agree to start our clock at the instant the earth is at perihelion (Jan 3) , then ψ(0) = 0 as indicated in Fig 1.3. We already assumed this in writing ψ = Ω0t inside the sine above. In this case, ψ(t) = Ω 0t + 2e sin( Ω0t) // radians (4.17) and using (1.9) and (1.10) we may express the angle φ' of Fig 1.3 this way, φ'(t) = Ω 0t + 2e sin( Ω0t) + φper - π/ 2 ( 4 . 1 8 ) where Ω o = 2π rad/year. During a year the sine phase goes ar ound once with a zero average value, so that <ψ(t) – Ω 0t > = 0 . This plot of ψ(t) - Ω 0t shows how ψ(t) in degrees deviates from straight-line clock time for t = (0,1) yr: Section 4: Eccentricity effect on ψ(t) 40 Fig 4.2 The above curve will reappear in Section 6 (e) in th e discussion of the Equation of Time for the "black earth". The preceding discussion addresses the variation in ψ due to eccentricity in the orbit of the earth around the sun. As shown in Fig 3.4, this same ψ applies to the projected s un's orbit around the ecliptic circle on the celestial sphere. However, when this red circle gets projected down onto the equator of the celestial sphere (Fig 3.5), a "tilt effect" arises which fu rther affects azimuthal timing on that equator. This will be clarified in Section 5 below. (b) A systematic solution for ψ(t) and estimate of error in using (4.17) The exact (non-linear) differential equation for ψ shown in (4.8) is this d ψ/dt = C [1 + e cos ψ ] 2 = C [1 + 2ecos ψ + e2cos2ψ ] . (4.8) One can expand ψ in a power series in smallness parameter e ≈ .0167 , ψ (t) ≈ ψ (0)(t) + e ψ(1)(t) + e2 ψ(2)(t) + .... (4.19) where we require that for each order ψ (n)(0) = 0, so that ψ(0) = 0. Now using cos(x+ ε) ≈ cos(x) – ε sin(x) + O( ε2) ( 4 . 2 0 ) we have cosψ = cos[ψ (0)(t) + e ψ(1)(t)] = cos[ ψ(0)(t)] - e ψ(1)(t) sin[ψ(0)(t)] + O(e2) cos 2ψ = cos2[ψ(0)( t ) ] + O ( e ) . ( 4 . 2 1 ) Defining c ≡ cos[ψ (0)(t)] and s ≡ sin[ψ(0)(t)] this pair of equations becomes Section 4: Eccentricity effect on ψ(t) 41 cosψ = c - e s ψ(1)(t) + O(e2) cos 2ψ = c2 + O ( e ) . ( 4 . 2 2 ) Retaining only terms through order e2, the right side of (4.8) becomes = C [1+2e{ c - e s ψ(1)(t)} + e2c2] = C [ 1 +2ec - 2e 2s ψ(1)(t) + e2c2 ] = C + e [2Cc] + e 2[C(c2-2s ψ(1)( t ) ] . ( 4 . 2 3 ) Balancing powers of e with the left side of (4.8) give n by differentiating (4.19), we get then this set of equations, dψ (0)/dt = C dψ (1)/dt = 2Cc = 2C cos[ ψ(0)(t)] dψ (2)/dt = Ccos2[ψ(0)(t)]- 2Csin[ ψ(0)(t)] ψ(1)(t) . (4.24) The solution to the first equation is just ψ (0)(t) = Ct and therefore C = Ω 0 = 2π radians/year, ψ (0)(t) = Ω0t . // solution for a perfectly circular orbit. e = 0 (4.25) The remaining two equations are then dψ (1)/dt = 2Cc = 2 Ω0cos(Ω0t) dψ (2)/dt = Ω0cos2(Ω0t)- 2 Ω0sin(Ω0t) ψ(1)(t) (4.26) The ψ(1) equation can be trivially integrated to give ψ(1)(t) = 2 sin( Ω0t ) ( 4 . 2 7 ) so the third equation becomes dψ (2)/dt = Ω0[ cos2(Ω0t)- 4sin2(Ω0t) ] = Ω 0[ (1/2)(1+cos(2 Ω0t)) - 4(1/2)(1 – cos(2 Ω0t)) ] = ( 1 / 2 ) Ω 0[ (1+cos(2 Ω0t)) - 4 (1 – cos(2 Ω0t)) ] Section 4: Eccentricity effect on ψ(t) 42 = ( 1 / 2 ) Ω0[ (-3+5cos(2 Ω0t ) ] . ( 4 . 2 8 ) The integral of interest in determining ψ(2) is ∫0 t [-3+5cos(2 Ω0t') ] dt' = -3t + 5/(2 Ω0)sin(2Ω0t ) ( 4 . 2 9 ) so that ψ (2)(t) = –(3/2) Ω0t + (5/4) sin(2 Ω0t ) . ( 4 . 3 0 ) Therefore the solution ψ(t) through second order in eccentricity e is ψ(t) = Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) + order(e3) . (4.31) The e2 contribution is very small, being on the order of (5/4)e2 = (5/4) (.0167)2 = 3.5 x 10–4 . Converting this from radians to degrees makes this term on the order of 1/50th of one degree. Here is a plot of ψ(t) - Ω0t in degrees, with (red) and without (black) the e2 term, Fig 4.3 For our purposes in this document, such as plotti ng analemmas, there is no need to keep the e2 term and we can use this simple approximation for ψ(t) ψ(t) = Ω 0t + 2e sin( Ω0t) which is the same as (4.17) obtained by our ad hoc method. The result (4.31) is sometimes written this way, δ ≡ ψ(t) – Ω 0t ≈ 2e sin(Ω0t) + (5/4)e2 sin(2Ω0t) , (4.32) Section 4: Eccentricity effect on ψ(t) 43 where δ is then the difference between the true anomaly ψ and the mean anomaly Ω0t as discussed at the end of Section 1. This expression appears in the code shown in Ref [4] where the time scale is chosen so that a year is 2 π units long, so t' = 2 πt = Ω0t. That code reads (kappa = e) As a check on (4.31), let's use it to verify our earlier calculation of t Me shown in (4.11) by setting t = t Me in (4.31): Thus, when we install t = t Me = .2084458064 yrs from (4.11), we recover our ψMe = α as shown at the end of (4.11) to 6 decimal places accuracy. One could go on and find the e3 and higher terms in the expansion (4.31) by the method outlined above. Section 5: How tilt angle θt varies in time 44 5. How the Earth's Tilt Angle θt Varies with Time In this section we think of θt as the angle of tilt of the earth's rotation axis toward or away from the sun. This angle varies between - θtmax and θtmax where θtmax is the earth's obliquity 23.4o. This tilt θt is described below in a special frame of reference that rotates around the sun with the earth so that the sun and earth maintain their positions (earth on left, s un on right). Later in Section 6 (b) we show that θt is also the declination of the sun on the celestial sphere , its more conventional interpretation. This fact is implied by the inset in Fig 5.1 below. We first replicate Fig 1.3 since it is closely related to the drawing we are about to present: Fig 1.3 In Fig 1.3 angle φ' denotes the azimuthal position of the earth with respect to the March equinox. This angle φ' appears in the following picture, along with angle φ~ = π-φ' and θ~ = θtmax : F i g 5 . 1 The earth is represented by the black origin dot on the le ft, and its equatorial plane by the gray disk. The red arrow ω is the angular momentum vector of the rotating earth. This vector's tip precesses clockwise Section 5: How tilt angle θt varies in time 45 around the pink circle as the year progresses. The sun is on the right, "looking at the earth". The reader is invited to correlate this drawing and its ma rkings with Fig 1.3. For example, when φ' = 0, in Fig 1.3 the sun sees the earth tilting to the right. The red vector on th e cone in Fig 5.1 also ti lts to the right as viewed from the sun when φ ' = 0 (March equinox). For the situation shown in Fig 5.1, we want to determine the location of the sun's path projected on the earth. We know that this path is a circle of cons tant latitude because the path is created by the earth rotating (see blue circles in Fig 3.1). We only need to find one point on this circle. The line from sun to earth-center penetrates the earth's surface at a point wh ich must lie on this blue latitude circle of interest, since the circle is a projection of the sun's path (see Fig 5.1 inset in which ω^ is tilted toward the viewer as on the left). First consider this standard expression for a spheri cal-coordinates radial unit vector in terms of the Cartesian basis vectors, r^ = sinθ~cosφ~ x^ + sinθ~sinφ~ y^ + cosθ~ z^ . ( 5 . 1 ) We can think of r^ = ω^, the direction of the red arrow. We know that ω^ is a unit vector normal to the equatorial plane of the earth. We know that y^ points toward the sun. Therefore, the colatitude θ of our penetration point of interest is given by cos θ = ω^ • y^ as shown in the Fig 5.1 inset . This colatitude θ is just the polar angle measured down from the north pole of the tilted earth to the y^ vector. The corresponding latitude we shall call θ t = π/2 - θ . Now, cosθ = cos(π /2-θt) = sinθt . ( 5 . 2 ) Therefore, sinθt = cosθ = ω^ • y^ = sinθ~sinφ~ = sinθtmax sin(π-φ') = sinθtmax sinφ' (5.3) and we obtain the key result which is sinθ t(t) = sinθtmax sinφ'(t) where θtmax = 23o. (5.4) To check the result, at June solstice φ' = π /2 so sinθ t = sinθtmax* 1 => θt = θtmax . Things check similarly at Dec solstice and the equinoxes. From Fig 3.4, φ' = ψ - α , so using ψ from (4.31) we get φ'(t) = ψ - α = Ω 0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) - α // α = π/2 - φper so (5.4) becomes θ t(t) = sin-1 { sinθtmax sin(ψ-α)} = s i n -1 { sinθtmax sin[Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) - α] } (5.5) where t = 0 at perihelion Jan 3. The following black plot of θ t(t) shows (5.5) with only its first term Ω0t, while the red plot shows (5.5) with the first two terms. The e2 term makes no visible difference so we don't bother with it. Section 5: How tilt angle θt varies in time 46 Fig 5.2 The extremes of the plot are θt ≈ ± 23o at the solstices, while the ze ro crossings are the equinoxes θt= 0. At t = 0 Jan 3 (left edge) we are just after December solstice where the sun was at - 23o. A geometric interpretation of θ t. It is possible to interpret θ t as "the amount of tilt of the earth toward the sun" according to the following somewhat elaborate construction : F i g 5 . 3 The distance from the z axis to the tip of the red unit vector is sin θ~ as indicated. This distance can then be projected onto the gray wall (the yz plane) to get a distance sin θ~sinφ~. We then create a black unit vector on the gray wall which has this distance sin θ~sinφ~ as its y component. If we denote the angle of this black unit vector relative to the z axis as α, then sin α = sinθ~sinφ~. But in (5.3) we had sin θ~sinφ~ = sinθt. Therefore, α = θt and then we can interpret θt as the "amount of tilt of the earth toward the sun". Note that the black arrow is not the projection of the red one onto the wall, because the projected arrow would be too short to be a unit vector. Section 6: Coordinates of the Sun on the Celestial Sphere 47 6. Coordinates of the Sun on the Celestial Sphere (a) Cartesian coordinates x,y,z and velocity of the sun in Frame S Consider this picture showing the tilted red ecliptic circle and the black equatorial circle on the celestial sphere (compare to Fig 3.5 above). Fig 6.1 The black equatorial circle is associated with Frame S with x,y,z axes, while the ecliptic red circle is associated with Frame S' with x',y',z' axes. Staring at the figure one sees that z^' = Rx(θtmax) z^ and x^' = x^ => x' = x y^' = Rx(θtmax) y^ ( 6 . 1 ) where θtmax = 23o. According to the rule that vectors rotate backwards from the unit vectors (see Appendix B), we then have in 2D that r = Rx(θtmax)r' => ⎝⎛ ⎠⎞ y z = ⎝⎛ ⎠⎞ cosθtmax -sinθtmax sinθtmax cosθtmax ⎝⎛ ⎠⎞ y' z' => y = cos θtmaxy' – sinθtmaxz' z = sin θtmaxy' + cosθtmaxz ' . ( 6 . 2 ) When θtmax = 0 this is clearly valid. When θ tmax = π/2 we get y = -z' and z = y' which agrees with the picture -- just a validity check on the sign of the sine. Now our interest is not so much in arbitrary point s, but in points which lie on the red circle. For these points we have z' = 0 so the above becomes Section 6: Coordinates of the Sun on the Celestial Sphere 48 x = x' y = cos θtmax y' z = sin θtmax y ' . ( 6 . 3 ) In Frame S' we can define spherical coordinates r', θ',φ' related to x',y,'z' according to the usual equations, x' = r'sin θ'cosφ' y' = r'sin θ'sinφ' z' = r'cos θ' . But again, for points on the red circle, θ' = π/2 and r' = R so x' = Rcos φ' y' = Rsin φ' z' = 0 . Inserting these expressions into (6.3) then gives x = Rcos φ' y = Rcos θ tmax sinφ' z = Rsin θtmax sinφ' . ( 6 . 4 ) Equations (6.4) give the Cartesian coordinat es of the sun in Frame S. Recall that angle φ' appearing in (6.4) is given by, φ' = ψ - α (1.10) ψ = Ω 0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) // + order e3 which we ignore (4.31) α = π/2 - φper (1.9) (6.5) where t = 0 is perihelion (Jan 3) and t = 1 is the next perihelion. The velocity and acceleration of the sun's motion on the celestial sphere (as seen by our Frame S observer) can be obtained by differentiating (6.4). For example, x• = –R φ•' sinφ' y• = R φ•' cosθ tmax cosφ' z• = R φ•' sinθtmax cosφ' ( 6 . 6 ) where from (4.18) φ•' = ψ• = Ω0 [1 + 2e cos( Ω0t) + (5/2) cos(2 Ω0t)] . (b) Spherical coordinates θ, φ of the sun in Frame S In Frame S, for poi nts on the celestial sphere, we can define spherical coordinates R, θ,φ related to x,y,z according to the usual relations, x = Rsin θcosφ y = Rsin θsinφ z = Rcos θ . ( 6 . 7 ) Section 6: Coordinates of the Sun on the Celestial Sphere 49 Comparison of equations (6.7) with (6.4) shows that, sinθcosφ = cosφ' sinθsinφ = cosθ tmax sinφ' cosθ = sinθtmax sinφ' . ( 6 . 8 ) Dividing the second equation by the first gives tanφ = cosθ tmax tanφ' . ( 6 . 9 ) We then have from the last line in (6.8) and (6.9), θ(t) = cos -1[sinθtmax sinφ'(t)] φ(t) = tan-1[cosθtmax tanφ' ( t ) ] ( 6 . 1 0 ) r(t) = R where φ' given by (6.5). Equation (6.10) gives the official Frame S spherical coordinates of the projection of the sun on the celestial sphere. We now claim that our polar angle θ(t) is related to the tilt angle θ t(t) introduced in Section 5 (which was the angle of tilt of the earth toward the sun at some time t) by, θ(t) = π/2 - θt(t). => sin θt = cosθ . ( 6 . 1 1 ) If this were true, the last line of (6.8) would say sinθ t = sinθtmax sinφ' . But this is the same as (5.4). Ther efore (6.11) is correct and, since θ is the polar angle or colatitude of the sun on the celestial sphere, θt must be the latitude or "declination" of the sun on the celestial sphere, as shown in Fig 6.1. Equations (6.10) can then be wr itten as (the first equation then agrees with (5.5)) θ t(t) = sin-1[sinθtmax sinφ'(t)] // latitude (declination) of sun on celestial sphere φ(t) = tan-1[cosθtmax tanφ'(t)] // longitude (right ascension) of sun on celestial sphere (6.12) where the sun azimuth φ', ψ and α as given by (6.5). Thus, θ t(t) = sin-1[sinθtmax sin(ψ-α)] φ(t) = tan-1[cosθtmaxtan(ψ-α) ] . ψ = Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) // + order e3 which we ignore (4.31) α = π/2 - φper (1.9) (6.13) These equations gives the declination θ t and right ascension φ of the sun on the celestial sphere during the year t = (0,1) where t = 0 is perihelion, Jan 3. Section 6: Coordinates of the Sun on the Celestial Sphere 50 It is useful to think of a fictitious sun which orbits the earth in a perfectly circular orbit (e=0) and whose ecliptic lies in the equatorial plane of the celestial sphere ( θtmax = 0). For this sun the above equations are θt,f(t) = 0 φf(t) = tan-1[tan(ψf - α)] = ψf - α = Ω0t - α ψf = Ω0t α = π/2 - φper // fictitious sun (1.9) (6.14) Interpretation of the second equation of (6.12): The equation in question and its time derivative are, φ(t) = tan -1[cosθtmax tanφ' ( t ) ] (6.12) φ•(t) = { cos θtmax / [cos2φ' + cos2θtmaxsin2φ']} φ•' . ( 6 . 1 2 ) dot Near the equinoxes φ' = 0 and π one has cos2φ' ≈ 1 and sin2φ' ≈ 0 so then φ• ≈ cosθtmax φ•'. The angular velocity φ• at the equinoxes is slower than φ•' by factor cosθ tmax . The simple geometric reason for this, as Fig 6.1 shows, is that the two circular paths intersect at angle θtmax at the equinoxes, and the projection of speed on the ecliptic to speed on the eq uatorial circle is reduced by this factor at the equinoxes. Near the solstices φ' = π /2 and 3π /2 one has cos2φ' ≈ 0 and sin2φ' ≈ 1 so then φ• ≈ φ•'/ cosθtmax . At these two points φ• runs faster than φ•' by factor 1/cos θtmax . The geometric reason for this is slightly less obvious: at solstice an "orange slice" of the celestial sphere is thinner at the ecliptic than at the equator so the projection of a speed on the ecliptic gets magnified by factor 1/ cos θtmax . Both these situations are illustrated on the left of the following figure, where the red arrows indicate speed along the ecliptic and green arrows speed along the equatorial circle. The first graph shows crudely how the ratio of angular speeds must therefore vary during the year where τ = 0 is March equinox (Me). In the case that φ•' = constant, φ must be generally as shown on the far right, since the derivative of this graph gives the first graph minus 1. This situation applies to "the blue earth" below where e = 0 so φ•' = Ω 0. The actual plot is the blue curve in Fig 6.5 below, which begins at t = 0 which is the Jan 3 perihelion. Thus we have a physical interpretation of "the tilt effect" on azimuthal velocity of the sun. More generally , this interpretation explains why the full equation of time shown as the red curve in Fig 6.5 has two "cycles" during the year . The lower amplitude cycle and the higher amplitude cycle create the small and large lobes of the analemma in Section 11. Roughly the difference between the red curve and the blue curve in Fig 6.5 is the black curve which accounts for "the eccentricity effect" on azimuthal speed as studied in Section 4. Section 6: Coordinates of the Sun on the Celestial Sphere 51 (c) The Equation of Time The difference between the azi muth (right ascension) of the sun and that of the fictitious sun is known as "the Equation of Time" and is given by (6.13) and (6.14) above (using ψ and α as in (6.13) above) Δφ(t) ≡ φ - φf = tan-1[cosθtmax tan(ψ-α)] – (Ω0t-α) (6.15a) = tan -1[cosθtmax tan(ψ-α)] – tan-1[ tan(Ω0t-α)] (6.15b) = tan-1 [cosθtmax tan(ψ-α)– tan(Ω0t-α) 1 + cosθtmax tan(ψ-α) tan(Ω0t-α) ] (6.15c) where the last line follows from this identity from p 121 of Ref [5] . (6.16) There are other variations of (6.15). One can multiply top and bottom by cos( ψ-α)cos(Ω0t-α) to get Δφ(t) = tan-1 [– sin(Ω0t-α) cos(ψ-α) + cosθtmax sin(ψ-α) cos(Ω0t-α) cos(ψ -α)cos(Ω0t-α) + cosθtmax sin(ψ-α) sin(Ω0t-α) ] . (6.15d) A somewhat peculiar form is given in equation (21) of Ref [3] as –Δφ(t) = tan -1[–sin(ψ-Ω0t) + tan2(θtmax/2) sin(ψ + Ω0t -2α) cos(ψ -Ω0t) + tan2(θtmax/2) cos(ψ + Ω0t -2α) ] ≡ μ(t) . (6.15e) To verify (6.15e), comparing to (6.15c) we have to show that (-1) cosθtmax tan(ψ-α)– tan(Ω0t-α) 1 + cosθtmax tan(ψ-α) tan(Ω0t-α) = –sin(ψ-Ω0t) + tan2(θtmax/2) sin(ψ + Ω0t -2α) cos(ψ -Ω0t) + tan2(θtmax/2) cos(ψ + Ω0t -2α) ? Let θtmax = x, ψ = y, α = a, Ω0t = c, then we have to show that cosx tan(y-a)– tan(c-a) 1 + cosx tan(y-a) tan(c-a) = (-1) –sin(y-c) + tan2(x/2) sin(y + c -2a) cos(y-c) + tan2(x/2) cos(y + c-2a) ? A B = – C D ? AD + BC = 0 ? This is a job for Maple, and the conclusion is that we can remove all the question marks above : Section 6: Coordinates of the Sun on the Celestial Sphere 52 Wiki provides a nice plot of the equation of time (see Ref [4] or Fig 6.8 below) which is based on the form (6.15e) above. Here is some of that code, where Ω0 = 1 (year = 2 π time units) , δ ≡ ψ-t, ε ≡ θtmax, κ ≡ e. This δ is the same as our (4.32). [ German Zähler = numerator, Nenner = denominator ] Fig 6.2 (d) Plots of φ (right ascension of the sun) and Δφ (equation of t ime) Note : For all plots in this section and the next, the left edge of the plot is t = 0, perihelion, Jan 3. Defining q and s by q ≡ 2e sin(Ω 0t) + (5/4)e2 sin(2Ω0t) // ψ - Ω0t s ≡ Ω0t + φper - π/2 // Ω0t - α ; s+q = ψ-α (6.17) we shall make four vers ions of equation (6.13), φ(t) = tan -1[cosθtmaxtan(ψ - α)], and we associate a color with each version (note that φred = φ = right ascension of the actual sun) φ red(t) = tan-1[cosθtmax tan(s+q)] // the actual sun φblue(t) = tan-1[cosθtmax tan(s)] // e = 0 so blue means tilt only φblack (t) = tan-1[tan(s+q)] = s+q // θtmax = 0 so black is eccentricity only φgreen (t) = tan-1[tan(s)] = s = φ f(t) // e = 0 AND θtmax = 0 (green is fictitious sun) (6.18) Our first task is to enter the four functions with 2012 parameters, Section 6: Coordinates of the Sun on the Celestial Sphere 53 The function arctan2Pis is a special version of tan-1 which causes the arctangent to be a continuous function even if it goes out of the nominal range (0,2 π). It keeps internal track of a certain winding number as discussed in Appendix A. The scale of the data is radians. The full ψ function is used in q, even though the e2 term really makes no significant difference (see Section 4). Here then is a plot of the four functions, radians F i g 6 . 3 It is far more interesting to plot the e quations of time which are these differences, Δφ red(t) ≡ φred(t) - φgreen (t) Δφblue(t) ≡ φblue(t) - φgreen (t) // e = 0 so only tilt is active Δφblack (t) ≡ φblack (t) - φgreen (t) // θtmax = 0 so only eccentricity is active Δφgreen (t) ≡ φgreen (t) - φgreen (t) = 0 (6.19) Δφpink(t) ≡ Δφblue(t) + Δφblack (t) = φblue(t) + φblack (t) – 2 φgreen (t) (6.20) Section 6: Coordinates of the Sun on the Celestial Sphere 54 where recall φgreen is for the fictitious sun: radians F i g 6 . 4 Comment : The same plot is obtained with this single command (but see end of Appendix A) Looking at the vertical scale (radians), one sees that these Δφ differences are all quite small compared to the previous plot. Since the pink curve visibly devi ates somewhat from the red curve, we see that the often assumed fact that Δφ blue(t) + Δφblack (t) = Δφred(t) is only approximately true, and this is due to the relative smallness of angle θtmax. But qualitatively, we see that the red curve is roughly the sum of the blue curve and the black curve, so the effects of orbital tilt (blue curve) and eccentricity (black curve) are roughly additive. We shall dispense with the pink curve from now on. By adding a factor of 180/ π we can replot in degrees, Section 6: Coordinates of the Sun on the Celestial Sphere 55 degrees F i g 6 . 5 Since 15 degrees of azimuth corresponds to one hour of earth rotation time (1 degree = 4 minutes), one often sees the equation of time plotted with minutes as the vertical scale ( factor now 4(180/ π) ) minutes F i g 6 . 6 In these plots the left edge t = 0 corresponds to perihelion or Jan 3. The actual expressions for the three Δφ differences are Δφ red(t) ≡ φred(t) - φgreen (t) = tan-1[cosθtmax tan(s+q)] – s Δφblue(t) ≡ φblue(t) - φgreen (t) = tan-1[cosθtmax tan(s)] – s Δφblack (t) ≡ φblack (t) - φgreen (t) = (s+q) - s = q (6.21) q ≡ 2e sin(Ω 0t) + (5/4)e2 sin(2Ω0t) // = ψ - Ω0t s ≡ Ω0t + φper - π/2 // = Ω0t - α ; s+q = ψ-α (6.17) Section 6: Coordinates of the Sun on the Celestial Sphere 56 It seems clear that the mean value of Δφblack (t) = 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) will be zero over the year since it is an odd function of t (think of the year as running t = -1/2 to +1/2) and is periodic. The mean of Δφblue(t) is also zero, but that fact is less obvious. One can write, using the trig identity (6.16) quoted above, Δφblue(t) = tan-1[cosθtmax tan(s)] - tan-1[tan(s)] = tan-1 [(cosθtmax – 1)tan(s) 1 + cosθtmaxtan2(s) ] , (6.22) where s = Ω0t - α = Ω0t', where t' is a shifted time. Then the function is odd in t' and also periodic in t', so the blue mean is in fact 0. The mean of Δφ pink(t) ≡ Δφblue(t) + Δφblack (t) is then exactly 0 as well. The above argument fails on the red mean, Δφred(t) = tan-1[cosθtmax tan(s+q)] - tan-1[tan(s)] = tan-1 [cosθtmaxtan(s+q) – tan(s) 1 + cosθtmaxtan(s)tan(s+q) ] . (6.23) Nevertheless, the red mean is very small and we can regard it as being zero to a few parts in a million. Here are the means computed by Maple by approxima ting each integral as a sum of diy = 365 rectangles (Maple cannot integrate a function which uses arctan2Pis) Now replacing φ red by its original symbol φ, and with φ f = Ω0t - α as in (6.14), we have red mean = 0 = < Δφ red > = <Δφ> = <φ > - <φf> so <φ f> = <φ> . This last equation says that during the year, the mean va lue of the azimuth of the fictitious sun is the same as the mean value of the azimuth of the actual sun. For this reason, the fictitious sun is often referred to as the mean fictitious sun . This fictitious sun is represented by the horizontal green line in the equation of time plots above, and the distance between the red curve and the green line is then the azimuthal distance between the actual sun and the mean fictitious sun at any instant of time t. As the plots show, the fictitious sun's right ascension matches that of the actual sun exactly four times per year (places where Δφ = 0). Section 6: Coordinates of the Sun on the Celestial Sphere 57 Frequently the equation of time plot is shown in minutes with the vertical axis flipped (- 4 (180/ π) ), minutes Fig 6.7 As a check on our plotting methods, if we temporarily set our orbital parameters to those quoted above in Ref [4] , e = .016722 θ tmax = 23.45 φper = 11.5o α = π/2-φper = 78.5o we can replot the flipped Fig 6.7 and then supe rpose that on this wiki plot of Ref [4], F i g 6 . 8 Section 6: Coordinates of the Sun on the Celestial Sphere 58 Here is the superposed result after axes ar e adjusted to match in both directions, F i g 6 . 9 and the plots are basically identical. Their parame ters seem to match the year 1919 according to the NASA source http://aom.giss.nasa.gov/srorbpar.html , Table 6.1 In Section 8 we shall investigate the significance of these Δφ plots in terms of sundials and analemmas. Section 6: Coordinates of the Sun on the Celestial Sphere 59 (e) Plots of θt (declination of the sun) Using the same q and s shown in (6.17) q ≡ 2e sin(Ω 0t) + (5/4)e2 sin(2Ω0t) // ψ - Ω0t s ≡ Ω0t + φper - π/2 // Ω0t - α ; s+q = ψ-α (6.17) we can write the four "color" versions of the θ t equation in (6.13) as θ t,red (t) = sin-1[sinθtmax sin(s+q)] // the actual sun θt,blue (t) = sin-1[sinθtmax sin(s)] // e = 0 so blue means tilt only θt,black (t) = 0 // θtmax = 0 so black is eccentricity only θt,green (t) = 0 // e = 0 AND θtmax = 0 (fictitious sun) (6.24) Having Maple plot the red and blue functi ons in the same manner as done above, we find that θt,red and θt,blue are nearly identical over the year, with the red being slightly above the blue near t = 75 and 275 days, degrees Fig 6.10 Section 6: Coordinates of the Sun on the Celestial Sphere 60 We can plot θt,red - θt,blue to see the small effect that eccentricity adds to the tilt-only curve: degrees F i g 6 . 1 1 Since θ t,blue (t) is an odd function of t and is periodic in t, we get < θt,blue (t)> = 0. Since the difference plotted above is mostly positive, it must be that the mean of < θt,red (t)> is positive. Computing the means we find Section 7: Coordinates of the Sun viewed from Frame S' 61 7. Coordinates of the Sun in an Arbitrary Frame S' Glued to the Earth (a) Setup: θ 1, φ1, φb, φb0, τ, tMe Frame S of Fig 6.1 continues as Frame S in this new dr awing, a frame fixed relative to the celestial sphere with e1 = x^, e2 = y^, e3 = z^. We discard Frame S' of Fig 6.1 and re place it with the one shown here, a new Frame S' glued to the rotating earth ( "east" basis vector e '1 points mostly away from the viewer), Fig 7.1 Star 1 is some star one sees looking from the earth through the March equinox point on the celestial sphere. In 2012 this happens to be somewhere in th e Pisces constellation (but it is still called "the First Point of Aries" E, where it was 2000 years ago). Recall (Section 3(d)) that the equinox moves 1/72 degree along the ecliptic each year due to astronomical precession, causing it to also move on the celestial equator, so a new "Star 1" has to be found each year. Thus, Frame S very slowly rotates, a fact we ignore below. A system of spherical coordinates for Frame S has these unit basis vectors at the observer's location b on the earth's surface, r^ = e' 3 θ^ = - e'2 φ^ = e'1 , (7.1) where we associate each of these unit basis vectors with a Frame S' unit basis vector. Again, the earth is rotating counterclockwise in the above figure which is a snapshot at time t. The red Frame S' rotates with the earth, the black Fram e S is essentially fixed relative to the stars. Section 7: Coordinates of the Sun viewed from Frame S' 62 The observation site is placed on the earth at polar coordinates ( θ1, φ1) where as usual θ1 is measured down from the north pole (colatitude), and φ1 is the longitude measured east of the Greenwich meridian. It is convenient now to defi ne a shifted time coordinate τ ≡ t - tMe where t Me is the time from perihelion to March equinox as computed in (4.11). Then τ = 0 at March equinox which is the basic azimuthal reference time in Fig 7.1. This distinction is very important: t = time ranging 0 to 1 year where t = 0 at perihelion τ = time ranging 0 to 1 year where τ = 0 at March equinox τ = t - t Me . ( 7 . 2 ) The time t was convenient for the work done in Section 4 concerning orbital eccentricity, while τ is convenient in our current context because τ = 0 when the sun is at 0o right ascension at March equinox, see Fig 3.5. Let φb0 be the right ascension that vector b had at τ = 0. Then at any later time, φ b = φb0 + ωsτ = φbo + ωst – ωs tMe. ( 7 . 3 ) We use the sidereal rate ω s because φ b is relative to the celestial sphere which is at rest with respect to the stars. We use + ωs because the earth is rotating counterclockwise causing φb to increase. Therefore, the location of the observation point is given in Frame S spherical coordinates by b (τ) = (RE,θ1,φb(τ)) . (7.4) The sun in Fig 7.1 is taken to be the projection of the sun on the celestial sphere. This sphere has a very large radius R >> r S-E >>> RE . The position of this sun projection in Frame S spherical coordinates is given by r = (R, θ, φ) = (R, π/2- θ t, φ) ( 7 . 5 ) where the two angles φ and θ t are given by (6.13) and appear in Figs 6.1 and 7.1. Here φ is the right ascension and θt is the declination of the sun on the celestial sphere. In Cartesian coordinates the sun position (7.5) is x = R sin( π/2-θ t) cosφ = R cos θt cosφ y = R sin( π/2-θt) sinφ = R cos θt sinφ z = R cos(π /2-θt) = R sin θt ( 7 . 6 ) r = x x^ + y y^ + z z^ = x e1 + y e2 + z e3 where we have used sin( π/2-β) = cosβ and cos( π/2-β) = sinβ. Equations (7.6) are the same as equations (6.7) since sin θ = cosθt and cosθ = sinθt. Section 7: Coordinates of the Sun viewed from Frame S' 63 (b) Obtaining a value for φb0 in terms of φ G0 Let φG(τ) be the right ascension of the projected Greenwich meridian at time τ. We know that at any time φ b(τ) – φG(τ) = φ1 where it is understood that everything is mod 2 π. Therefore φ b(0) – φG(0) = φ1 or φ b0 = φ1 + φG0 . ( 7 . 7 ) Here φ G0 is the right ascension of Greenwich on the last March equinox. This can be obtained from φ G0 = [ 12:00 – (Greenwich mean time at March equinox) ] * 150 . (7.8) For example, if at the March equinox (our τ = 0) the Greenwich clock read exactly noon = 12:00, then the Greenwich meridian was exactly facing the sun so its projection onto the celestial sphere had right ascension 0o. If the Greenwich clock instead read 11:00 = 11AM at equinox, then the projection of the Greenwich meridian onto the celestial sphere lay at right ascension + 15o at τ = 0. Here is some data from http://aom.giss.nasa.gov/srver4x3.html , T a b l e 7 . 1 For example, in 2012 the March equinox occurred at 05:20 GMT so we find φ G0 = [ 12:00 – 05:20 ] * 150 = (12 - 5.33) * 150 = 100.05o = 1.75 radians Section 7: Coordinates of the Sun viewed from Frame S' 64 and then φb0 = φ1 + φG0 gives the constant φb0 which appears in the various equations below. (c) Cartesian coordinates x',y',z' and velocity of the sun in Frame S' We wish to know r ' , the position of the celestial sphere sun as viewed from rotating Frame S' . Our general relation is (see Fig 7.1) r = r' + b . Since R >> b = R E we have r' ≈ r . ( 7 . 9 ) As shown in Ref [6] Appendix A (A.13.c), at some arbitrary point r, θ,φ in Frame S the Cartesian unit vectors are given in terms of the spherical unit vectors according to x^ = cosφsinθ r^ + cosφcosθ θ^ - sinφ φ^ y^ = sinφsinθ r^ + sinφ cosθ θ^ + cosφ φ^ z^ = cosθ r^ – sinθ θ^ . (7.10) We wish to apply these equations at our observation point b = (R E,θ1,φb), and we wish to convert to the Frame S' basis vectors e 'n, so we make these edits (using 7.1), θ → θ1 φ → φb r^ → e'3 θ^ → - e '2 φ^ → e'1 to get x^ = cosφ b sinθ1 e'3 – cosφb cosθ1 e'2 – sinφb e'1 y^ = sinφb sinθ1 e'3 – sinφb cosθ1 e'2 + cosφb e'1 z^ = cosθ1 e'3 + sinθ1 e'2 . ( 7 . 1 1 ) Therefore r' ≈ r = x x^ + y y^ + z z^ = x (cos φ bsinθ1 e'3 – cosφbcosθ1 e'2 – sinφb e'1) + y ( s i n φbsinθ1 e'3 – sinφbcosθ1 e'2 + cosφb e'1) + z ( c o s θ1 e'3 + sinθ1 e'2) = [ – x s i n φ b + y cosφb]e'1 + [ – x cosφbcosθ1 – ysinφbcosθ1 + z sinθ1]e'2 + [ x c o s φbsinθ1 + y sinφbsinθ1 + z cosθ 1]e'3 ( 7 . 1 2 ) Section 7: Coordinates of the Sun viewed from Frame S' 65 = x ' e'1 + y' e'2 + z' e'3 so that the Cartesian sun location c oordinates in Frame S' are given by x' = – x sinφ b + y cosφb y' = – x cosφbcosθ1 – ysinφbcosθ1 + z sinθ1 z' = x cos φbsinθ1 + y sinφbsinθ1+ z cosθ 1 . ( 7 . 1 3 ) Next, into these three expressions we insert our previous results (7.6), x = R cos θ t cosφ y = R cos θt sinφ z = R sin θt (7.6) to get: : x' = - x sin φ b + y cosφb = - R c o s θt cosφ sinφb + R cosθt sinφ cosφb = R c o s θt[sinφ cosφb - cosφ sinφb] = R c o s θt sin(φ –φb) y' = – x cosφ b cosθ1 – y sinφb cosθ1 + z sinθ1 = – R cosθt cosφ cosφb cosθ1 – R cosθt sinφ sinφb cosθ1 + R sinθt sinθ1 = R [ s i n θt sinθ1 – cosθtcosθ1(cosφ cosφb + sinφ sinφb) ] = R [ s i n θt sinθ1 – cosθt cosθ1 cos(φ –φb) ] z' = x cos φ b sinθ1 + y sinφb sinθ1+ z cosθ 1 = R c o s θt cosφ cosφb sinθ1 + R cosθ t sinφ sinφb sinθ1 + R sinθt cosθ1 = R [ c o s θt sinθ1 (cosφ cosφb + sinφ sinφb ) + sinθt cosθ1 = R [ c o s θt sinθ1 cos(φ -φb) + sinθt cosθ1] = R [ s i n θt cosθ1 + cosθt sinθ1 cos(φ-φb)] The last three results may be summarized ( r' ≈ x' e' 1 + y' e'2 + z' e'3 ) , x' = R cos θ t sin(φ –φb) y' = R [sin θt sinθ1 – cosθt cosθ1 cos(φ –φb) ] z' = R [sin θt cosθ1 + cosθt sinθ1 cos(φ-φb) ] ( 7 . 1 4 ) where θ t and φ are given by (6.13) and φb by (7.3). Equations (7.14) give the Cartesian coordinates of the sun's location on the celestial sphere in Frame S' which rotates with the earth. As shown in Fig 7.1, the phase (φ b–φ) takes all values in the range 0,2 π during a 24 hour day, with φb doing most of the moving (in the positive direction). If we are only interested in using (7.14) to plot the daily trajectory of the sun, we can treat p ≡ (φ b–φ) as a generic parametric argument p which runs 0 to 2 π Section 7: Coordinates of the Sun viewed from Frame S' 66 during the day. Such an argument can then be written p = ωt as t runs 0 to 24 hours( ω = 2π/24 rad/hr). We then obtain this version of (7.14) suitable for plotting the daily trajectory of the sun : x' = – R cos θt sin(ωt) // valid only for making trajectory plots y' = R [sinθ t sinθ1 – cosθt cosθ1 cos(ω t) ] z' = R [sinθ t cosθ1 + cosθt sinθ1 cos(ωt ) ] . ( 7 . 1 5 ) One might wonder now about the location of the mean fictitious sun in Frame S'. Recall from (7.3) the right ascension of the observer site (projected onto the celestial sphere) φ b(τ) = φb0 + ωsτ . (7.3) For the mean fictitious sun we have from (6.14) that φ f = Ω0t - α (6.14) but we want this in terms of time τ, so using (7.2) and t Me as in (4.12), φ f = Ω0(τ + tMe) - α = Ω0τ + [Ω0 tMe - α] = Ω0τ + k k ≡ [Ω0 tMe - α] = -.032437114 rad (2012) . (7.16) Angle k is thus the right ascension of the mean fictitious sun at March equinox τ = 0. Then the fictitious sun is described by φf(τ) = Ω0τ + k // radians k = -.032437114 rad (2012) θt,f(τ) = 0 // θ tmax = 0 ( 7 . 1 7 ) so that, using the fact (1.6) that ω s - Ω0 = ω, for the mean fictitious sun we get (φ f-φb) = Ω0τ + k – [φb0 + ωsτ] = –ωτ + (k - φb0) ( 7 . 1 8 ) and then (7.14) for the fictitious sun position in Frame S' reads, x'f = – R sin( ωτ – k + φb0) y'f = R [–cos θ1 cos(ωτ – k + φb0) ] z'f = R [ sin θ1 cos(ωτ – k + φb0)] . // mean fictitious sun (7.19) Viewing the fictitious sun from the north pole one gets x' f = – R sin ( ωτ – k + φb0) y'f = – R cos( ωτ– k + φb0) ] z'f = 0 // mean fictitious sun, North Pole (7.20) Section 7: Coordinates of the Sun viewed from Frame S' 67 As τ increases, these equations describe the fictitious sun rotating clockwise around the north pole, so as usual the sun rises in the east, Fig 7.2 The velocity and acceleration of both the actual an d fictitious sun's apparent motions on the celestial sphere (as seen by our Frame S' observer) can be obtained by differentiating (7.14) and (7.19). The velocity requires expressions for θ• t and φ• (φ•b = ωs) obtained from (6.13), and then acceleration requires θ•• t and φ••. The results are easily obtained in Maple and then approximations can be made as needed. (d) Spherical coordinates θ',φ' of the sun in Frame S' Warning: The angle φ' here is totally unrelated to the angle called φ' which appears in Figs 1.3, 3.4 and 6.1 and which is used throughout Section 6. We can convert the above equations to a spherical coordinate system (r', θ',φ') which has its origin at the Frame S' observation point on the surface of the earth with the z' axis pointing up: Fig 7.3 In this system, azimuth φ' = 0 is the east, and increases positively to the north. The noon sun always has φ' = ± π /2. At equinox, the sun rises at φ' = 0 and sets at φ' = π. The general plan for converting to spherical coordinates starts as follows : Section 7: Coordinates of the Sun viewed from Frame S' 68 ρ'2 = x'2+ y'2 r'2 = x'2+ y'2 + z'2 . ( 7 . 2 1 ) Since r' lies on a celestial sphere of radius R, we know ah ead of time that r' = R and can confirm that in Maple using the x',y',z' expressions gi ven in (7.14). On the other hand, ρ' is a bit messy, (ρ'/R)2 = (x'/R)2 + (y'/R)2 = [cosθt sin(φ –φb)]2 + [sinθt sinθ1 – cosθt cosθ1 cos(φ –φb)]2 = c o s 2θtsin2(φ –φb) + cos2γ ( 7 . 2 2 ) where cosγ ≡ sinθ t sinθ1 – cosθt cosθ1 cos(φ –φb) . ( 7 . 2 3 ) The relation between Cartesian and spherical coordinates is x' = r'sin θ'cosφ' y' = r'sin θ'sinφ' z' = r'cos θ' where sin θ' is always positive since θ ' lies in (0, π). Meanwhile, ρ' 2 = x'2+ y'2 = r'2sin2θ' => ρ' = r' sinθ' , ( 7 . 2 4 ) so that x' = ρ'cosφ' y' = ρ'sinφ' z' = r'cos θ' (7.25) and we end up with cosθ' = z'/r' cos φ' = x'/ρ' tan φ' = y'/x' sinθ' = ρ'/r' sin φ' = y'/ρ' r ' = R ( 7 . 2 6 ) Installing the expressions (7.14) we then find that cosθ' = [sinθ t cosθ1 + cosθt sinθ1 cos(φ-φb)] sinθ' = ρ'/r' = cos2θtsin2(φ –φb) + cos2γ ( 7 . 2 7 ) cosφ' = x'/ρ' = [cosθt sin(φ –φb)] / cos2θtsin2(φ –φb) + cos2γ sinφ' = y'/ρ' = [sinθt sinθ1 – cosθt cosθ1 cos(φ –φb) ]/ cos2θtsin2(φ –φb) + cos2γ tanφ' = y'/x' = [sin θt sinθ1 – cosθt cosθ1 cos(φ –φb) ]/ [cosθt sin(φ –φb)] (7.28) Section 7: Coordinates of the Sun viewed from Frame S' 69 Therefore, we may write θ' = cos-1[sinθt cosθ1 + cosθt sinθ1 cos(φ-φb)] φ' = tan-1 ( [sinθt sinθ1 – cosθt cosθ1 cos(φ –φb) ] / [cosθt sin(φ –φb)] ) . (7.29) In the above equations, the following objects appear: θt = sin-1[sinθtmax sin(ψ - α) ] (6.13) φ = tan-1[cosθtmax tan(ψ - α) ] . (6.13) ψ = Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) // + order e3 which we ignore (4.31) α = π/2 - φper (1.9)) t = τ + tMe (7.2) φb = φb0 + ωsτ = φ1 + φG0 + ωsτ . (7.3),(7.7) As the sun moves in the sky, equations (7.29) gi ve the position of the sun in Frame S' spherical coordinates based at an observation point located at ( θ1,φ1), where θ1 is the polar angle from the north pole (colatitude) and φ1 is the longitude east of the Greenwich meridian. Angle φG0 is the right ascension of Greenwich at the last March equinox. Section 8: The Equation of Time Revisited 70 8. The Equation of Time Revisited (a) Definition of new azimuths Φ ' and Φ 's Recall the depiction of the Frame S' coordinate system given in Fig 7.3 which we replicate here with the addition of two new azimuthal angles Φ' and Φ's, Fig 8.1 For an observer in the southern hemisphere below latitude - 23o, the azimuth φ' of the sun follows a regular pattern all year long: the sun rises generally in the east at φ' = 0 ± Δ , reaches noon in the north at φ' = π/2, sets generally in the west at φ' = π ∓ Δ, and hits midnight at φ' = 3π/2 = -π/2. The size of Δ varies with latitude and the time of year and is ze ro at equinoxes and maximal at solstices, as shown in Fig 9.5 below. For such an observer we might define a shifted azimuth in order to place midnight at the discontinuity point where 2 π → 0, Φ's ≡ π/2 + φ' . ( 8 . 1 ) The sun rises in the east at Φ' s = π/2 ± Δ, reaches noon in the north at Φ's = π, sets in the west at Φ 's = 3π/2 ∓ Δ, and hits midnight in the north at Φ's = 2π = 0. For an observer in the northern hemisphere above latitude 23 o, the appropriate sun-tracking azimuth would be this, Φ' ≡ π/2 – φ ' . // Φ's = π – Φ' ( 8 . 2 ) The sun rises in the east at Φ' = π/2 ∓ Δ , reaches noon in the south at Φ' = π, sets in the west at Φ' = 3π/2 ± Δ, and hits midnight in the north at Φ' = 2π = 0. For viewers in the tropic region between -23o and +23o, the "sun tracking azimuth" is Φ's for part of the year, and Φ' for another part of the year. For example, the Honolulu image in Fig 3.12 shows that Φ' is Section 8: The Equation of Time Revisited 71 appropriate for most of the year when the sun is in the south, but for part of the year near June solstice the Φ's azimuth is appropriate. polar = 70 lat = 20 Honolulu, Mexico City, San Juan Fig 8.2 As we shall see in Section 11, it is not very useful to study the analemma for a tropical viewer using Frame S' polar coordinates θ',φ' using any of the three azimuths, but the direct calculation of the analemma "on film" using Cartesian c oordinates x',y',z' always works. We continue now assuming a north ern hemisphere observer above 23 o latitude, mainly because we shall be interested in the limiting case of the Nort h Pole, and because we shall be studying analemmas photographed by such observers. Using equation (7.29), φ' = tan -1 [sinθt sinθ1 – cosθt cosθ1 cos(φ –φb) ]/ [cosθt sin(φ –φb)] , (7.29) we may relate Φ' to the position of the actual sun on the ecliptic this way, using (8.2), Φ' = π/2 – tan -1 ( [sinθt sinθ1 – cosθt cosθ1 cos(φb-φ) ]/ [–cosθt sin(φb-φ)] ) . Then using the identity (A.12), π/2 – tan -1([y]/[x]) = π + tan-1( [ - x ] / [ - y ] ) (A.12) this becomes, for the actual and fictitious suns, Φ' = π + tan -1 ([cosθt sin(φb-φ)] / [– sinθt sinθ1 + cosθt cosθ1 cos(φb-φ) ] ) Φ'f = π + tan-1 ( [sin(φb-φf)] / [cosθ1 cos(φb-φf) ] ) . (8.3) Using the fact that Φ' s = π – Φ' , and the second line of (A.4) inverted, tan -1( [y]/[x]) = π - tan([y]/[-x]), (A.4) and the fact that - π and π represent the same angle, we can write for southern observers below -23 o, Section 8: The Equation of Time Revisited 72 Φ's = π + tan-1 ([cosθt sin(φb-φ)] / [sinθt sinθ1 – cosθt cosθ1 cos(φb-φ) ] ) Φ's,f = π + tan-1 ( [sin(φb-φf)] / [– cos θ1 cos(φb-φf) ] ) (8.4) where the difference from (8.3) is that both bottom brackets are negated. At the north pole (8.3) becomes (on the second line, tan-1( tan(φb-φf) ) = (φb-φf) ) Φ'NP = π + tan-1 ([cosθt sin(φb-φ)] / [cosθt cos(φb-φ) ] ) = π + (φb-φ) Φ'fNP = π + tan-1 ( [sin(φb-φf)] / [cos(φb-φf) ] ) = π + (φb-φf) (8.5) Taking the difference, we find that Φ' NP – Φf'NP = – (φ – φf) o r ΔΦ'NP = – Δφ . => Φ'NP = Φf'NP – Δφ (8.6) where Δφ is the "equation of time" shown in (6.15a), Δφ(t) ≡ φ - φ f = tan-1[cosθtmax tan(ψ-α)] – (Ω0t-α) (6.15a) ψ = Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) // + order e3 which we ignore . (4.31) Recalling now that φ f(τ) = Ω0τ + k k = -.032437114 rad (2012) (7.17) φb(τ) = φb0 + ωsτ . (7.3) We can write the second line of (8.5) as Φ' fNP = π + (φb-φf) = π + [φb0 + ωsτ] - [Ω0τ + k] = π + (ωs- Ω0)τ + (φb0-k) = π + ωτ + (φ b0-k) // using (1.6) (8.7) (b) Local Time, Sundial Time and the Equation of Time GMT. For our purposes here, we make no distinction be tween the various kinds of precision global time indicated by the letters GMT, UT, UT0, UT1, UTC and so on, so we will just use GMT. We are not concerned about tiny differences between these times due to things like the Chandler wobble mentioned at the end of Section 1 (a). GMT means Greenwich Mean Time which is a 24 hour time system set so that the mean fictitious sun always passes over the Green wich meridian (the prime meridian, longitude = 0 0) at GMT = 12:00. At this instant, the entire world experiences the same "day", perhaps Thursday. At GMT = 12:01 the world is still doing Thursday, except a tiny strip on the meridian opposite Greenwich has moved to Friday, and then Friday "grows" as time moves on. This meridian at φ1 = 180o is the GMT international date line, though Section 8: The Equation of Time Revisited 73 the actual date line has a deviation in the southern hemi sphere. Thus it is that Samoa rings in the new year Jan 1 while it is still GMT noon Dec 31 in London (and 4AM in California). LMT. The corresponding time at an arbitrary Frame S' observer site is LMT meaning Local Mean Time, and it is given by LMT = [ GMT + ( φ1/15) ] mod 24 . (8.8) where φ 1 is the site's longitude east of Greenwich, meas ured in degrees. At our observer site, the mean fictitious sun always crosses the local meridian at LMT = 12:00. So GMT is just LMT for Greenwich. This LMT should not be confused with time zone times which are binned every 15o (more or less). Our observer can set a local clock to his LMT using the above equation and obtaining live GMT from a website displaying GMT, from a GMT "clock app", or from a WWV broadcast (where GMT is called "coordinated universal time" or UTC). We know from (8.7) that Φ' fNP advances at the linear rate ωτ , so it can be related to LMT which advances at that same linear rate, albeit scaled to hours instead of degrees. We claim in fact that LMT = 12(Φ ' fNP/π) / / L o c a l M e a n T i m e ( 8 . 9 ) We know that Φ' fNP = π when the fictitious sun crosses the noon meridian, and there LMT = 12 = 12:00. Midnight occurs when Φ'fNP = 0 which is then LMT = 0. Since Φ 'f is linear, (8.9) is thus proved. Let us now define another time LST based on the azimuth of the actual sun, LST ≡ 12(Φ'NP/π) // Local Sundial Time (8.10) This LST clock reads 12:00 when the actual sun crosses the local meridian, Φ ' NP = π. This is the time when a local sundial vane shows no shadow to eith er side, hence the name "local sundial time". Building an electronic LST clock: One could imagine constructing an electronic clock to display LST by having a feedback servo system which very slowly slews the clock rate so that it stays in daily sync with the event of the actual sun crossing the meridian every noon time (LST = 12:00). The design might be sophisticated in that it can remain in "lock" through some number of completely overcast days due to prediction based on recent past history. The clock rate is uneven because, during the year, the actual sun slowly wanders back and forth in right ascension relati ve to the mean fictitious sun by about ±15 minutes as indicated by the Equation of Time shown as the red curve in Fig 6.6. Our electronic LST clock is "better" than an actual sundial only because it displays time at all times, day and night, and even works on overcast days. Using (8.9), (8.10) and (8.6), we find that LST – LMT = 12( ΔΦ' NP/π) = – 12(Δφ /π) . // hours and radians (8.11) Section 8: The Equation of Time Revisited 74 Here are four ways to write this last equation [LST – LMT]hrs = – 12(Δφrad/π) ( 8 . 1 2 ) [LST – LMT] hrs = – 12(Δφdeg/180) (8.13) [LST – LMT] min = – 60*12(Δφrad/π) = – (720/π )Δφ ( 8 . 1 4 ) [LST – LMT] min = –60*12(Δφdeg/180) = – 4 Δφdeg = – Δφmin . (8.15) It is most common to express the equation of time in minutes, and from (8.14) and (6.15a) we write EOT(t) ≡ [LST – LMT]min = – (720/ π)Δφrad = - 4 (180/π ) Δφ(6.15a) = – (720/π ) { tan-1[cosθtmax tan(ψ-α)] – (Ω0t-α) } . (8.16) A plot of this EOT(t) ~ – Δφ = φf-φ was already shown in Fig 6.7. Here we replicate that plot adding some horizontal gridlines F i g 8 . 3 The EOT(t) plot for the actual sun is of course the re d curve, while the green line represents the fictitious sun. Recall that the left edge corresponds to t = 0, Jan 3, perihelion. On the upper half of this plot, we have – Δφ > 0 so LST > LMT which means the sundial tim e is ahead of LMT, the sundial is "running fast", and the actual sun is ahead of the fictitious sun on the celestial sphere, φf > φ . The opposite of course is true for the lower half of the plot. Section 9: Noon and Midnight, Sunrise and Sunset 75 9. Noon and Midnight, Sunrise and Sunset Angle θ1 continues here to be the colatitude of the Frame S' observation site. Letting θLAT be the corresponding latitude , so that θLAT = π/2 - θ1, one can write all the equations below with these simple changes, sinθ1 → cosθLAT cosθ1 → sinθLAT cotθ 1 → cotθLAT . (a) Conditions for Noon and Midnight The path of t he sun in local Frame S' coordinates is given by θ' = cos-1[sinθt cosθ1 + cosθt sinθ1 cos(φ-φb)] φ' = tan-1 [sinθt sinθ1 – cosθt cosθ1 cos(φ –φb) ]/ [cosθt sin(φ –φb)] (7.29) where θ t = sin-1[sinθtmax sin(Ω0t + 2e sin( Ω0t) - π/2 + φper)] φ = tan-1[cosθtmax tan(Ω0t + 2e sin( Ω0t) - π/2 + φper) ] (6.13) φb = φb0 + ωsτ = φ1 + φG0 + ωsτ with τ = t - tMe (7.3), (7.7), (7.2) and here is the picture showing the obs erver site Frame S' coordinates, Fig 7.3 Noon occurs when the local polar angle θ' assumes its smallest value during the day; the sun is as much "overhead" as it will be for that day; it is at its closest to the zenith z' axis. Looking at the first of equations (7.29) quoted above, cosθ' = [sinθ t cosθ1 + cosθt sinθ1 cos(φ-φb)] , (9.1) Section 9: Noon and Midnight, Sunrise and Sunset 76 we naively ask: what is the value of ( φ-φb) at noon? The quantities cos θt and sinθ 1 are always positive, so cosθ' is maximum when cos( φ-φb) = 1. Pondering the shape of cos θ' for θ' in 0,π (its legal range), one concludes that θ' is minimal when cos θ' is maximal. Therefore: noon: cos(φ -φb) = 1 φ = φb φ = (6.13) quoted above midnight: cos( φ-φb) = -1 φ = φb + π (9.2) These both imply that sin( φ-φ b) = 0. Looking at the second of equations (7.29), tan φ' = [sinθ t sinθ1 – cosθt cosθ1 cos(φ –φb) ]/ [cosθt sin(φ –φb)] (7.29) we conclude that at noon, tanφ ' = ±∞ so φ' = ±π/2. According to the figure above, this means that at noon the sun lies on the local meridian , a north-south line of longitude passing through the observation site. Having "done the math", we now peruse Fig 7.1 where things are completely obvious : Fig 7.1 When φb = φ, the Frame S' origin on the earth and the sun on the celestial sphere have the same azimuth angle, so this is noon. Midnight occurs when φb = φ + π and therefore cos(φ -φb) = -1. The conditions for noon and midnight of the mean fictitious sun are given by fictitious noon: cos(φ f-φb) = 1 φf = φb fictitious midnight: cos(φ f-φb) = -1 φf = φb + π (9.3) where Section 9: Noon and Midnight, Sunrise and Sunset 77 φf = Ω0t - α = Ω0t - π/2 + φper = Ωoτ + k with k = [ Ω 0tMe - α ] (6.14),(7.16) φb = φb0 + ωsτ = φ1 + φG0 + ωsτ with τ = t - tMe (7.3), (7.7), (7.2) (b) Local time of Noon and Midnight : approximations From (8.12) we know that local mean time and local sundial time are related by LMT = [LST + (12/π )Δφ] mod 24 , (9.4) where Δφ is the Equation of Time (in the form (6.15a)), Δφ = -φf + φ = – Ω0t +α + tan-1[cosθtmax tan(ψ-α) ] . (6.15a) Noon means the sundial has no side shadow and so LST = 12:00. From (9.4) this occurs at local time, LMT(noon) = [12:00 + (12/ π)Δφ] mod 24 (9.5) and we know this LMT will lie in the general range 11:45 to 12:15. Similarly, LMT(midnight) = [00:00 + (12/ π)Δφ)] mod 24 . (9.6) These results can be converted to GMT using (8.8) written as GMT = [ LMT – ( φ 1/15)] mod 24 (9.7) where φ 1 is the observer's longitude east of Greenwich in degrees. For 2012 parameters, here is a plot of the LMT(noon) from (9.5) and (6.15a) above Section 9: Noon and Midnight, Sunrise and Sunset 78 F i g 9 . 1 The vertical lines were placed using days computed in (4.11). Note that 0.1 hour = 6 minutes. This is just a restatement of the equation of time Fi g 6.6 with some LMT clock times added on the left. On September equinox, at solar noon the LMT clock r eads 11:52 AM so "the sundial is running fast". Precision calculation of solar noon. In the discussion just above, on e can roughly look up on the graph Δφ for a day of interest and use that to compute LMT( noon) using (9.5). To get a more precise result, one would use equation (6.15a) above to compute Δφ, but the question arises: exactly what t should one use in (6.15a)? If one knew the exact t of solar noon on our day of interest, that would be the right t to use, but of course one doesn't know the exact t of solar n oon! We then have to make an approximation. In Δφ(t=τ+tMe) we might set τ = 24d hours for day "d" of the year starting at March equinox. This approach corresponds to the left drawing below, Fig 9.2 where the red line shows how we approximate Δφ as being a constant for the entire day d. In this approach, things are a little hazy at the left end of the graph near τ = 0, but since Δφ is a slowly varying function over a range of a few days, this approximation will work pretty well. A better approximation (right drawing) is to set τ in Δφ(t=τ+t Me) to be the time of noon for the mean fictitious sun, which we know is going to be ± 15 min from the noon of the actual sun. This occurs according to (9.3) when φf = φb + 2πd which from (7.1) tell us – ωτ + (k - φ b0) = 2πd. This then gives Section 9: Noon and Midnight, Sunrise and Sunset 79 the time of fictitious sun noon to be τ = 12 (k - φb0)/π + 24d hours, where φb0 = φ1 + φG0 as in (7.7) and constant k is given in (7.16) as -.03243 radians for year 2012. Quantity φ G0 is the right ascension of Greenwich at τ = 0. Setting d = 0, we find that the time of fictitious sun noon on day 0 occurs at τ0 = 12(k - φb0)/π. For the picture as drawn, this would be a positive number. If the τ = 0 line falls in the right half of day d = 0, then this τ0 is a negative number. Finally, one could find a numerically exact solution for t of actual solar noon on day d by solving a transcendental equation for a value of t on each day. That equation is just (9.2), φ = φb (9.2) Replacing φ from (6.13) and φ b from (7.3), and with t Me from (4.12), the above becomes tan -1[cosθtmax tan{ Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) -π/2+φper}] = φb0 + ωst – ωstMe or cosθ tmax tan{ Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t) -π/2+φper} = tan(φb0 + ωst – ωstMe) where one sees t appearing on both sides ( t = τ + t Me from (7.2)). (c) Another approach to noon As noted above, noon occurs when ( φ-φb) = 0. We can write an expression for ( φ-φb) as follows, using the facts that φ = φf + Δφ, φb = φb0 + ωsτ (7.3) and φf = Ω0τ + k (7.17) : (φ-φ b) = φf + Δφ - ωsτ - φb0 = Ω0τ + k + Δφ - ωsτ - φb0 = - ωτ +k+ Δφ - φ b0 // since ωs = ω + Ω0 as in (1.6) = - ωτ - [φ b0 - k – Δφ]] = - ωτ- κω // κ ≡ [φ b0 - k – Δφ]/ ω = - ω(τ+κ) = - ωt ' ( 9 . 8 ) where we define a shifted time t' according to t' = τ + κ κ = [φ b0 - k – Δφ ]/ω . ( 9 . 9 ) In this shifted time t', which is specific to our selected day (since Δφ is specific), noon occurs at t ' = 0 / / n o o n ( 9 . 1 0 ) since this makes φ = φ b. This fact is not very exciting, but becomes useful in the next section. Section 9: Noon and Midnight, Sunrise and Sunset 80 (d) Sunrise and Sunset, Length of Day, Noon Elevation of the Sun Sunrise and S unset Looking at Fig 7.3 above, the condition for nominal sunrise and sunset on a perfectly smooth spherical earth is θ' = π/2 . => cos θ' = 0 . ( 9 . 1 1 ) Nominal means that the center of the sun passes through the horizon plane θ' = π/2. Looking then at the first equation of (7.29) we obtain th e condition (sr = sunrise, ss = sunset) [sinθ t cosθ1 + cosθt sinθ1 cos(φ-φb)sr,ss ] = 0 ( 9 . 1 2 ) so that cos(φ-φb)sr,ss = – sinθ t cosθ1 / (cosθt sinθ1) = – tanθ t cotθ1 . (9.13) Since φ and θ t are both complicated functions of time t, as shown in (6.13), this is a transcendental equation one would have to solve numerically for the times of sunsets and sunrises. We can, however, use either of the approximation methods discussed in section (b) above. The first method would just set τ = 24d to evaluate θt(t=τ+tMe). As for the second method, using the fact that (9.13) says cos(φf-φb)sr,ss = 0 for the fictitious sun, and the fact that ( φf-φb) = –ωτ + (k - φb0) from (7.18), one finds that ( φf-φb)sr,ss = ± π/2 - 2πd where d = days, and then τ± = 12(k - φ b0)/π ∓ 6 + 24d where the - is for sunrise and the + for sunset. This then just says τ± = τnoon ∓ 6 hours, which seems pretty reasonable (all for the fictitious sun). So, for our selected day, we compute θ t(t=τ+tMe) from (6.13) and then sunrise and sunset occur when cos(φ-φ b)sr,ss = – tanθt,sr,sr cotθ1 (9.13) or, using (9.8), cos(ωt' sr,ss ) = – tanθ t,sr,ss cotθ1 . ( 9 . 1 4 ) If we assume the simpler approximation τ = 24d, then θt is the same for sunrise and sunset so that t' sr = – (1/ω ) cos-1[– tanθt cotθ1] < 0 // sunrise θ1 = colatitude t'ss = + (1/ω) cos-1[– tanθt cotθ1] > 0 // sunset o r ( 9 . 1 5 ) t' sr = – (1/ω ) cos-1[– tanθt tanθLAT] < 0 // sunrise θLAT = latitude t'ss = + (1/ω) cos-1[– tanθt tanθLAT] > 0 // sunset Section 9: Noon and Midnight, Sunrise and Sunset 81 where ω = 2π rad/day = ( π/12) rad/hr. These times are relative to t' = 0 as noon. We know the LMT of noon from (9.5) and then (9.15) tells us the displaceme nt from that noon time of sunrise and sunset. In an exact numerical solution (or using the more accurate approximation discussed above), one would find that these displacements are not quite the same because θt changes slightly between sunrise and sunset, but we ignore this fact in our simple approximation as illustrated in Fig 9.2. Looking at (9.14), one sees there is no solu tion t' for sunrise or sunset when | tan θt cotθ1 | > 1. This occurs in the northern hemisphere for observers with θ1 < θt . For them the sun is either up all the time or down all the time so there are no sunrises or sunsets. The corresponding condition in the southern hemisphere is θ1 > π-θt . Here is a plot of t' ss for a selection of northern hemisphere latitudes which uses (9.15) for t' ss and (6.24) for θt [tp_ss = t' ss, phip_sr = φ'sr for making Fig 9.4 below ] Fig 9.3 Section 9: Noon and Midnight, Sunrise and Sunset 82 For the southern hemisphere, negate latitude labels in Fig 9.3 and refl ect all curves vertically through the horizontal line at t ss = 6 hours. This follows from (9.15) and the rule cos-1(-x) = π - cos-1(x). Example: In London ( θLAT = 51.5o) at June solstice, t' ss ≈ 8.2 = 8:12. From Fig 9.1, noon occurs at LMT = 12:02, so nominal sunset would then be at about 8:14 PM. Since London is on daylight savings time, clock time of nominal sunset would then be 9:14PM. The official value is 9:21 PM, about 7 minutes later. Sunrise time would be 12:02 - 8:12 = 3:50 AM, but DST moves this to 4:50 AM. The official number is 4:43 AM, about 7 minutes earlier (see comments in next section). http://www.sunrisesunset.com/England/ Length of Day The nominal length of day is given by, t day = 2 t'ss = (24/π ) cos-1[- tanθt cotθ1] . ( 9 . 1 6 ) Here we have assumed that "day" is twice the noon-to-s unset time, where sunset means the point at which the center of the sun passes through the horizon. Various adjustments are usually made to this result: (1) sunrise is usually defined as occurring when the upper limb of the sun first appears, and sunset when the upper limb disappears, since this controls "daylight". This correction is a function of the angle the sun makes with vertical at the horizon when it rises and sets (see below); (2) when the sun is just below the horizon, it appears above the horizon due to refraction in the atmosphere. The diamet er of the sun is about a half degree ( ~ 30 arc minutes) and the refraction eff ect is about the same size . Both these effects make the official "daylight" longer than our nominal value t day. As a ballpark estimate, these effects might add 2 degrees to daylight or 8 minutes to the nominal value given above. Example: At summer solstice, an observer at θLAT = 40o north latitude would have, t day = (24/π )cos-1[- tan(+23.4382o) tan(40o)] = 14.844 hours = 14 hr 51 min tdaylight ≈ tday + 8 min ≈ 15 hrs This would be the longest daylight of the year for this observation site. Sunrise azimuth We can install cos( φ-φb)sr of (9.13) into (7.29) to find the sun' s azimuth at sunrise. From (9.8) we know that (φ-φb)sr = -ωt'sr > 0 so sin( φ –φb)sr > 0. Then, Section 9: Noon and Midnight, Sunrise and Sunset 83 tanφ'sr = [sinθt sinθ1 – cosθt cosθ1 cos(φ –φb)sr ] / [cosθt sin(φ –φb)sr] = [ s i n θt sinθ1 + cosθt cosθ1 tanθt cotθ1 ] / [cosθt 1 - tan2θt cot2θ1 ] = [ s i n θ t sinθ1 + cosθt cosθ1 (sinθt/cosθt) (cosθ1/sinθ1)] / [cosθt 1 - tan2θt cot2θ1 ] = [ s i n θ t sinθ1 + sinθt cos2θ1/sinθ1] / [cosθt 1 - tan2θt cot2θ1 ] = t a n θt[sinθ1 + cos2θ1/sinθ1] / 1 - tan2θt cot2θ1 = (tan θt/sinθ1) / 1 - tan2θt cot2θ1 = (tan θ t secθLAT) / 1 - tan2θt tan2θLAT θLAT = latitude = π/2-θ1 (9.17) Notice that φ 'sr is an odd function of θt. This fact appears in the red circle pictures of Fig 3.12 where the azimuth of sunrise at winter and summer solstice are equal and opposite. For example, with θLAT = 40o, F i g 9 . 4 In this particular case the φ' sunrise azimuth at summer solstice is +31.3o north of east, and therefore the sunrise azimuth at winter solstice is - 31.3 o. Here is plot of sunrise azimuth as a function of tim e of year for a set of northern hemisphere latitudes, using equation (9.17) for φ'sr and (6.13) for θ t. The code is as shown above Fig 9.3 with the addition of these lines [ phip_sr = φ'sr ], Section 9: Noon and Midnight, Sunrise and Sunset 84 F i g 9 . 5 The curve for θ LAT = 400 latitude peaks at ~31o at June solstice, in agreement with the calculation just done above. For latitudes in the arctic, the curves end where there is no daily sunrise. In the plotting mechanism, the curves end because the function goes complex. The graphs are exactly the same for southern hemisp here latitudes, just negate the latitude labels. This follows from (9.17) which is even under θ LAT → - θLAT. Sunset azimuths are the same numbers, but are measured north of west instead of north of east. Section 9: Noon and Midnight, Sunrise and Sunset 85 Noon Elevation of the Sun The elevation of the sun above the horizon in terms of the sun's polar angle is given by θ el = π/2 - θ ' . ( 9 . 1 8 ) Then using (9.1) with (9.2), φ = φ b, one finds sin(θ el noon) = cosθ'noon = [sinθt cosθ1 + cosθt sinθ1] = sin(θt+θ1). (9.19) The correct solution of sin( θ el noon) = sin(θt+θ1) for θel noon in this situation is as follows: θ 1+θt range: 0 π/2 π θel noon = : no noon θ1+θt π - (θ1+θt) no noon (9.20) For example, if θ 1+θt < 0 or θ1+θt > π there is no noon because the sun never rises. If θ1+θt lies between 0 and π/2, then θel noon = θ1+θt. Examples: θ1 = 50o (40o north latitude) θ1+ θt θel noon June solstice ( θt= 23.4o) 73.4o 73.4o high in sky equinoxes ( θt= 0) 50o 50o Dec solstice ( θt= -23.4o) 26.6o 26.6o low in sky θ 1 = 130o (40o south latitude) θ1+ θt θel noon June solstice ( θt= 23.4o) 153.4o 26.6o low in sky equinoxes ( θt= 0) 130o 50o Dec solstice ( θt= -23.4o) 106.6o 73.4o high in sky Section 10: Insolation 86 10. Insolation (a) The t ss function and hours of daylight In (9.15) the time t ss from noon to sunset is given by ( we drop the prime on t ss) t ss = (12/π ) cos-1[- tanθt cotθ1] ( 1 0 . 1 ) where cos -1 is in radians and t ss in hours. Here θt is the declination (on the celestial sphere) of the sun at some time of year as shown in Fig 6.1, and θ1 is the colatitude (polar angle) of the observer site. Latitude θLAT = π/2 - θ1 is often used below: sin θ1 = cosθLAT, cosθ1 = sinθLAT, cotθ1 = tanθLAT. Comment : Although cos(-x) = cos(x), it happens that cos-1(-x) = π - cos-1(x), so one cannot "pull out" the minus sign in (10.1). On the other hand, both tan(x) and tan-1(x) are odd functions of x. This function t ss(θt,θ1) plays an important role in the discussion below, so we want to write a precise expression for it, including the regions in which it is pi nned to 12 or 0 hours, meaning the sun is up all the time or down all the time. Here is that e xpression (see Section 3 (f) and equation (3.2)), north pole pin south pole pin t ss = 12 if: ( θt > 0 AND θ1 < θt) OR (θt < 0 AND θ1 > π - |θt| ) south pole pin north pole pin t ss = 0 if: ( θt > 0 AND θ1 > π - θt) OR (θt < 0 AND θ1 < |θt| ) t ss = ω-1 cos-1[- tanθt cotθ1] = ω-1 cos-1[- tanθt tanθLAT] otherwise . (10.2) Eq. (10.2) is implemented as t ss(θt,lat) (θt radians, lat degrees) by this chunk of Maple code, where th1 = θ1, tht = θt, thLAT = θLAT, and lat = latitude : A 3D plot of this t ss function is rather interesting, Section 10: Insolation 87 Fig 10.1 This plot shows the regions where t ss is pinned to 12 or 0. When t ss = 12, the sun never sets and "day" is a full 24 hours long. When t ss = 0 the sun never rises. By doubling the vertical axis, the above plot shows nominal daylight as a function of θt= tilt angle and θ LAT = latitude. Recall equation (5.5) for θt, where we now ignore the eccentricity of the earth's orbit, θ t(t) ≈ sin-1 { sinθtmax sin[Ω0t + φper - π/2] } . // assumes e = 0 (5.5) In this equation, t = 0 corresponds to perihelion and Ω 0 = 2π rad/year. The argument may be written [Ω 0t + φper - π/2] = Ω0(τ + tMe) - α // (7.2) and (1.9) = Ωτ + (Ω 0 tMe - α) = Ωτ // (4.12) and our equation above can be expressed in terms of τ as [θ t(t) and θt(τ) are different functions ] θ t(τ) ≈ sin-1 { sinθtmax sin(Ω0τ) } // e = 0, τ = 0 is Me (10.3) where τ = 0 is March equinox. Note that θ t = 0 at τ = 0 as one would expect from Fig 3.5. For insolation purposes, we want τ = 0 to be the June solstice, the time when θt is maximal. With this new time origin, the above becomes θt(τ) ≈ sin-1 { sinθtmax cos[Ω0τ] } // e = 0, τ = 0 is Js (10.4) or θ t(τ) ≈ sin-1 { sinθtmax cos[2π (d/365.25)] } (10.5) Section 10: Insolation 88 where d is the days elapsed since the June solstice. Using (10.4) for θt in (10.2), we obtain a 3D plot of tss where the θt axis is replaced by a days d axis, Fig 10.2 Below is a slice of 3D Fig 10.2 taken in mid plot parallel to the left face (June solstice) Fig 10.3 The curve is of course continuous but has serious kinks at the pinning points. The plot below gives cross sections of Fig 10.2 parallel to the back face (but doubled in height). Section 10: Insolation 89 F i g 1 0 . 4 The rectangular curve is for latitude = 90 o (north pole, back wall of Fig 10.2). This curve is rectangular because at the north pole, the (nominal) sun is either up all day or down all day. The horizontal line for latitude = 0o reminds us that the day is always 12 hours l ong at the equator. These facts can be visualized by putting an observer on the north pole or equator in Fig 3.13 or 3.14 and remembering that the earth should be shrunk down to a small dot. Only the curves of latitude above the Arctic Circle pin at 24 hours on the top and pin at 0 hours on the bottom. In the southern hemisphere the above plot is th e same but one should nega te all latitude labels. For some external verification of Fig 10.4, here is a web plot showing latitudes 70,50,34,23 and 0 : Section 10: Insolation 90 Fig 10.5 http://www.geog.ucsb.edu/ideas/Insolation.html In the above we have assumed that "day" is twice the sunrise time, where sunrise means the point at which the center of the sun passes through the horizon (same for sunset). The ~ 8 minute adder correction for solar diameter and refraction was discussed in Section 9 (d). These effects and the small time distortion effect of orbital eccentricity are ignored in this entire Section. (b) Daily Insolation versus latitude and day of year The Fra me S' polar angle θ' of the sun is given by (7.29) where we set ( φ-φb) = -ωt' as in (9.8), where t' is a shifted time which causes noon to occur at t' = 0, cosθ' = [sinθt cosθ1 + cosθt sinθ1 cos(ωt')] . (10.6) The solar power falling on 1 square meter of earth surface is proportional to cos θ'. This is simply because the area presented by 1 square meter of ground to the sun at angle θ' from vertical is only cos θ' square meters. For example, if θ' = π /2 at sunrise, then zero power is deli vered to that square meter of ground. Our patch of ground does not care about the azimuth of the sun, only its elevation ( π/2-θ'). If we assume a surface solar power flux of S watt/m2, then the solar energy delivered into 1 square meter of ground during one day (assume the ground is perfectly black and reflects nothing) is given by insolation/S = 2 ∫0 tss cosθ' dt = 2 ∫0 tss [sinθt cosθ1 + cosθt sinθ1 cos(ωt)] dt = s i n θt cosθ1 (2tss) + 2 cosθt sinθ1 ∫0 tss cos(ω t)dt = s i n θt cosθ1 (2tss) + 2 cosθt sinθ1 (1/ω) sin(ωt)|tss 0 = s i n θ t cosθ1 (2tss) + (2/ω) cosθt sinθ1 sin(ωtss) (10.7) Section 10: Insolation 91 where tss is given in (10.2). In terms of site latitude θLAT = π/2-θ1, insolation = S [ sin θt sinθLAT (2tss) + (2/ω) cosθt cosθLAT sin(ωtss) ] watts-hrs/m2/day (10.8) Daily Insolation plotted versus latitude The following is a plot of this daily insolation vers us latitude at three times of the year (S = 1 watt/m2): Fig 10.6 where green, red and blue are for December solstice, equinox and June solstice. The kinks in the solstice curves are due to the kinks shown in t ss at the pinning points in Fig 10.3. The blue June curve shows that as one moves north from the equator, insolation peaks around 41o, then starts to drop, but then rises again in the all-daylight region to have a final peak at the north pole (where the sun is up all the time). The units shown on this graph are watt-hours/m2/day assuming S = 1 watt/m2 from the sun. As was noted in Section 1 (a), the solar flux is in fact 6.9% greater at summer solstice in the southern- hemisphere compared to the northern due to the elliptical earth orbit. The annual average value of S is about 1.367 kW/m2, so we can then make this more practical version of the above plot where the vertical axis now shows kW-hrs/m2 per day: Section 10: Insolation 92 F i g 1 0 . 7 which may be compared with this plot found on the web, Fig 10.8 http://www.applet-magic.com/insolation.htm Section 10: Insolation 93 The 3D Daily Insolation function plo tted versus latitude and day of year Here is a 3D plot showing the shape of the insolation function versus latitude and day of year (back to S = 1). Two views of the same plot are shown. F i g 1 0 . 9 The front face of the cube on the left is the blue curve above, while a slice ha lfway back gives the green curve. Daily Insolation plotted versus day of year Slices parallel to the left face of the view on the left above (front face of the view on the right) give daily insolation versus day of year for various latitudes, wh ich we shall now plot. Rather than express the result in watt-hrs/m2 per day, it is useful to express things in terms of watt-hrs/m2 per average hour of the day, so the units are then watt-hrs/m2/hr = watts/m2. In other words, this is the average power/m2 during the day and the number is then easily compared with the incident power flux 1367 watts/m2. Section 10: Insolation 94 F i g 1 0 . 1 0 An interesting fact here is that within 10o of the equator, daily insolation is not at a maximum at summer solstice. In fact, at the equator in solation is the same at summer and winter solstice and is at its lowest level of the year; it is maximal at the equinoxes. Again, Fig 3.13 and 3.14 show why this is so. Negate the latitude numbers for the southern hemisphere. Below is a web plot showing the same information for latitude 0,30,60 and 90. There is some small discrepancy in the vertical scale; th e author does not say what value of S was used to make the plot. Fig 10.11 http://www.physicalgeography.net/fundamentals/6i.html Section 10: Insolation 95 (c) Annual Insolation plotted versus latitude Maple is happ y to do some numerical integrations over the year of the daily insolation at selected latitudes. These integrals can be thought of as areas under the curves shown in Fig 10.10. The numbers are in kW-hrs/m2 per year: From this data we can create a normalized plot of annual insolation versus latit ude. First we construct a left side and right side sequence of points, Then we plot the curve with gridlines, Section 10: Insolation 96 F i g 1 0 . 1 2 The curve peaks at the equator as one might e xpect. (Mean temperature has a similar shape.) The above graph may be compared to this web graph which was formed by averaging insolation at the two solstices and the two equinoxes: Fig 10.13 http://www.applet-magic.com/insolation.htm Section 11: The Solar Analemma 97 11. The Solar Analemma (a) The Sun Sphere and the Analemma In this sectio n the word "ecliptic" means the path of the sun on some surface. On the celestial sphere, the ecliptic is the red circle shown in Fig 3.5. We imag ine this red circle is inscribed or burned onto the celestial sphere by the moving sun, and then we can late r view the entire circle. In this context, we shall show below that "the solar analemma is the sun's ecliptic on the sun sphere". What we call "the sun sphere" is a mental construct similar to, and aligned with, the celestial sphere of Section 3. The celestial sphere, enclosing the rate- ω s rotating earth, is at rest relative to the stars. The sun sphere, being a similar exterior cage affair having lines of longitude and latitude, rotates once per year (rate Ω0) relative to the distant stars (and the celestial s phere). Thus, the annual paths of the distant stars on the sun sphere are in fact white circles of latitude (declination). The fictitious sun, however, is at rest on the sun sphere! Recall that this fictitious sun executes a uniform annual path around the equator of the celestial sphere, having therefore an ecliptic equal to that equatorial circle. But the sun sphere rotates in the same direction and at the same speed as the fictitious sun, causing that sun to be at rest on it. If we now hop on board the sun sphere as a frame of reference, here is what we see: (1) the fictitious sun at rest at a point on the sun sphere equator; (2) the earth rotating inside the sun sphere at rate ω = ω s – Ω0, which is to say, one rotation each 24 hours. (3) we also see the actual sun following some closed path on the surface of the sun sphere during the year, gene rally in the neighborhood of the fixed point of the fictitious sun. That closed path is called "the solar analemma"; the path is the actual sun's ecliptic on the sun sphere. We can and should imagine this analemma to be in scribed on the surface of the sun sphere. We could photograph this inscribed trail from the earth at any time after it was created. In practice, however, the sun leaves no such trail on the sun sphere, so in orde r to see the analemma we can take a series of 365 multiple-exposure, stroboscopic, heavily filtered photo s of the sun. These exposures are taken at the exact same local mean time every day from a rigidly mounted camera on the earth that is aimed generally in the known direction of the fictitious sun at that time of day. The camera film will then contain 365 white dots which trace out the continuous analemma locus. The fi lm is then printed to obtain a photograph. A good representation of the analemma could be obt ained by taking a photo perhaps every 10 th day. Of course there are cloudy days and snowstorms and othe r disturbances (earthquakes!), so photographing the analemma is a challenging yet doable task. In the strobe light of this camera stroboscope, the sun sphere is at rest, so in effect we are taking pictures of a fixed patch on the sun sphere. From a reference frame attached to the earth, the sun sphere rotates around the earth at rate - ω. The sun then appears to be executing a spiral path similar to that shown in Fig 3.7. The strobe camera samples this path once per rotation, and the strobed image is the analemma. As we shall see, the analemma has r oughly the shape of a skinny bowling pin Section 11: The Solar Analemma 98 Bowling pin Analemma Fig 11.1 Sometimes analemmas are pl otted with the angle axes having different scales, but in the plot above the two scales are the same. The vertical extent of the analemma is caused by the earth's tilt. Looking at the actual sun ecliptic on the celestial sphere in Fig 3.5 or Fig 6. 1, one sees that the sun varies between -23 o to 23o relative to the celestial equator, so we expect the angular height of the bowling pin analemma on the sun sphere to be 46o, which is a fairly large angle. The horizontal extent is due to the azimuthal vari ation of the actual sun's position from that of the fictitious sun. This difference is the E quation of Time (6.15). Here we plot Δφ in degrees using the code shown above (6.4) but with a scaling factor (180/ π) and with n = t Me to diy + t Me with tMe from (4.11), Fig 11.2 The plot shows a peak-to-peak deviation of about 7. 5 degrees (green line is the fictitious sun). The left end of this graph is τ = 0, the March equinox, and corresponds to th e gap we have intentionally left in the analemma drawing above. The analemma path goes up to the northeast from this gap. The narrower upper part of the bowling pin corresponds to the two smalle r bumps on the left of the Equation of Time in the first half of the year, and then the wider bottom part is due to the two larger bumps which end the year. The shape of the equation of time curve results from a combination of the effect of the earth's tilt on the sun's azimuthal speed, and the effect of the earth's orbital eccentricity on that same speed, as was shown in Fig 6.5. Since each pair of bumps is not perfectly symmetric, the analemma is not left/right symmetric though it appears that way above (see Figs 11.18 and 11.26 below). Section 11: The Solar Analemma 99 (b) Qualitative appearance of the inscribed analemma in photos As just explained, the an alemma is aligned vertically along that sun-sphere line of longitude which contains the fictitious sun. In this snapshot, taken from a foggy and snow less north pole, we have turned off the sun, and we are viewing the previously "i nscribed" analemma on the sun sphere. We can only see that part of the analemma which lies above the horiz on plane, but we show the full analemma anyway (Maple, Stellarium Visio and Irfan) Fig 11.3 The gridlines are those of the sun sphere. Notice that all the lines of longitude meet the horizon plane vertically, but bend toward the center going up since they must meet directly overhead of the viewer. The above picture will look the same regardless of the lo cal mean time of the photo. One could wait an hour, the earth would rotate 15 o, and then one could take another photo aiming the camera again at the inscribed analemma. The analemma will look exactly the same as it did in the first photo. In other words, the inscribed analemma just moves in azimuth as the earth rotates. Now what does this sa me analemma look like around 8AM local mean time when viewed from 40o north latitude on the earth, instead of from the north pole? The picture below is approximate, but demonstrates the basic facts. Section 11: The Solar Analemma 100 Fig 11.4 For this viewer, the "cage" of the sun sphere has latitude and longitude lines which appear as above. The sun-sphere equator is marked top center (it hits the letter E). The north pole of this cage is at the North star, behind the viewer and up in the sky. The anale mma is still lying along its same sun-sphere line of longitude, but that line of longitude is tipped to the left due to the new observer location and the ~8 AM time of day. If the observer waits 4 hours till noon, the inscribed an alemma will move to the right and will become vertical at noon, lining up with the local meridian. At 4 PM in the afternoon the analemma will be off to the right and will be tipped to the right, a mirror imag e of the above drawing, reflected through the right edge of the picture which is the meridian. Here is a genuine 34-shot strobed analemma which we will use below in section (g) as a case study (photographer's latitude 44.7 o ; shot at 8 AM local mean time; credit and reference given there) Fig 11.5 Section 11: The Solar Analemma 101 For the observer in Fig 11.4, the analemma can never be horizontal because its longitude line is not visible when horizontal. If the viewer moves to the e quator, however, then at 6 AM local mean time the analemma will be exactly horizontal, F i g 1 1 . 6 The qualitative images of the imaginary inscribed analemma shown above are those that would appear on the film of a camera. In the following section, we shall compute the exact shape of the analemma in "angle space" and not in "camera film space". In angl e space, we shall find that the analemma pattern bends away from the center as we move away from local noon. Later, however, we shall compute the image of the analemma as it appears on film ("screen space") and this bending of the angle-space analemma goes away, and the approximate images above are realized. The exact shape and tilt of the film analemma depends, due to perspective effects, on the exact direction in which the camera is pointed. (c) Calculation of the analemma in angle space In this section we develop the equations used in our Maple program to compute the analemma. The math is described here, while the program itself is described in section (e) below. The notion of local mean time LMT was presented above in Section 8 (b). Section 11: The Solar Analemma 102 1. Location of the actual sun and the mean fictitious sun The location of the actual sun for our Frame S' observer at polar angle θ1 is given by θ' = cos -1[sinθt cosθ1 + cosθt sinθ1 cos(φb-φ) ] (7.29) φ' = tan-1 ( [sinθt sinθ1 – cosθt cosθ1 cos(φb –φ) ] / [– cos θt sin(φb –φ)] ) (7.29) Φ' = π + tan-1 ([cosθt sin(φb-φ)] / [– sinθt sinθ1 + cosθt cosθ1 cos(φb-φ) ] ) (8.3) θ t = sin-1[sinθtmax sin(ψ-α) ] (6.13) φ = tan-1[cosθtmax tan(ψ-α) ] (6.13) ψ = Ω0t + 2e sin( Ω0t) + (5/4)e2 sin(2Ω0t), t = τ + tMe (6.13), (7.2) φb = φb0 + ωsτ with φb0 = φ1 + φG0 . (7.3),(7.7) (11.1) Angles θ',φ' are the spherical angles of the actual sun in Frame S' as in Fig 7.3, while Φ' = π/2-φ' as in (8.2). Angles θt and φ are the declination and right ascension of the actual sun in Frame S as in Fig 7.1. Angle φ1 is the east longitude of the observer site relative to Greenwich, and φ G0 is the right ascension that Greenwich had at the last March equinox, τ = 0. Angle θtmax is the famous 23o tilt (obliquity). Rate Ω0 = 2π radians per year, and ω s = ω + Ω0 is the sidereal rotation rate of the earth as in (1.6), where ω = 2π/24 radians per hour is the solar rotation rate of the earth. The dimensionless e is the eccentricity of the earth's orbit, while α = π/2 - φper where φper is the 13o perihelion angle shown in Fig 1.3. Time t Me is the time of March equinox relative to perihelion t = 0 as in (4.12). For the fictitious sun, having e = 0 and θtmax = θt = 0, these equations reduce to θ' f = cos-1[sinθ1 cos(φb-φf)] φ'f = tan-1 ( [–cosθ1 cos(φb –φf) ] / [–sin( φb –φf)] ) Φ'f = π + tan-1 ( [sin(φb-φf)] / [cosθ1 cos(φb-φf) ] ) . (8.3) θ t,f = 0 φf = Ω0t - α = Ω0τ + k (7.17) k = [Ω0tMe - α] (7.16) φb = φb0 + ωsτ with φb0 = φ1 + φG0 . (7.3),(7.7) (φb - φf) = (φb0 + ωsτ) – (Ω0τ + k) = ωτ + (φb0- k) (1.6) (11.2) The difference Δφ ≡ φ - φ f is then obtained by subtracting φf above from φ in (6.13) above, and the result matches (6.15a), Δφ = – Ω0τ – k + tan-1[cosθtmax tan(ψ-α) ] (6.15a) (11.3) This function is called "the equati on of time" and has this characteristic shape in degrees instead of radians (see Fig 6.5 which has t = 0 as its left edge; the figure here has τ = 0 as its left edge), Section 11: The Solar Analemma 103 Fig 11.2 According to (11.2) and the fact that Δφ ≡ φ - φf we get (φ b-φf) = ωτ + (φb0- k) (φb-φ) = (φb-φf) - Δφ = ωτ + (φb0- k) - Δφ ( 1 1 . 4 ) and these phases may be inserted into expressions above as needed. 2. Relationship between local time LMT and master time τ : shooting tim es τN In this section, ω = π/12 rad/year, and time τ is measured in hours. Recall from Section 8 that : LMT = 12(Φ ' fNP/π) (8.9) where this fictitious-sun-tracking No rth Pole azimuth has the form ωτ + constant, Φf'NP = { ωτ + π + φb0 - k } mod 2 π . (8.7) Here τ is our "master time" which is 0 at the last March equinox. Then, LMT = (12/π ) Φ'fNP = (12/π) { ωτ + π + φb0 - k } mod 2 π = { (12/π ) [ωτ + π + φ b0 - k] } mod [(12/π )2π] = { (12/π ) [(π/12)τ + π + φ b0 - k] } mod 24 = { τ + 12 + (12/ π)(φ b0 - k) } mod 24 (11.5) If τ N are the times of analemma daily phot o shoots, we can solve (11.5) for τN in terms of day N and the LMT at which those shoots are taken, { τN + 12 + (12/ π)(φb0 - k } = LMT + 24N, N = integer τ N = [ LMT – 12 – (12/ π)(φb0 - k)] + 24N, N = 0,1,2.... Section 11: The Solar Analemma 104 which is to say τN = τ0 + 24N, N = 0,1,2.... τ0 = [ LMT – 12 – (12/ π)(φb0- k)] (11.6) In other words, N corresponds to the number of the day of the N th photograph. N = 0 marks a day in the vicinity of the March equinox. Notice that, with ω = π/12 rad/hr, ωτ N = ωτ0 + N(2π) => ωτN mod 2π = ωτ0 mod 2π sin(ωτN + a) = sin( ωτ0 + a) cos( ωτN + a) = cos(ωτ 0 + a) (11.7) 3. Location of the actual and fictitious strobed suns in terms of LMT The phase appearing in the first three equations of (11.1) may be written, (φb-φ)N = (ωτN + (φb0- k) - Δφ(τN)) = (ωτ0 + (φb0- k) - Δφ (τN)) // (11.2) and (11.7) = (ω [ LMT – 12 – (12/ π)(φb0-k)] + (φb0- k) - Δφ (tN)) // (11.6) and ω = π/12 = (ω [ LMT – 12] – ( φ b0-k) + (φb0- k) - Δφ(τN)) = (ω [LMT- 12] - Δφ N) . ( 1 1 . 8 ) We then arrive at this strobed ve rsion of (11.1) and (11.3), with ω = 12/π and τ in hours, θ'N = cos-1[sinθtN cosθ1 + cosθtN sinθ1 cos(ω[LMT- 12] - ΔφN)] Φ'N = π + tan-1( [cosθtN sin(ω[LMT- 12] - ΔφN)] [– sinθtN sinθ1 + cosθtN cosθ1 cos(ω [LMT- 12] - ΔφN) ] ) θtN = sin-1[sinθtmax sin(ψN-α)] ψN = Ω0tN + 2e sin( Ω0tN) + (5/4)e2 sin(2Ω0tN), tN = τN(yrs) + tMe (6.13) Δφ N = – Ω0τN(yrs) – k + tan-1[cosθtmax tan(ψN-α)] τ N(yrs) = [τ0 + 24N] / (diy*24) N = 0,1,2...364 τ0 = [ LMT – 12 – (12/ π)(φb0- k ) ] ( 1 1 . 9 ) and for the fictitious sun, this strobed version of (11.2), θ' fN = cos-1[sinθ1 cos (ω[LMT- 12] )] Φ'fN = π + tan-1 ([sin (ω[LMT- 12])] / [cos θ1 cos (ω[LMT- 12]) ]) . (11.10) As expected, the strobed fictitious sun position is independent of strobe time τ N. Its analemma is a dot. Equations (11.9) lend themselves to a hand calculati on of the analemma, but in our Maple code below we shall instead use the continuous-time equations ( 11.1) and create data arrays by sampling (strobing) Section 11: The Solar Analemma 105 these equations with τ = τN from (11.9). In that code ω is expressed in radians/year, so the τN strobing times appearing there are τ N(yrs). 4. View of the continuous θ ',Φ' data before strobin g One should remember that during the year, the functions θ' and Φ' of (11.1) [ the position of the sun as seen from Frame S' ] oscillate every day. Here is a typical view of θ ' during a year, with a zoomed section below (the dark red areas in the upper image are just an aliasing artifact) Fig 11.7 The daily strobing picks out a slow sinusoidal curve which tracks the envelope. In these plots, the value θ' = 1.57 = π/2 is that of the points of daily sunrise and sunset. In the first plot, the sun reaches maximum elevation (smallest θ', lowest excursion of the bottom of the envelope) at June solstice located at τ = tJs - tMe = .4623 - .2084 = .2539, using numbers from (4.11). Meanwhile, Φ ' is even more violent in its activity, Section 11: The Solar Analemma 106 Fig 11.8 During each day, angle Φ' of the observed sun increases from 0 to 2 π and then jerks back to 0 again because the value is always shifted to the range (0,2π ). The jerk point is Φ' = 0 which is midnight for a northern hemisphere observer, so we don't need to worry about our solar photographer strobing right on those vertical lines. (d) Calculation of the analemma in film-space The equation s developed in section (c) above are used in our Maple program to compute the angle-space analemma. Here we show a method (installed at th e end of the same Maple program) for computing the analemma as it appears on film. 1. Overview In co mputer graphics language, one has "the scene" whose point s have Cartesian coordinates r = x,y,z in some Frame S. As viewed by our Frame S' observer, points in the scene have Cartesian coordinates r' = x',y',z' and spherical coordinates R, θ',φ'. Our analemma observer has a cam era which he points in some direction generally not aligned with the axes of Frame S'. That camera then defines a Frame S c in which point r' now has coordinates rc= xc, yc, zc . This Frame S c is often called "the eye coordinate system" and the transformation from r' to rc is called "the eye transformation", but we shall replace the word eye with the word camera. The unit vector z^c is the direction in which the camera is pointing, x^c is "to the right", and y^c is "up". We can imagine the camera viewfinder providing a "rectangular view" of the scene. This view is a chunk of solid angle bordered by four planes. The intersection of this solid angle with a Section 11: The Solar Analemma 107 plane whose normal is z^c and which is some arbitrary distance out in front of the camera is a rectangle called "the viewport". Each 3D point in the scene r c= xc, yc, zc maps to some 2D point on the viewport, rv = xv,yv. The final step is to transform these viewport coor dinates to coordinates on our piece of film, or to coordinates on a photo printed from this film. This last step is called "the screen transformation" and we end up with rs = xs, ys. In the case of a camera, the lens might introduce some distortion in going from the viewport to screen coordinates, but we shall assume this does not happen. The process just described runs on millions of electronic devices on a daily basis. We now carry out the steps of the program outlined above. 2. The transformation from Frame S to Frame S' The first step is already done. Equation (7.14) provides the Cartesian coordinates of points in Frame S', x' = R cos θt sin(φ –φb) y' = R [sin θt sinθ1 – cosθt cosθ1 cos(φ –φb) ] z' = R [sin θt cosθ1 + cosθt sinθ1 cos(φ-φb) ] (7.14) where (φ b-φ) = (φb-φf) - Δφ = ωτ + (φb0- k) - Δφ , (11.2) and R is the radius of the celestial sphere. These points are of course positions of the sun on that sphere. 3. The transformation from Frame S' to Frame S c of the camera (eye transformation) The next step is geared to a northern hemisphere observer, and we recall Fig 7.3 to the stand, Fig 7.3 For a southern hemisphere observer, we replace do south ↔north, east↔west, and Φ' → Φ's. Section 11: The Solar Analemma 108 Imagine a camera at the origin of the above figure which is facing south. For this camera we have z^c ≡ straight ahead to the south = - e'2 = - y^' y^c ≡ up = e'3 = z^' x^c ≡ to the right = west = - e'1 = - x^' . ( 1 1 . 1 1 ) Imagine that the point marked "sun" is the location of the fictitious sun at some LMT for the observer. We want this point to appear at the center of the camera image. Rather than rotate the camera, we will move the sun to achieve our goal. We have to do this set of rotations on the sun point to get the sun at camera image-finder center: r' new = Re'1(π/2-θ'f ) Re'3( -[π-Φ'f]) r'fictsun . (11.12) For reasons to be given in a moment, let us define these new variables αt ≡ π/2 - θ 'f αp ≡ Φ'f - π ( 1 1 . 1 3 ) so that the sun transformation (11.12) becomes [ r' new = (x'new, y'new, z'new) ] r' new = R r'fictsun R ≡ R1(αt) R3(αp) . ( 1 1 . 1 4 ) Then, given the usual active rotation matrices, R1(θ) = ⎝⎜⎛ ⎠⎟⎞ 1 0 0 0 cosθ -sinθ 0 sinθ cosθ R2(θ) = ⎝⎜⎛ ⎠⎟⎞ cosθ 0 sinθ 0 1 0 -sinθ 0 cosθ R3(θ) = ⎝⎜⎛ ⎠⎟⎞ cosθ -sinθ 0 sinθ cosθ 0 0 0 1 , (11.15) we find that, with help from Maple, R ≡ R 1(αt) R3(αp) = ⎝⎜⎜⎛ ⎠⎟⎟⎞ 1 0 0 0 cosαt -sinαt 0 sinαt cosαt ⎝⎜⎜⎛ ⎠⎟⎟⎞ cosαp -sinαp 0 sinαp cosαp 0 0 0 1 = ⎝⎜⎜⎛ ⎠⎟⎟⎞ cosαp -sin αp 0 cosαtsinαp cosαtcosαp -sinαt sinαtsinαp sinαtcosαp cosαt . ( 1 1 . 1 6 ) Then according to the first line of (11.14) we have x' new = cosαpx' - sinαpy' y'new = cosαtsinαp x' + cosαtcosαpy' -sinαt z' z'new = sinαtsinαpx' + sinαtcosαpy' + cosαt z' . (11.17) Now consider this vector r'new expanded two ways Section 11: The Solar Analemma 109 r'new = x'newx^' + y'newy^' + z'newz^' = x cx^c + ycy^c + zcz^c ( 1 1 . 1 8 ) On the first line we replace the unit vectors as shown in (11.11) to get r' new = - x'newx^c - y'newz^c + z'newy^c ( 1 1 . 1 9 ) Comparing the last two equations we must have x c = - x'new, yc = z'new and zc = - y'new. Thus, x c = - cosαpx' + sinαpy' yc = sinαtsinαpx' + sinαtcosαpy' + cosαt z' zc = -cosαtsinαp x' - cosαtcosαpy' + sinαt z' (11.20) and this then is the "eye transformation" taking points from Frame S' to Frame S c of the camera. Active vs Passive: camera pan and tilt In the above development we moved the sun to a ca mera-centered position using rotation R. As shown in Appendix B, the same result is obtained (sun at came ra center) if the camera (and its rigidly attached axes) are instead rotated by R -1 and the sun is left where it is in the sky. The camera rotation is then R-1 = R3(-αp) R1(-αt) = Rz'(-αp) Rx'(-αt) . (11.21) Looking at Fig 7.3 just above, and using the right hand ru le, this says that the camera is first "tilted up" by angle αt , and it is then "panned" clockwise by αp . If we imagine the camera attached to some kind of polar angle mount, it is this mount that is in fact "panned", not the camera itself, but we shall still refer to this as a "pan" action. So αt and αp are the camera tilt and pan paramete rs. It is the equations (11.13) which cause the tilted then panned camera to point to th e sun at film center, but we could then imagine adjusting αt and αp slightly to make the camera point to some location near the sun, or perhaps to some other location altogether. In the particul ar case illustrated by Fig 7.3, perhaps αt = 50o and αp = -120o. This model for the operation of the camera platform h as implications for the way an analemma appears on film, as discussed in Appendix D. 4. The transformation from Frame S c to the camera film (viewport transformation) Since we have only one screen in mind (as opposed to a different screen for each different TV monitor resolution, say) we let the viewport be the screen and then we do the viewport transformation and dispense with the screen transformation. The viewport tr ansformation is easily obtained from this picture, Section 11: The Solar Analemma 110 F i g 1 1 . 9 where one point of interest rc is shown on the celestial sphere of radi us R, and where the "viewport" is positioned some arbitrary distance z 0 to the right of the camera marked "eye". The viewport transformation taking us to our screen coordinates is then (by similar triangles y c/zc = ys/z0), x s = z0(xc/zc) ys = z0(yc/zc) . ( 1 1 . 2 2 ) This shows the famous "perspective divide" (by z c for each point) that used to be very expensive to do in computer hardware, but which has become cheap thanks to Moore's Law. If angle α were 450, one would have yc = zc and then y s = z0. In the Maple code we set z 0 = 1. Comments : In the above picture, regardless of the tilt a nd pan of the camera relative to "south", the camera points to a scene-center point C, distance R away, which lies on a spherical surface (celestial sphere) whose local region is perpendicular to the came ra's pointing direction. The picture above is then obtained by "rolling" the camera about its axis so the point r c = (xc,yc) lies in the plane of paper defined by xc = 0. It is not assumed that the sphere is "flat" in the region of interest, so the screen transformation is accurate even if the angle α is large, as it in fact is in an analemma situation. 5. Putting the pieces together The Fram e S' positions of the sun are given by x' = R cos θt sin(φ –φb) y' = R [sin θt sinθ1 – cosθt cosθ1 cos(φ –φb) ] z' = R [sin θt cosθ1 + cosθt sinθ1 cos(φ-φb) ] (7.14) (φb-φ) = (φb-φf) - Δφ = ωτ + (φb0- k) - Δφ . (11.2) The camera coordinates of the sun position are x c = -cosαpx' + sinαpy' yc = sinαt sinαp x' + sinα t cosαp y' + cosαt z' zc = -cosαt sinαp x' - cosαt cosαp y' + sinαt z ' (11.20) Section 11: The Solar Analemma 111 and the film coordinates are x s = zo(xc/zc) ys = zo(yc/zc) . (11.22) We shall point the camera to the static location of the fictitious sun since this point is roughly at the center of the analemma. Recall that this point is given by (11.10) θ' f = cos-1[sinθ1 cos (ω[LMT- 12] )] Φ'f = π + tan-1 ([sin (ω[LMT- 12])] / [cos θ1 cos (ω [LMT- 12]) ]) (11.10) where ω = 2π/24 radians per hour. As noted in (11.13) and the discussion above, the parameters α t and αp in (11.20) are the camera tilt and pan parameters where αt = π/2 - θ 'f // tilt up to strobe-time elevation of the fictitious sun αp = Φ'f - π // pan clockwise to the strobe-time location of the fictitious sun (11.23) The Maple code allows for adjustments to these a ngles to accommodate a camera which is pointed at some other point near the analemma. As the battery of equations above suggests, the appearance of the analemma on the film will be dependent on the direction in which the camera is pointed. Section 11: The Solar Analemma 112 (e) A tour of Maple code which comput es angle-space and scre en-space analemmas The Maple code is freely available from the author upon request. The code appearing below is in screen clips to improve clarity and save time. In Maple, # is the comment character and % refers to "the last thing computed". A command ending with ; displays its output, one ending with : suppresses its output. The parameters show n here are for our 1999 Ukraine "case study" discussed in sec tion (g) below. In the first block of code, Maple is initialized, a plotting package is activated, a nd then constants are entered Next, some one-time calculations are carried out, Section 11: The Solar Analemma 113 A procedure to do the mod 2 π operation is specified, with some "noise reduction" installed, and then comes the arctan2Pi routine described in detail in Appendix A, Equation (4.10) is written as a function and then t Me and k area computed from that function, Section 11: The Solar Analemma 114 Equations for the position of the sun in Frame S are then entered ( θ and φ with θt = π/2 - θ ), Next, equations for the position of the sun in Frame S' are entered ( θ' and φ' with Φ' = π/2 – φ '), The photographer's local mean time LMT is then entered, and the function tstrobe(N) is defined, The analemma data arrays are then gene rated, sampling the continuous functions θ' and Φ' using tstrobe, Section 11: The Solar Analemma 115 The first two plots are θ' and Φ ' (left edges are τ = 0, March equinox) The next plot is the analemma in angle space, where Φ' is the horizontal axis and θEL = 90-θ ' the vertical: Fig 11.10 The final section of the code computes the anale mma in the film "screen space". as summarized in Section 11 d (5). The user can enter tilt and pan off sets to move the camera's pointing direction away from the point of the mean fictitious sun. Since phot ographers generally don't doc ument the exact point at which their analemma camera was pointed, and since this is sometimes further confused by cropping and/or compositing in a background photo after the fact, some trial and error might be needed to "fit" the analemma computed by the Maple program to a publ ished photo. For Mr. Rumyantsev's photo we use a pan adjust of -6o as indicated by the cross hair location on the plot below. This block of code allows entr y of these offsets and computes the αt and αp tilt and pan parameters, Section 11: The Solar Analemma 116 The final block of code then computes and plots the analemma as it appears in a photo, again as summarized in Section 11 d (5). The viewport distance is set z 0 = 1 unit from the camera : Analemma on the multiple-exposure film Fig 11.11 If this blue screen-space analemma (with no pan offset) is superposed on the red angle-space analemma without altering the aspect ratio of either image, we find Section 11: The Solar Analemma 117 Fig 11.12 The blue screen-space analemma reliably represents the near symmetry of th e analemma along its long axis, but the red angle-space analemma is bent as shown. Here is a very crude explanati on of this bending. The drawing be low shows a slice of the celestial sphere near the pole for emphasis : Fig 11.13 The blue region is a patch of physical space (a patc h of the sky) on the slice between two blue lines of right ascension. When this region is mapped to angle space (red rectangle), the physical blue region is stretched as indicated by the arrows. In the middle of th e blue region is a straight blue line, and this gets warped into the red curve. The bl ue line represents the nearly-strai ght analemma photo, while the red curve represents the angle-space analemma. Section 11: The Solar Analemma 118 (f) An Analemma Gallery 1. Analemmas for an observer at the North Pole These curves are made for the blue, black and red earths discussed in Section 6 (d). The red earth is the real earth, the blue earth's sun has a tilted but circul ar orbit, while the black earth's sun has an elliptical orbit but no tilt. At the north pole, θ' and Φ' are simply related to θt and Δφ as follows, θ' = θ = π/2-θt => θ'EL = π/2-θ' = θt // north pole Φ' = π + ωτ - (k-φb0) - Δφ . // north pole (11.24) The first equation θ' = θ is obvious since at the north pole, the polar angle of the sun is the same in Frame S' as in Frame S. It can also be shown from (7.29) with θ1 = 0. The second equation comes from (8.5) which says Φ 'NP = π + (φb-φ), from Δφ ≡ φ-φf, and then from the last line of (11.2). The main idea of this equation is that ωτN = ωτ0 mod 2π so under daily strobing, Φ'N = constant - ΔφN, so Φ 'N maps out the equation of time which gives the horizontal extent of the analemma of Fig 11.14 below. Blue Earth Analemmas at the North Pole (eccentricity = 0) When e = 0, (4.12) shows that t Me = α/Ω0 and then (7.16) says k = 0, so these plots have e = k = 0. 12 noon e = k = 0. Fig 11.14 This is the "prototype analemma" and we have much to say about it that won't be repeated for later analemmas. A blowup is shown on the right (see also labeling in Fig 11.1). Again, this is the "blue earth" situation of S ection 6 (d) and (e). The plots in the middle are θ' and Φ', where the left edge of the plots are at March equinox, τ = 0. For the analemma plot, Φ' is the horizontal axis in degrees, and elevation = θ EL = π/2 – θ' = θt is plotted vertically in degrees. Both axes are shown with the same scale. Looking south at noon every day for a year, the sun will trace out the blue analemma, but only the top half is visible above the horizon. The topmost point is the day of June solstice (the day elevation is the highest), the bottom is December solstice, and the horizon is the equinox. The path starts Section 11: The Solar Analemma 119 at March equinox on the horizontal axis. It then goes up and to the right because Φ' is moving positive. The days then trace out the "8" pattern like a car on a racetrack of this shape. This blue earth analemma is symmetric in both directions. Although this seems obvious from the plots shown of θ' and Φ', and although it perhaps seems intuitiv e, actually showing this symmetry from the analemma equations requires a bit of work, and that work has been relegated to Appendix C. The analemma is centered in the south at 180o because the photos are taken at noon and the fictitious sun is at Φ' = 180o every day at noon. The analemma's vertical extent above the horizon is 23o. Its horizontal extent is about 5o. If we take photos in the morning instead of at noon, the analemma retains its shape but shifts to the left, as is reasonable since it is morning. From (8.6) and (8.7) we know that Φ ' = π + ωτ - (k-φ b0) - Δφ , so if we do morning shots instead of noon shots, ωτ is a smaller number at each τ = τN and Φ'N moves to the left. Since Δφ is a slow function of τ, it does not change much in going to a morning strobe time. Here then is the 9 AM analemma, north pole 9 AM e = k = 0. Fig 11.15 Of course the afternoon-shot analemma will be shifted to the right of 180o instead of to the left. Black Earth Analemmas at the North Pole (tilt = 0) The black earth has no tilt, so the an alemma has no vertical extent. For example, at noon the analemma is a horizontal locus centered at 180 o. The width of this locus is about 4o (see black curve in Fig 6.5). north pole 12 noon θtmax= 0. Fig 11.16 Morning or afternoon photography again just shifts this locus left or right. Section 11: The Solar Analemma 120 Red Earth Analemmas at the North Pole (actual earth) We now have the full equation of time and here is the noon analemma (blowup on right) north pole 12 noon actual Fig 11.17 The equation of time (- Δφ) has less amplitude in the summer than in the winter as the red Φ' graph shows, so the analemma is narrower on the top than on the bottom. Again, only the upper part is visible. The vertical extent is still 23 0 on the visible part. The horizontal extent is now about 7o on the (non-visible) bottom and 2.50 on the top (see red curve Fig 6.5). Again, morning or afternoon photography agai n just shifts this locus left or right. Allowing the horizontal axis to have a different scale from the vertical, the plot on the right above has this appearance, where the left/right asymmetry is more visible, Fig 11.18 Section 11: The Solar Analemma 121 2. Analemmas for observers at other lati tudes in the northern hemisphere. Blue Earth Analemmas at the Latitude 70o (eccentricity = 0) Here is the noon analemma. The upper and lower l obes are different but not too different, and the analemma his lifted almost clear of the horizon: latitude = 70o 12 noon e = k = 0 Fig 11.19 Now a new effect appears. If the photos are shot in the morning or afternoon, the analemma shifts horizontally of course, but now the top tilts away from 180o latitude = 70 o 9 AM e = k = 0 Fig 11.20 Section 11: The Solar Analemma 122 latitude = 70o 3 PM e = k = 0 Fig 11.21 In these last cases, Appendix D Fig D.2 says the tipping angles away from vertical are about ± 15o for the film analemmas, which roughly align wi th the angle space ones shown above. Red Earth Analemmas at the Latitude 70 o (actual) latitude = 70 o 12 noon actual Fig 11.22 The analemma now looks like our bowling pin of Fig 11.1, it is more asymmetrical vertically. Morning and afternoon shooting causes horizontal shift and tilt of the top away from 180o just as in the blue case. Red Earth Analemmas at the Latitude 40o (actual) The change here is that the analemma is shifted to a higher elevation in the sky, Section 11: The Solar Analemma 123 latitude = 40o 12 noon actual Fig 11.23 latitude = 40o 9 AM actual Fig 11.24 latitude = 40o 3 PM actual Fig 11.25a Section 11: The Solar Analemma 124 Now the angle-space anale mmas in morning or afternoon are shifted, tilted, and bent away from the 1800 line. Figure D.2 gives the film analemma tipping angles as ± 40o away from vertical. Here is the three corresponding film analemmas: Fig 11.25b Red Earth Analemma at the Latitude 51.4791o Our computed analemma is on the left. On the right it is superposed on a red analemma which appears on the wiki page http://en.wikipedia.org/wiki/Analemma . F i g 1 1 . 2 6 Section 11: The Solar Analemma 125 Red Earth Analemmas at the Latitude 24.5o latitude = 24,5o 12 noon actual Fig 11.27 Now the top of the analemma in our ( Φ', 90-θ') angle space experiences a horizontal stretching at the top similar to the "Greenland Effect" on a Mercator projection map. The camera film analemma is shown in blue on the right and does not exhibit the angle spac e distortion. The distortion is worst at noon because that drives the upper end of the analemma closer to 90o elevation. Red Earth Analemmas at the Latitude 20o We have entered the tropical zone, and the a ngle-space analemma pulls apart at the top : latitude = 20o 12 noon actual Fig 11.28 The above plots are completely accurate in angle space. First look at the black θ' plot on the right. Starting at March equinox, the strobed sun position go es to higher elevation and lower polar angle θ' until it reaches polar angle 0. At that point it goes "over the top", the polar angle starts increasing again, and at Section 11: The Solar Analemma 126 some point the azimuth Φ' crosses our Fig 7.3 discontinuity between 0o and 360o. The camera film analemma on the right is normal. Here is the 20o latitude image from Fig 3.12 showing this over-the-top behavior, Whenever the "up" z' axis in Frame S' hits the interior of the analemma , this azimuthal discontinuity will be encountered as the year progresses. Red Earth Analemmas at the Latitude 5 o (actual) Now the angle-space analemma rips apart at a lower poi nt, things are more dramatic, but the camera film analemma on the right is just fine. We pass through the azimuthal discontinuity three times! latitude = 20o 12 noon actual Fig 11.29 Moving to 2 PM, the analemma lays over to the right, Section 11: The Solar Analemma 127 latitude = 20o 2 PM actual Fig 11.30 It might seem strange that this analemma tilts over so far only 2 hours after noon. This is the subject of Appendix D. Fig D.2 says that the blue analemma is in fact tilted 80o to the right from vertical. Section 11: The Solar Analemma 128 (g) Analemma case study from the Crimea (Ukraine) Mr. Vasilij Rum yantsev, who works at the Crimean As trophysical Observatory (CrAO) in the town of Nauchny, has a nice analemma picture up on the web at http://vrum.chat.ru/Photo/Astro/analema.htm which is reproduced on the left below (1998-9). His removable but rigidly mounted camera took a shot about every 10 days. Although the strobe sun photos were taken in the early morning, the background picture was taken in the evening. His web page provi des an angle-space representation of the analemma, shown on the right below. Probably this was a pl anning picture he obtained from his own analemma- generating program. Fig 11.31 His web page provides enough information for us to test out our Maple analemma generating code against his, and to then see how his actual photo compares to our screen space prediction. Note: As of Dec 2012 the Ukraine seems intent on "ter minating" the independent Crimean Astrophysical Observatory, founded 1945, proba bly for budgetary reasons, see http://www.crao.crimea.ua/ . The photos were taken close to the observatory's solar telescope which is close to the following coordinates, Section 11: The Solar Analemma 129 Therefore we use, φ1 = longitude = 34.000 + 0/60 + 57.1/3600 = 34.016o θLAT = latitude = 44.000 + 43/60 + 36/3600 = 44.727o . The analemma was taken Aug 1998 to Aug 1999, so we can use the March 1999 equinox data (source given in Section 7 (b)) so we take the GMT of March equinox as 01:41. Mr. Rumyantsev indicates his pictures were taken every 10 days (more or less) at exactly GMT = 05:45. From this we compute his LMT from (8.8), LMT(t) = [ GMT(t) + ( φ 1/15) ] mod 24 . = [05:45 + 34.016/15] mod 24 = 5.75+2.26773 = 8.01773 ≈ 8:01 AM . Finally, the orbital parameters in 1999 were (source given in Section 1 (a)) so θtmax = 23.4399, ε = .016704 and φ per = 12.878o. Entering this data into the Maple program, we ge t the red analemma shown in Fig 11.10 above. If we superpose that red analemma on his data and car efully align the axes, here is the result : Section 11: The Solar Analemma 130 Fig 11.32 The break in our red curve (above 3-16) shows the March equinox which is in agreement with his data. This seems to show that our angle-space Maple progr am agrees with whatever program he used to generate his planning analemma. We now superpose our blue screen-space analemma plot of Fig 11.11 on Mr. Rumyantsev's photo without altering the aspect ratio of either image. We had to set our pan adjust to -6 o to get things to line up (cross hairs are 6o to the left of the fictitious sun location). Section 11: The Solar Analemma 131 Fig 11.33 The fit is reasonable but unfortunately not perfect in the smaller summer lobe. As shown in the previous figure, our computed angle-space an alemma agrees exactly with that of Mr. Rumyantsev, so the photo mismatch is either due an error in our screen space calculation, or to a problem with the camera mount, the former seeming more likely since the winter lobe matches so well. The first known analemma photograph was taken by Dennis di Cicco in 1978-9, Ref [7] . Since that time, perhaps 10 people have put up analemma photos on the web. A notable effort is that of Greek photographer Anthony Ayiomamitis who uses multiple cameras to do simultaneous analemma photos during the same year. He does this with cameras mounted in his yard near Athens, each one strobing at a different hour of the day (circa 2002-2004). To add interest, he used Photoshop to composite in photos of various Greek ruin s. We did Maple fits to several of his photos with results similar to that shown above. The first li nk below shows a collection of Mr. Ayiomamitis's analemma photos which demonstrates the effect of time of day, in agreement with our analemma gallery comments above. The second link gives details of his efforts (five web pages) : http://www.perseus.gr/Astro-Solar-Analemma.htm http://www.astrosurf.com/luxorion/analemma.htm Appendix A: The arctangent function 132 Appendix A: A Note on Various Arctangent Functions When band ying about equations, as we have done freely in this document, one can be cavalier about the meaning of functions like tan-1. But when it comes time to actually compute something, or to make a proper plot of some function, subtle problems arise which can be very annoying. This appendix describes issues involving the arctangent function in various fo rms. These issues are of course present not only in Maple implementations, but in any calculational method of doing something at the nuts and bolts level. (a) arctan2Pi(y,x) and comparison with arctan(y/x) and arctan(y,x) Consider these two equations, x = cosφ y = sin φ where φ runs through the range (0,2 π) so that r = (x,y) moves in a circle. Given x and y, how do we find φ? Consider this picture, Fig A.1 From the signs of x and y, we can certainly determ ine the quadrant of (x,y). So here is how we find φ : If (x,y) is in Quadrant I, then φ = α where α ≡ tan-1(|y|/|x|) If (x,y) is in Quadrant II, then φ = π – α If (x,y) is in Quadrant III, then φ = π + α If (x,y) is in Quadrant IV, then φ = 2π - α . ( A . 1 ) This algorithm is implemented in the following Maple function arctan2Pi(y,x), where we put the arguments in order y,x since it reminds us of tan-1(y/x) : Appendix A: The arctangent function 133 Here are two 3D plots of arctan2Pi(y, x), showing its scale invariant shape, F i g A . 2 As one takes a spiral hike around the qu adrants going I, II, III, IV the value φ = arctan2Pi(y,x) increases smoothly. When using this functi on, we stick to the range 0 ≤ φ < 2π to avoid the discontinuity. NOTE : All equations below are followed by an implicit "mod 2 π". When checking an equation of the form a = b in Maple, one must show that mod(a-b,2 π) = 0. In particular notice that π = -π in this world since mod( π - (-π ),2π) = mod(2 π,2π) = 0. In general, the function arctan2Pi( y,x) has these properties, as one can show by exhausting all the cases in Fig A.1, Appendix A: The arctangent function 134 arctan2Pi (-y,x) = 2 π – arctan2Pi (y,x) // negate y arctan2Pi (y,-x) = π – arctan2Pi (y,x) // negate x arctan2Pi (-y,-x) = π + arctan2Pi (y,x) // negate x and y . (A.2) These rules differ considerably from the principle branch rule that tan-1(-z) = - tan-1(z) . To verify (A.2) we run this Maple code [ modf(a,b) computes a mod b, see Sec 11 (e) ] and all three plots come out looking like this (two views of the same 3D plot) Fig A.3 Rather than write arctan2Pi everywhe re, we have adopted this notation, tan-1( [y]/[x] ) = tan-1( [y] [x] ) ≡ arctan2Pi(y,x) . (A.3) One must not "move signs around" carelessly, thinking of tan-1(-z) = - tan-1(z) . In fact, from (A.2), tan -1( [–y]/[x] ) = 2 π – tan-1( [y]/[x] ) tan-1( [y]/[–x] ) = π – tan-1( [y]/[x] ) tan-1( [–y]/[-x] ) = π + tan-1( [y]/[x] ) . (A.4) Here are some examples of (A.4) : Appendix A: The arctangent function 135 tan-1( [sinφ]/[cosφ] ) = φ tan-1( [-sinφ]/[cosφ] ) = 2π – φ tan-1( [sinφ]/[-cosφ] ) = π – φ tan-1( [-sinφ]/[-cosφ] ) = π + φ . ( A . 5 ) Sometimes we use this abbreviated notation tan -1(a tanβ) ≡ tan-1( [a sinβ]/[cosβ] ) if a ≥ 0 ( A . 6 ) as in (6.12) and (6.13). Of course if a = 1, then tan -1(tanβ) = β . Similarly, tan -1( [c]/[ab] ) = tan-1( [a-1c]/[b] ) if a > 0 (A.7) with corollary tan-1( [ac]/[ab] ) = tan-1( [c]/[b] ) if a > 0 . (A.8) The important fact in such identities is that the quadrant does not change. The cotangent can be handled similarly. We first define cot -1([x]/[y]) ≡ tan-1([y]/[x]) = arctan2Pi(y,x) . (A.9) We claim that, just as for the regular tan-1 and cot-1 functions, cot-1([x]/[y]) = π /2 – tan-1([x]/[y]) . (A.10) This can be verified for (x,y) in any quadrant by showing that the following is true, arctan2Pi(y,x) = π/2 – arctan2Pi(x,y) . Testing this claim, we find that the code produces the same null plot shown in Fig A.3. Combining (A.10) with (A.9) then gives tan-1([y]/[x]) = π /2 – tan-1( [ x ] / [ y ] ) . ( A . 1 1 ) From (A.11) and (A.4) one may show that π/2 – tan -1([y]/[x]) = π + tan-1([-x]/[-y]) (A.12) which we once again verify in Maple. We find that this code, Appendix A: The arctangent function 136 produces the usual null plot of Fig A.3. In the special case (x,y) = (1,y), we know that the normal tan-1(y/x) = tan-1(y) lies in -π /2 to π/2. For our situation this is Quadrant I or Quadrant IV and we then obtain ( don't forget the implicit mod 2 π) tan -1( [y]/[1]) = tan-1( y ) ( A . 1 3 ) then using the second equation of (A.4) we get tan -1( [y]/[-1]) = π – tan-1( y ) . ( A . 1 4 ) Other arctan functions in Maple Maple contains a built-in function arctan(y,x) with tw o arguments. Here is a plot of this function where the axes are the same as in the left side of Fig A.2 for arctan2Pi(y,x) : Fig A.4 The output range of this function is (- π,π) rather than (0,2 π), which is not what we want for our purposes. It turns out that we can simulate arctan2Pi(y,x) by arctan(-y,-x) + π, Appendix A: The arctangent function 137 Fig A.5 but it seemed better to write our own clearly define d function arctan2Pi to be sure what it does. The following three plots using the three arctangent functions demonstrate our motivation in using the arctan2Pi function for our (0,2π ) region of interest: F i g A . 6 (b) arctan2Pis(y,x) First, we de monstrate what this "s" version of arctan2Pi does: Fig A.7 Appendix A: The arctangent function 138 As a (x,y) vector winds around the origin multiple times, the arctangent increases smoothly so that the result tan-1(tanz) = z is realized for z over a larger range than 2π . Here is the inelegant code for this function : The function remembers the x,y quadrant of the previous call to the function and compares that with the x,y quadrant of the current call. If the 4 →1 boundary is crossed, a winding number "wind" is upticked by 1, and if it crosses the other way it is downticked. Then 2 π*wind is added to the arctan2Pi(y,x) result. This function must be initialized before it is used to evaluate a function for some sweep of a variable, As an illustration of the motivation for using this arct an2Pis function, consider the seemingly simple plot of four curves (as in Fig 6.3) which lie very clo se together. Here is the code using arctan2Pis : Appendix A: The arctangent function 139 Fig A.8 When these curves are plotted with the other three arctangent functi ons, we get these annoying results, ( the cases are arctan2Pi(y,x), then arctan(y,x), then arctan(y/x) ) F i g A . 9 F i g A . 1 0 Appendix A: The arctangent function 140 F i g A . 1 The following drawing shows how our four different arctangent functions use different branches of the infinitely-branched general arctangent function: F i g A . 1 2 The normal function φ = tan-1(z) uses only the principle branch and returns φ in -π/2 to π /2. Technical Code Note: The code which generates Fig 6.4 is stat ed above that figure, and the functions like f_red are defined in terms of the arctan2Pis functio n. Since this function is designed only for use in a clean "run" from left to right in a plot of a si ngle function, where it internally manages its winding number "wind", it seems odd that there are no "spikes" in the plots like Fig 6.4. There, the plot algorithm evaluates d_red = f_red - f_green which means arctan2Pis is applied first to f_red and then to f_green, and then after the next plot increment arctan2Pis is again applied to f_red and then to f_green. So the use of arctan2Pis alternates between the red and green functions, viola ting the coder's intent. It nevertheless turns out that things really do work even though the red and green func tions cross the axis near t = 0.2 (76 days) at different times, and this is re gardless of how fine the plotting steps are. Appendix B: Active and Passive Rotations 141 Appendix B: The Active and Passive Views of Vector Rotation (a) Active view Let Frame S have orthon ormal basis vectors en for n = 1,2,3. Let V be some vector in Frame S which has components V n. Here is the expansion of vector V on these basis vectors, V = ΣnVnen where V n = V • en . (B.1) Suppose we create a new vector in Frame S, called V ', by applying some rotation R to vector V, V' = R V meaning V' i = ΣjRijVj . (B.2) We refer to this as " the active rotation of vector V into vector V', all in Frame S ". We can expand this new vector V' on the same basis vectors en as follows, V' = ΣnV'nen where V' n = V' • en ( B . 3 ) where the numbers V' n are the components of vector V' in Frame S. The interpretation of the numbers V' n as just described is "the active point of view of a rotation". Here is a simple example of such an active rotation where R = R z(ψ). Fig B.1 (b) Passive view There is another way to interpret the numbers V' n. Suppose we create a new reference Frame S' which has unit basis vectors e'n which are related to those of Frame S as follows: e'n = R-1en n = 1 , 2 , 3 ( B . 4 ) That is to say, these vectors are rotated "backwards" compared to what is shown in (B.2). The new interpretation of the numbers V' n is (we claim) that these numbers are the components of the vector V when it is observed from Frame S'. That is to say, we claim that the expansion of vector V in Frame S' looks like this, where V n' are as given in (B.3). Appendix B: Active and Passive Rotations 142 V = ΣnV'ne'n where V' n = V • e'n . ( B . 5 ) Dotting both sides of V = ΣnV'ne'n into e'm shows that V' n = V • e'm (since e'n• e'm = δn,m), so basically we want to show that V' n = V • e'n of (B.5) is the same as V' n = V' • en of (B.3) . Proof : V' n = V' • en // from (B.3) = [ R V] • [Re' n] // from (B.2), and (B.4) inverted = [ R V ]T [Re'n] // the above just written in matrix notation, T = transpose vector = VTRTRe'n // rule for applying T to product of matrices/vectors = V T e'n // RTR = 1 because a rotation is a real orthogonal matrix = V • e'n // going back to dot product notation = V ' n // from (B.5) Thus the numbers V' n appearing in (B.3) and (B.5) are in fact the same numbers so our claim is verified. Here is a simple example of such a p assive rotation. Notice that the vector V stays put, and the unit vectors are rotated backwards relative to the ac tive rotation shown in the previous picture, F i g B . 2 (c) Summary The num bers V'n can be interpreted in two ways, active and passive: (1) The V' n are the components of a vector V' in Frame S (with en) obtained from V' = R V . (2) The V' n are components of the vector V in Frame S', where e'n = R-1en for Frame S'. Suppose we know that there is some matrix R -1 such that the basis vectors of Frame S' are related to those of Frame S in this way Appendix B: Active and Passive Rotations 143 e'n = R-1en . If we want to compute the components of a vector V in Frame S', we may use this formula, V' = R V . The simple rule is that "vectors rotate backwards from the unit vectors". This rule is invoked just after equation (6.1). Appendix C: Shape of the e=0 analemma 144 Appendix C: Shape of the north pole analemma when e = 0 (a) Parametric equation for the analemma We assu me that the orbit of the earth is circular (e = 0, "blue earth" of Section 6 (d) ). When e = 0, we have the following set of facts: φ' = Ω 0τ // since circular orbit, φ' as in Fig 3.4 measured from March equinox ψ = Ω0t // since circular orbit, ψ as in Fig 3.4 measured from perihelion τ = t - tMe // valid for any e (7.2) tMe = α/Ω0 // α ≡ (π/2-φper) as in (1.9) (4.12) k = [Ω0 tMe - α] = 0 (7.17) φ' = ψ - α // valid for any e (1.10) From the first and last lines we get this e = 0 fact, ψ - α = Ω0τ . ( C . 1 ) The units are taken to be hours for t and τ and t Me and radians per hour for Ω 0. Now, at the north pole the continuous functions one samples to get the two ax es of the analemma have this simple form, as was shown in (11.24), vertical = θ' EL= θt = sin-1[sinθtmax sin(Ω0τ)] // (11.24), (11.1) horizontal = Φ' = π + ωτ + φb0 - Δφ // (11.24) Δφ = – Ω0τ – k + tan-1[cosθtmax tan(Ω0τ)] // (11.3) (C.2) Using s = ψ-α = Ω0τ from (6.17), one can write Δφ in the alternate form given in (6.22), Δφ = tan-1 [(cosθtmax – 1)tan(Ω0τ) 1 + cosθtmaxtan2(Ω0τ) ] . ( C . 3 ) Strobing Φ' at LMT on day N gives Φ'N = π + ωτN + φb0 - ΔφN τN = [τ0 + 2 4 N ] ; (11.9) τ0 = (LMT – 12) – ( φb0/ω) . // since k = 0, ω = π/12 rad/hr (11.9) Since ω = π/12, and since Φ' N is really mod 2 π, we can replace ωτN → ωτ0 in Φ 'N. That is to say, ωτN = ω[τ0 + 24N] = ωτ0 + (π/12)24N = ωτ0 + 2πN . (C.4) Then Appendix C: Shape of the e=0 analemma 145 Φ'N = π + ωτN + φb0 - ΔφN = π + ωτ0 + φb0 - ΔφN = π + ω[(LMT – 12) – ( φb0/ω)] + φb0 - ΔφN = π + ω(LMT – 12) - ΔφN . ( C . 5 ) This result can also be obtained directly from the Φ'N equation in (11.9) by setting θ1 = 0. Strobing θt in (C.2) and Δφ in (C.3) gives, along with (C.5), θ t,N = sin-1[sinθtmax sin(Ω0τN)] Φ'N = π + ω(LMT – 12) - ΔφN ΔφN = tan-1 [(cosθtmax – 1)tan(Ω0τN) 1 + cosθtmaxtan2(Ω0τN) ] τN = [τ0 + 24N]; τ0 = [ LMT – 12 – ( φb0/ω) ] . ( C . 6 ) To reduce clutter, we now define these symbols, c ≡ cosθ tmax s ≡ sinθtmax β ≡ Ω0τN x ≡ Φ'N - π - ω(LMT – 12) = - ΔφN y ≡ θt,N , ( C . 7 ) so the first three equations in (C.6) may then be written y = sin -1[s sin β] x = - ΔφN ΔφN = tan-1 [(c – 1)tan( β) 1 + ctan2(β) ] , ( C . 8 ) which we restate one more time as x = - tan-1 [(c – 1)tan( β) 1 + ctan2(β) ] y = sin-1[s sin β ] . ( C . 9 ) As the parameter β ranges from 0 to 2 π ( treating β now as a continuous parameter), (x,y) traces out the continuous blue-earth analemma in angle space. Just to make sure we are on track, we can plot this locus (with unequal axis scales each in radians), Appendix C: Shape of the e=0 analemma 146 Fig C.1 Plotting for a θtmax = 10,20...90 shows that the analemmas get larger as θtmax increases, and also change shape. The triangular analemma is for θtmax = 90. 1 0o through 400 300 through 900 Fig C.2 (b) Locus of the analemma in angle space: symmetry The curves above certainly appear to b e symmetric in both directions. We can rewrite (C.9) as tanx = - [(c – 1)tan( β) 1 + ctan2(β) ] siny = s sin β . ( C . 1 0 ) Appendix C: Shape of the e=0 analemma 147 From the second of (C.10) one ge ts (the goal is to eliminate β), sin2β = sin2y/s2 cos2β = 1 – sin2y/s2 => tan2β = sin2y/s2 / [1 – sin2y/s2] = sin2y / [s2 – sin2y] . (C.11) The first equation of (C.10) can be squared to yield tan 2x = (c-1)2tan2β / (1 + c tan2β)2 . ( C . 1 2 ) Using (C.11) in (C.12), tan 2x = { (c-1)2(sin2y / [s2 – sin2y]) } / {1 + c(sin2y / [s2 – sin2y]) }2 Multiply top and bottom by [s2 – sin2y]2 to get tan2x = { (c-1)2sin2y [s2 – sin2y]} / {[s2 – sin2y] + c(sin2y ) }2 or tan2x = { (c-1)2sin2y [s2 – sin2y]} / [s2 – sin2y + csin2y ]2 or tan 2x = (c-1)2sin2y (s2 – sin2y) / [s2 + (c-1)sin2y ]2 . (C.13) Thus β has been eliminated and we have obtained the lo cus of the analemma. We can verify (C.13) by doing a scanning "implicit plot" in Maple as follows: Fig C.3 The locus (C.13) is extremely non-lin ear in (x,y) space, but does expose the symmetry we have been seeking: if x → –x or y →–y (or both), the locus is unchanged. Appendix C: Shape of the e=0 analemma 148 (c) How β va ries along the analemma Recall the earlier parametric equation for the analemma, x(β) = - tan-1 [(c – 1)tan( β) 1 + ctan2(β) ] y(β) = sin-1[s sin β] . (C.9) Using these facts tan(β+π) = tanβ tan(2 π-β) = -tanβ tan( π-β) = -tan(β) sin(β+π) = -sinβ sin(2 π-β) = -sinβ sin( π-β) = sin(β) one finds x(β+π) = x(β ) x(2 π-β) = -x(β) x( π-β) = -x(β) y(β+π) = -y(β) y(2 π-β) = -y(β) y(π -β) = y(β) or x,y → x,-y x,y → -x,-y x,y → -x,y mirror in x=0 invert mirror in y=0 Here then is the situation with parameter β in (0,2π) for a typical e=0 analemma F i g C . 4 Note that continuous β = Ω0τ = φ', the azimuth angle of Fig 1.3 and Fig 3.4. Thus, β=φ' 0 March equinox π/2 June solstice π September equinox 3π/2 December solstice 2π March equinox Appendix C: Shape of the e=0 analemma 149 (d) Limit of analemma shape for small θtmax : Gerono! We have seen in Fig C.2 how the analemma approaches a pair of triangles in the limit θtmax → 900. In the opposite limit that θtmax → 0, there is also an interesting limiting shape. Recall again, x = - tan-1 [(c – 1)tan( β) 1 + ctan2(β) ] c = cos θtmax y = sin-1[s sin β ] s = sin θtmax . (C.9) When θtmax is small, s is small, and so y is small. Since c → 1, x is also small. Looking at (C.13) tan 2x = (c-1)2sin2y (s2 – sin2y) / [s2 + (c-1)sin2y ]2 (C.13) we then have, for very small θ tmax, x2 ≈ (c-1)2y2 (s2 – y2) / [s2 + (c-1)y2 ]2 . Now insert c ≈ 1-s 2/2 so that (c-1) ≈ - s2/2 to get x 2 ≈ ( s4/4)y2 (s2 – y2) / [s2 - (s2/2)y2 ]2 x 2 ≈ ( 1/4)y2 (s2 – y2) / [1 - (1/2)y2 ]2 x 2 ≈ y2 (s2 – y2) / [2 - y2 ]2 x 2 ≈ y2 (s2 – y2)/4 4x2 = s2y2- y4 . // in limit θtmax → 0 Rescale to new variables Y = (1/s) y X = (2/s 2)x and the small- θ tmax analemma becomes Y4 – Y2 + X2 = 0 . This closed quartic curve is known as the lemniscate of Gerono (not of Bernoulli), Appendix C: Shape of the e=0 analemma 150 Fig C.5 Wiki shows this curve with the axes swapped, x 4 – y2 + x2 = 0 http://en.wikipedia.org/wiki/Lemniscate_of_Gerono F i g C . 6 Appendix D: Analemma tilt on the film 151 Appendix D: Calculation of how much the analemma is tipped on film (a) Calculation of the tip angle from vertical We shall co mpute this analemma tip angle as a function of colatitude θ1 of the observer site, and LMT, the local mean time of the strobe photos. The plan here is to create a vector Δr' which represents the long direction of an analemma. We make this arrow point from the December solstice end of th e analemma to the June solstice end. The arrow is then shrunk in length and mapped in to screen space, where it becomes ( Δxs,Δys). We then compute the tip angle of this vector on the film away from vertical as follows tip ≡ angle from vertical = tan -1(Δxs/Δys) . (D.1) For an arrow with Δx s = 0, tip = 0o, a vertical analemma. For an arrow with Δys = Δxs, tip = tan-1(1) = 45o, tipped to the right. For an arrow with Δys = - Δxs, tan-1(-1) = -45o, tipped to the left. To simplify matters we shall use the blue earth analemma which has e = 0. As shown in Appendix C (b), this analemma is symmetric both vertically and horizontally. The vertical symmetry implies that Δφblue = 0 at both solstices, but here is a direct calculation of this claim. We first assemble a few facts, φ' = Ω0τ // since circular orbit, φ' as in Fig 3.4 measured from March equinox ψ = Ω0t // since circular orbit, ψ as in Fig 3.4 measured from perihelion τ = t - tMe (7.2) tMe = α/Ω0 . (4.12) Then from (6.22), Δφ blue(t) = tan-1 [(cosθtmax – 1)tan(s) 1 + cosθtmaxtan2(s) ] = tan-1 [(cosθtmax – 1)sin(s)cos(s) cos2(s) + cosθtmaxsin2(s) ] (6.22) where s = ψ-α since q = 0, see (6.17). But from the facts collected above, s = ψ-α = Ω0t- α = Ω0(τ + tMe) - α = Ω0τ + Ω0 tMe - α = Ω0τ + α - α = Ω0τ = φ' . At the two solstices, as show n in Fig 3.4 we have s = φ' = π/2 and 3π/2. For either value cos(s) = 0 and therefore Δφblue(t) = 0. In this case we conclude from (11.8) that, (φ b-φ)N = (ω[LMT- 12]) . N = day of either solstice strobe time, since ΔφN = 0 (11.8) Recalling then (7.14) which gives the Frame S' Cartesian coordinates of the sun, Appendix D: Analemma tilt on the film 152 x' = R cos θt sin(φ –φb) y' = R [sin θt sinθ1 – cosθt cosθ1 cos(φ –φb) ] z' = R [sin θt cosθ1 + cosθt sinθ1 cos(φ-φb) ] , (7.14) we can install (11.8) in three places and then evaluate the result at the two solstices, θt = ± θtmax, x' = – R cos θ tmax sin (ω[LMT- 12]) y' = R [sin θtmax sinθ1 – cosθtmax cosθ1 cos(ω [LMT- 12]) ] z' = R [sin θtmax cosθ1 + cosθtmax sinθ1 cos(ω[LMT- 12])] // June solstice (D.2) x' = – R cos θ tmax sin (ω[LMT- 12]) y' = R [-sin θtmax sinθ1 – cosθtmax cosθ1 cos(ω [LMT- 12]) ] z' = R [-sin θtmax cosθ1 + cosθtmax sinθ1 cos(ω[LMT- 12])] . // Dec solstice (D.3) Our difference arrow of interest is then Δx' = 0 Δy' = 2 R sin θ tmax sinθ1 Δz' = 2 R sinθ tmax cosθ1 ( D . 4 ) which is independent of LMT. This may seem stra nge, but a little pondering with a globe prop and one's right hand to model the Frame S' axes shows it to be true. Recall next the eye transformation to camera coordinates of (11.20), x c = -cosαpx' + sinαpy' yc = sinαt sinαp x' + sinα t cosαp y' + cosαt z' zc = -cosαt sinαp x' - cosαt cosαp y' + sinαt z ' . (11.20) Applying this to the vector Δr' of (D.4) gives (matrix transformations are linear so this is allowed) Δxc = -cosαpΔx' + sinαpΔy' Δyc = sinαt sinαp Δx' + sinαt cosαp Δy' + cosαt Δz' Δzc = -cosαt sinαp Δx' - cosαt cosαp Δy' + sinαt Δz ' (11.20) or Δx c = sinαp2 R sinθtmax sinθ1 Δyc = sinαt cosαp 2 R sinθtmax sinθ1+ cosαt 2 R sinθtmax cosθ1 Δzc = - cosαt cosαp 2 R sinθtmax sinθ1+ sinαt 2 R sinθtmax cosθ1 or Δx c = (2Rsin θtmax) sinαpsinθ1 Δyc = (2Rsin θtmax)[ sinαt cosαp sinθ1+ cosαt cosθ1 ] Δzc = (2Rsin θtmax)[- cosαt cosαp sinθ1+ sinαt cosθ1 ] . (D.5) Now we scale down our arrow so it is very short, and we imagine that the tail of the arrow is now sitting at the central point of the blue analemma (the fictiti ous sun location), still pointin g to the June solstice end Appendix D: Analemma tilt on the film 153 of the analemma. That central point has x c = yc = 0 since the camera is pointing right at this point. So we think of this small (differential) arrow as Δxc = ε sinαpsinθ1 Δyc = ε [ sinαt cosαp sinθ1+ cosαt cosθ1 ] Δzc = ε [- cosαt cosαp sinθ1+ sinαt cosθ1 ] ε << 1 . (D.6) The reason for this scaling down is the next processing step, the perspective transformation, x s = zo(xc/zc) ys = zo(yc/zc) , (11.22) which is a nonlinear transformation, so we can't just claim Δ x s = zo(Δxc/Δzc). Instead we differentiate (11.22) to get Δxs = zo(zcΔxc- xcΔzc)/zc2 = zo(zcΔxc)/zc2 Δys = zo(zcΔyc- ycΔzc)/zc2 = zo(zcΔyc)/zc2 ( D . 7 ) where we have used the fact that x c = yc = 0 for the tail end of our differential arrow. Then the tip angle is given by tip ≡ tan-1(Δxs/Δys) = tan-1(Δxc/Δyc) = tan -1( [sinαpsinθ1] / [ sinαt cosαp sinθ1+ cosαt cosθ1 ]) . (D.8) [ It happens that one gets this same result assuming Δx s = zo(Δxc/Δzc) which was dismissed above.] From (11.10) and (11.23) we find that α t = π/2 - θ 'f θ'f = cos-1[sinθ1 cos (ω[LMT- 12] )] αp = Φ'f - π Φ'f = π + tan-1 ([sin (ω[LMT- 12])] / [cos θ1 cos (ω [LMT- 12]) ]) which implies that α t = sin-1[sinθ1 cos (ω[LMT- 12] )] αp = tan-1 ([sin (ω[LMT- 12])] / [cos θ1 cos (ω[LMT- 12]) ]) . (D.9) The final result for the analemma tip angle is then this : tip(θ1, LMT) = tan-1( [sinαpsinθ1] / [ sinαt cosαp sinθ1+ cosαt cosθ1 ]) αt = sin-1[sinθ1 cos (ω[LMT- 12] )] αp = tan-1 ([sin (ω[LMT- 12])] / [cos θ1 cos (ω[LMT- 12]) ]) . (D.10) The tip angle is a function of θ 1 and LMT, as one would expect. Appendix D: Analemma tilt on the film 154 We duly enter equations (D.10) into Maple, and make a 3D plot of tip( θ1, LMT) Fig D.1 What the plot shows is that near the equator θ 1 = π/2 ≈ 1.57, the analemma is vertical at LMT = 12, but as soon as one moves to 12.01 the analemma instantly tips over to its right, tip = π/2, and at 11.99 it tips to its left, tip = - π/2. Moving up a bit from the equator, say to latitude 5o, the tip is still fast, but not instantaneous. Probably this situation is clearer in the following plot, Appendix D: Analemma tilt on the film 155 Fig D.2 Again, at the equator the analemma tips over instantly as one moves away from noon. At 2o latitude the tipping over is fairly quick, being perhaps 75o laid over to the right for an LMT = 12:30 PM analemma. We confirm the graph with direct calculation, On the other hand, at 40 o latitude, the analemma tips 45o to the right at about 4 PM, For the southern hemisphere, just nega te the latitude labels in Fig D.2. (b) Physical explanation of the tip angle behavior An explanation of the fast tip-over near the equator is as follows. Imagine sitting in a swivel chair looking due south from a site at the equator, camera in hands . At noon the fictitious sun is exactly overhead, so one must tilt the camera up 90o to make it point at the fictitious sun which marks the center of the analemma. There is no need for panning, the analemma appears vertical on the film, since it runs north- south on the celestial sphere (see Fig 11.3). At LM T = 12.001 PM, the fictitious sun has moved a little to the west, so one must swivel that chair 90o to the west in order to point the camera at the new fictitious Appendix D: Analemma tilt on the film 156 sun position. But then the analemma is horizontal on th e film since the viewer is pointing up and west but the analemma is still north-south! It "tipped over" instantaneously on the film as we moved away from noon. Similarly at LMT = 11.999 one must instead swivel the chair 90o to the east, so the analemma is horizontal the other way. In our film projection method, we are only allowed to tilt and then "pan" about the vertical axis to match the polar angles of the fictitious sun; we are not allowed to "roll" the camera about its z c axis, nor are we allowed to truly pan it about the y c camera axis. Now we move say to 2o north latitude. Some fast chair swive ling is still required to track the sun through noon with the camera, but the swiveling is not quite as violent. At 12:30 PM the fictitious sun has moved 7.5o west (15o per hour), and it is 2o to the south. Here is a picture F i g D . 3 The swivel angle is then tan-1(7.5/2) = 75o , which is a lot of swivel (though not 90o) so the analemma still being north-south on the celestial s phere is tipped to the right by 75o. It only took a half hour to get it tipped down this great amount. This is in ag reement with the above graph and calculation. One can simulate this perspective effect by making a rectangular viewport with two hands and pretending it is the camera. Find a room with a swivel chair which has parallel lines of some sort running north-south on the ceiling. A directly overhead line is the noon analemma, while a parallel line to the west is the analemma slightly after noon. References 157 References Various web links appear in-line in this document a nd most are not repeated here. Links for web pages indicated in-line or below may become invalid, but a search on title or keywords can reveal their new homes. [1] H. Goldstein, Classical Mechanics (Addison-Wesley, Boston, 1950). [2] I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, 7th Ed. (Academic Press/ Elsevier, Amsterdam, 2007). [3] N. Hungerbühler , "Another Design for Sundials" (9p), http://homeweb2.unifr.ch/hunge rbu/pub/sonnenuhr/s undial.ps [4] Our Fig 6.8 plot of the Equation of Time appears here in the Wikipedia Commons, http://en.wikipedia.org/wiki/File:Zeitgleichung.png . The web page also contains the actual code used to generate the plot. It gives two references, one of which is item [3] above. The other reference of Blatte r is a dead link, but the "approximation by Blatter" referred to in the code is our equation (4.32). The code uses φ per = 11.50 (α = 78.5o) and θtmax = 23.45o which were correct for year 1919. [5] F.W.J. Olver et. al ., NIST Handbook of Mathematical Functions (National Institute of Standards and Technology, Cambridge University Press, Cambridge , 2010). This excellent reference book is the follow- on to the famous 1964 “Abramowitz and Stegun” handbook. [6] P. Lucht, Rotating Frames of Reference (http://user.xmission.com/~rimrock/ , 2012) . This paper uses names Frame S and S', basis vectors en and e'n, origin-connecting vector b, and other notations which have been utilized in the current paper. Its A ppendix A has information on spherical coordinates. [7] D. di Cicco, "Exposing the Analemma" (Sky and Telescope, June 1979, pp536-540). D. di Cicco, "Photographing the Analemma" (Sky and Telescope, March 2000, pp135-140). Dennis di Cicco apparently took the first-ever so lar analemma strobe photo. His imaginary club of analemma photographers is slowly increasing its membership. [8] V. Rumyantsev, http://vrum.chat.ru/Photo/Astro/analema.htm , analemma photo and details. [9] A. Ayiomamitis, "The Solar Analemma", http://www.astrosurf.com/luxorion/analemma.htm . This is a long description of the analemma adventures of Mr. Ayiomamitis. His photos are displayed on another site, "Solar Image Gallery - Analemma", http://www.perseus.gr/Astro-Solar-Analemma.htm . [10] NASA, "Determination of Earth's Orbital Para meters", as used in Section 1 and elsewhere. http://aom.giss.nasa.gov/srorbpar.html . One can obtain actual numbers for e, θtmax, and φper. References 158 [11] NASA, "Time and Date of Vernal Equinox", gives times of solstices, equinoxes, perihelion and aphelion for the Earth, http://aom.giss.nasa.gov/srver4x3.html .