1====================
2The PyEphem Tutorial
3====================
4
5Version 2007.November.1, for PyEphem version 3.7.2.2 and later.
6
7The `PyEphem library`_ provides Python programmers
8with access to the astronomical routines
9used by the `XEphem`_ interactive astronomical ephemeris application;
10its author, `Elwood Charles Downey`_, has generously granted permission
11for PyEphem to be built upon his work.
12
13After installing the module,
14you can use it in Python with the statement:
15
16.. _PyEphem library: http://rhodesmill.org/pyephem/
17.. _XEphem: http://www.clearskyinstitute.com/xephem/
18.. _Elwood Charles Downey: http://www.clearskyinstitute.com/resumes/ecdowney/resume.html
19
20>>> import ephem
21
22This tutorial assumes that you are familiar with astronomy,
23but necessarily all of the issues surrounding astronomical calculation;
24those who find its discussions tedious
25will probably just want to read over its examples
26to quickly familiarize themselves with how PyEphem works.
27
28First Steps
29-----------
30
31PyEphem will compute the positions of celestial bodies on particular dates,
32and can determine where in the sky they appear from any location on earth.
33
34When using PyEphem,
35you will usually create instances of the bodies that interest you,
36compute their position for various dates and perhaps geographic locations,
37and finally print or display the result.
38For example,
39to determine the location and brightness of Uranus
40on the night it was discovered
41we simply create a ``Uranus`` object
42and ask where it was on the 13th of March, 1781:
43
44>>> u = ephem.Uranus()
45>>> u.compute('1781/3/13')
46>>> print('%s %s %s' % (u.ra, u.dec, u.mag))
475:35:45.28 23:32:54.1 5.6
48>>> print(ephem.constellation(u))
49('Tau', 'Taurus')
50
51Calling ``compute()`` sets many attributes of a body,
52beyond the right ascension, declination, and magnitude printed here;
53see the :doc:`quick`
54for other attributes that ``compute()`` sets.
55You see that measurements are formatted as an astronomer would expect:
56dates are expressed as year, month, and day, delimited by slashes;
57right ascension as hours of arc around the celestial equator;
58and declination as degrees north of the equator.
59The colons between the components of an angle are a compromise —
60the more traditional 21°46′15.66′′ is not possible
61with the symbols on a standard computer keyboard.
62
63The code above created and used only one instance of ``Uranus``,
64but you can also have several going at once.
65For example,
66to determine how close Neptune and Jupiter lay
67as Galileo famously observed them —
68he was busy watching the Jovian moons he had discovered two years earlier
69and, though Neptune had moved overnight, he dismissed it as a background star
70and left its discovery to wait another two hundred years —
71we create one instance of each planet and compare their positions:
72
73>>> j = ephem.Jupiter('1612/12/28')
74>>> n = ephem.Neptune('1612/12/28')
75>>> print("%s %s %s" % (j.ra, j.dec, j.mag))
7611:48:20.52 2:41:13.6 -1.96
77>>> print("%s %s %s" % (n.ra, n.dec, n.mag))
7811:49:15.77 2:37:04.5 7.92
79>>> print(ephem.separation(j, n))
800:14:24.6
81
82Notice that, while in our first example
83we created our ``Uranus`` instance
84and called its ``compute()`` method as two separate steps,
85we have here taken the shortcut of providing dates
86as we created each planet;
87in general, any arguments you provide when creating a planet
88are used to ``compute()`` its initial position.
89The ``separation()`` function
90computes the angle in degrees between two bodies,
91as measured by their right ascension and declination.
92In this case,
93the separation of 0°14′
94was small enough to place both planets in Galileo's field of view.
95
96You can even create several instances of the same body.
97Let's compare how far Mars moves in one day at perihelion versus aphelion,
98and verify that its speed is greater when closer to the Sun:
99
100>>> def hpos(body): return body.hlon, body.hlat
101>>> ma0 = ephem.Mars('1976/05/21')    # ma: mars near aphelion
102>>> ma1 = ephem.Mars('1976/05/22')
103>>> print(ephem.separation(hpos(ma0), hpos(ma1)))
1040:26:11.4
105>>> mp0 = ephem.Mars('1975/06/13')    # mp: mars near perihelion
106>>> mp1 = ephem.Mars('1975/06/14')
107>>> print(ephem.separation(hpos(mp0), hpos(mp1)))
1080:38:05.2
109
110Here we wanted to measure the motion of Mars around the Sun,
111but ``separation()`` normally compares
112the right ascension and declination of two bodies —
113which would measure the motion of Mars across the sky
114of the moving platform of our earth.
115So instead of giving ``separation()`` the Mars instances themselves,
116we specifically provided
117the heliocentric longitude and latitude of each instance,
118revealing how far Mars moved around the Sun
119regardless of how this motion appeared from earth.
120
121In general ``separation()`` can measure the angle
122between any pair of spherical coordinates,
123so long as the elements of each coordinate are spherical longitude
124(angle around the sphere)
125followed by spherical latitude
126(angle above or below its equator).
127Each pair should be provided as a two-item sequence like a tuple or list.
128Appropriate coordinate pairs include right ascension and declination;
129heliocentric longitude and latitude;
130azimuth and altitude;
131and even the geographic longitude and latitude of two locations on earth.
132
133Computing With Angles
134---------------------
135
136Sometimes you may want to perform computations with times and angles.
137Strings like ``'7:45:45.15'`` are attractive when printed,
138but cumbersome to add and multiply;
139so PyEphem also makes times and angles available as floating-point numbers
140for more convenient use in mathematical formulae.
141
142All angles returned by PyEphem are actually measured in radians.
143Let us return to our first example above,
144and examine the results in more detail:
145
146>>> u = ephem.Uranus('1871/3/13')
147>>> print(str(u.dec))
14822:04:47.4
149>>> print('%.12f' % float(u.dec))
1500.385365877213
151>>> print('%.11f' % (u.dec + 1))
1521.38536587721
153
154The rule is that angles become strings when printed or given to ``str()``,
155but otherwise act like Python floating point numbers.
156The same thing happens when you set an angle:
157a string is interpreted as degrees or hours
158(hours if you are setting right ascension, degrees everywhere else);
159while a float is interpreted as radians.
160
161>>> print(ephem.degrees('90.0'))
16290:00:00.0
163>>> print(ephem.degrees(3.141593))
164180:00:00.1
165
166Note that the format operator ``%`` can return either value,
167depending on whether you use ``%s`` or one of the numeric formats:
168
169>>> print("as a string: %s, as a float: %f" % (u.dec, u.dec))
170as a string: 22:04:47.4, as a float: 0.385366
171
172As an example computation,
173we can verify Kepler's Second Law of planetary motion —
174that a line drawn from a planet to the sun
175will sweep out equal areas over equal periods of time.
176We have already computed two positions for Mars near its aphelion
177that are one day apart
178(and defined a helpful ``hpos()`` function; see above).
179We can estimate the actual distance it moved in space that day
180by multiplying its angular motion in radians by its distance from the Sun:
181
182>>> aph_angle = ephem.separation(hpos(ma0), hpos(ma1))
183>>> aph_distance = aph_angle * ma0.sun_distance
184>>> print('%.13f' % aph_distance)
1850.0126911122281
186
187So, it moved nearly 0.013 AU in a single day (about 1.9 million kilometers).
188A line drawn between it and the sun would have, roughly,
189filled in a triangle whose base is 0.013 AU,
190whose height is the distance to the Sun,
191and whose area is therefore:
192
193>>> aph_area = aph_distance * ma0.sun_distance / 2.
194>>> print('%.13f' % aph_area)
1950.0105710807908
196
197According to Kepler our results should be the same
198for any other one-day period for which we compute this;
199we can try using the two Mars positions from near perihelion:
200
201>>> peri_angle = ephem.separation(hpos(mp0), hpos(mp1))
202>>> peri_distance = peri_angle * mp0.sun_distance
203>>> peri_area = peri_distance * mp0.sun_distance / 2.
204>>> print('%.13f' % peri_area)    # the area, to high precision, is the same!
2050.0105712665517
206
207Despite the fact that Mars moves twenty percent faster at perihelion,
208the area swept out — to quite high precision — is identical,
209just as Kepler predicted.
210Some of the tiny difference between the two numbers we got
211results from our having approximated sectors of its orbit as triangles;
212the rest comes from the pertubations of other planets
213and other small sources of irregularity in its motion.
214
215When you use an angle in mathematical operations,
216Python will return normal floats that lack the special power
217of printing themselves as degrees or hours or arc.
218To turn radian measures back into printable angles,
219PyEphem supplies both a ``degrees()`` and an ``hours()`` function.
220For example:
221
222>>> print('%.13f' % (peri_angle * 2))
2230.0221584026149
224>>> print(ephem.degrees(peri_angle * 2))
2251:16:10.5
226
227You may find that your angle arithmetic often returns angles
228that are less than zero or that exceed twice pi.
229You can access the ``norm`` attribute of an angle
230to force it into this range:
231
232>>> deg = ephem.degrees
233>>> print(deg(deg('270') + deg('180')))
234450:00:00.0
235>>> print(deg(deg('270') + deg('180')).norm)
23690:00:00.0
237
238Computing With Dates
239--------------------
240
241PyEphem only processes and returns dates that are in Universal Time (UT),
242which is simliar to Standard Time in Greenwich, England,
243on the Earth's Prime Meridian.
244If you need to display a PyEphem time in your own timezone,
245use the ``localtime()`` function,
246which returns a Python ``datetime`` object:
247
248>>> d = ephem.Date('1984/12/21 15:00')
249>>> ephem.localtime(d)
250datetime.datetime(1984, 12, 21, 10, 0)
251>>> print(ephem.localtime(d).ctime())
252Fri Dec 21 10:00:00 1984
253
254As you can see from this result,
255I am writing this *Tutorial* in the Eastern Time zone,
256which in the winter is five hours earlier than the time in Greenwich.
257
258PyEphem actually represents dates
259as the number of days since noon on 1899 December 31.
260While you will probably not find
261the absolute value of this number very interesting,
262the fact that it is counted in days
263means you can move one day forward or backward
264by adding or subtracting one.
265The rules described above for angles hold for floats as well:
266you can create them with ``ephem.Date()``,
267but after doing arithmetic on them
268you must pass them back through ``ephem.Date()``
269to turn them back into dates:
270
271>>> d = ephem.Date('1950/2/28')
272>>> print(d + 1)
27318321.5
274>>> print(ephem.Date(d + 1))
2751950/3/1 00:00:00
276
277The ``ephem`` module provides three constants
278``hour``, ``minute``, and ``second``,
279which can be added or subtracted from dates
280to increment or decrement them by the desired amount.
281
282You can specify dates in several formats;
283not only can the strings that specify them
284use either floating point days or provide hours, minutes, and seconds,
285but you can also provide the components of the date in a tuple.
286The following assignments are all equivalent:
287
288>>> d = ephem.Date(34530.34375)
289>>> d = ephem.Date('1994/7/16.84375')
290>>> d = ephem.Date('1994/7/16 20:15')
291>>> d = ephem.Date((1994, 7, 16.84375))
292>>> d = ephem.Date((1994, 7, 16, 20, 15, 0))
293
294And to complement the fact that you can specify dates as a tuple,
295two methods are provided for extracting the date as a tuple:
296``triple()`` returns a year, month, and floating point day,
297while ``tuple()`` provides everything down to floating point seconds.
298After any of the above calls,
299the date can be examined as:
300
301>>> print('as a float: %f\nas a string: "%s"' % (d, d))
302as a float: 34530.343750
303as a string: "1994/7/16 20:15:00"
304>>> print(d.triple())
305(1994, 7, 16.84375)
306>>> print(d.tuple())
307(1994, 7, 16, 20, 15, 0.0)
308
309Any PyEphem function argument that requires an angle or date
310will accept any of the representations shown above;
311so you could, for instance,
312give a three-element tuple
313directly to ``compute()`` for the date,
314rather than having to pass the tuple through the
315``Date()`` function before using it
316(though the latter approach would also work).
317
318Computations for Particular Observers
319-------------------------------------
320
321The examples so far have determined
322the position of bodies against the background of stars,
323and their location in the solar system.
324But to observe a body we need to know more —
325whether it is visible from our latitude,
326when it rises and sets,
327and the height it achieves above our horizon.
328In return for this more detailed information,
329PyEphem quite reasonably demands to know our position on the earth's surface;
330we can provide this through an object called an ``Observer``:
331
332>>> gatech = ephem.Observer()
333>>> gatech.lon, gatech.lat = '-84.39733', '33.775867'
334
335When the ``Observer`` is provided to ``compute()``
336instead of a simple date and epoch,
337PyEphem has enough information
338to determine where in the sky the body appears.
339Fill in the ``date`` and ``epoch`` fields of the ``Observer``
340with the values you would otherwise provide to ``compute()``;
341the epoch defaults to the year 2000 if you do not set it yourself.
342As an example, we can examine the 1984 eclipse of the sun from Atlanta:
343
344>>> gatech.date = '1984/5/30 16:22:56'   # 12:22:56 EDT
345>>> sun, moon = ephem.Sun(), ephem.Moon()
346>>> sun.compute(gatech)
347>>> moon.compute(gatech)
348>>> print("%s %s" % (sun.alt, sun.az))
34970:08:39.2 122:11:26.4
350>>> print("%s %s" % (moon.alt, moon.az))
35170:08:39.5 122:11:26.0
352
353For those unfamiliar with azimuth and altitude:
354they describe position in the sky by measuring angle around the horizon,
355then angle above the horizon.
356To locate the Sun and Moon in this instance,
357you would begin by facing north and then turn right 122°,
358bringing you almost around to the southeast
359(which lies 125° around the sky from north);
360and by looking 70° above that point on the horizon —
361fairly high, given that 90° is directly overhead —
362you would find the Sun and Moon.
363
364Eclipses are classified as *partial*
365when the Moon merely takes a bite out of the Sun;
366*annular*
367when the Moon passes inside the disc of the sun
368to leave only a brilliant ring (Latin *annulus*) visible;
369and *total* when the moon is large enough to cover the Sun completely.
370To classify this eclipse we must compare the size of the Sun and Moon
371to the distance between them.
372Since each argument to ``separation()``
373can be an arbitrary measure of spherical longitude and latitude,
374we can provide azimuth and altitude:
375
376>>> print(ephem.separation((sun.az, sun.alt), (moon.az, moon.alt)))
3770:00:00.3
378>>> print("%.8f %.8f %.11f" % (sun.size, moon.size, sun.size - moon.size))
3791892.91210938 1891.85778809 1.05432128906
380
381The Sun's diameter is larger by 1.05′′,
382so placing the Moon at its center
383would leave an annulus of width
3841.05′′ / 2 = 0.52′′
385visible around the Moon's edge.
386But, in fact, the center of the Moon lies 0.48 arc seconds
387towards one edge of the sun —
388not enough to move its edge outside the sun and make a partial eclipse,
389but enough to make a quite lopsided annular eclipse,
390whose annulus is 0.52′′ + 0.48 = 1.00′′
391wide on one side
392and a scant 0.52′′ - 0.48 = 0.04′′ on the other.
393
394The sky positions computed by PyEphem
395take into account the refraction of the atmosphere,
396which bends upwards the images of bodies near the horizon.
397During sunset, for example, the descent of the sun appears to slow
398because the atmosphere bends its image upwards as it approaches the horizon:
399
400>>> gatech.date = '1984/5/31 00:00'   # 20:00 EDT
401>>> sun.compute(gatech)
402>>> for i in range(8):
403...     old_az, old_alt = sun.az, sun.alt
404...     gatech.date += ephem.minute * 5.
405...     sun.compute(gatech)
406...     sep = ephem.separation((old_az, old_alt), (sun.az, sun.alt))
407...     print("%s %s %s" % (gatech.date, sun.alt, sep))
4081984/5/31 00:05:00 6:17:36.8 1:08:48.1
4091984/5/31 00:10:00 5:21:15.6 1:08:36.3
4101984/5/31 00:15:00 4:25:31.6 1:08:20.0
4111984/5/31 00:20:00 3:30:34.2 1:07:56.5
4121984/5/31 00:25:00 2:36:37.8 1:07:22.7
4131984/5/31 00:30:00 1:44:04.6 1:06:32.2
4141984/5/31 00:35:00 0:53:28.7 1:05:17.0
4151984/5/31 00:40:00 0:05:37.8 1:03:28.3
416
417We see that the Sun's apparent angular speed
418indeed decreased as it approached the horizon,
419from around 1°08′ to barely 1°03′ each five minutes.
420
421Since atmospheric refraction varies with temperature and pressure,
422you can improve the accuracy of PyEphem
423by providing these values from a local forecast,
424or at least from average values for your location and season.
425By default an ``Observer`` uses 15°C and 1010 mB,
426the values for these parameters at sea level
427in the standard atmosphere model used in aviation.
428Setting the pressure to zero
429directs PyEphem to simply ignore atmospheric refraction.
430
431Once PyEphem knows your location it can also work out
432when bodies rise, cross your meridian, and set each day.
433These computations can be fairly involved,
434since planets continue their journey among the stars
435even as the rotation of the earth brings them across the sky;
436PyEphem has to internally re-compute their position several times
437before it finds the exact circumstances of rising or setting.
438But this is taken care of automatically,
439leaving you to simply ask:
440
441>>> print(gatech.next_setting(sun))
4421984/5/31 00:42:22
443>>> print("%s %s" % (sun.alt, sun.az))
444-0:15:45.8 297:20:43.7
445
446Functions also exist for finding risings, transits, and —
447just for completeness —
448the moment of “anti-transit” when the object lies along the meridian
449directly under your feet.
450See the section on :any:`transit-rising-setting`
451in the Quick Reference for more details.
452
453Loading Bodies From Catalogues
454------------------------------
455
456So far we have dealt with the planets, the Sun, and the Moon —
457major bodies whose orbits PyEphem already knows in great detail.
458But for minor bodies, like comets and asteroids,
459you must aquire and load the orbital parameters yourself.
460
461Understand that because the major planets constantly perturb
462the other bodies in the solar system, including each other,
463it requires great effort —
464years of observation yielding formulae with dozens or hundreds of terms —
465to predict the position of a body accurately over decades or centuries.
466For a comet or asteroid,
467astronomers find it more convenient
468to describe its orbit as perfect ellipse, parabola, or hyperbola,
469and then issue new orbital parameters as its orbit changes.
470
471The PyEphem home page provides links to several
472:doc:`catalogs` of orbital elements.
473Once you have obtained elements for a particular body,
474simply provide them to PyEphem's ``readdb()`` function
475in `ephem database format`_ and the resulting object is ready to use:
476
477>>> yh = ephem.readdb("C/2002 Y1 (Juels-Holvorcem),e,103.7816," +
478...    "166.2194,128.8232,242.5695,0.0002609,0.99705756,0.0000," +
479...    "04/13.2508/2003,2000,g  6.5,4.0")
480>>> yh.compute('2003/4/11')
481>>> print(yh.name)
482C/2002 Y1 (Juels-Holvorcem)
483>>> print("%s %s" % (yh.ra, yh.dec))
4840:22:44.58 26:49:48.1
485>>> print("%s %s" % (ephem.constellation(yh), yh.mag))
486('And', 'Andromeda') 5.96
487
488.. _ephem database format: http://www.clearskyinstitute.com/xephem/help/xephem.html#mozTocId468501
489
490(Unfortunately, the library upon which PyEphem is build
491truncates object names to twenty characters, as you can see.)
492Each call to ``readdb()`` returns an object appropriate
493for the orbit specified in the database entry;
494in this case it has returned an ``EllipticalBody``:
495
496>>> print(yh)  # doctest: +ELLIPSIS
497<ephem.EllipticalBody 'C/2002 Y1 (Juels-Holvorcem)' at 0x...>
498
499For objects for which you cannot find an entry in ephem database format,
500you can always create the appropriate kind of object
501and then fill in its orbital parameters yourself;
502see the :doc:`quick` for their names and meanings.
503By calling the ``writedb()`` function of a PyEphem object,
504you can even get it to generate its own database entry
505for archiving or distribution.
506
507There is one other database format with which PyEphem is familiar:
508the NORAD Two-Line Element format (TLE) used for earth satellites.
509Here are some recent elements for the International Space Station.
510
511>>> iss = ephem.readtle("ISS (ZARYA)",
512...  "1 25544U 98067A   03097.78853147  .00021906  00000-0  28403-3 0  8652",
513...  "2 25544  51.6361  13.7980 0004256  35.6671  59.2566 15.58778559250029")
514>>> gatech.date = '2003/3/23'
515>>> iss.compute(gatech)
516>>> print("%s %s %s" % (iss.rise_time, iss.transit_time, iss.set_time))
5172003/3/23 00:00:50 2003/3/23 00:03:26 2003/3/23 00:06:01
518
519The ``transit_time`` for an artificial satellite is actually
520defined in PyEphem as the moment at which it is at highest altitude,
521not the moment at which it crosses (transits) the local meridian.
522
523Note that earth satellites are fast movers —
524in this case rising and setting in less than six minutes!
525They can therefore have multiple risings and settings each day,
526and the particular ones you get from ``rise_time`` and ``set_time``
527depend on the particular time of day for which you ask.
528Repeating the above query eight hours later gives complete different results:
529
530>>> gatech.date = '2003/3/23 8:00'
531>>> iss.compute(gatech)
532>>> print("%s %s %s" % (iss.rise_time, iss.transit_time, iss.set_time))
5332003/3/23 08:03:40 2003/3/23 08:08:25 2003/3/23 08:13:10
534
535When calling ``compute()`` for an earth satellite
536you should provide an ``Observer``,
537and not simply a date and epoch,
538since its location is entirely dependent
539upon the location from which you are observing.
540PyEphem provides extra information about earth satellites,
541beyond the ones available for other objects;
542again, see the :doc:`quick` for details.
543
544Fixed Objects, Precession, and Epochs
545-------------------------------------
546
547The simplest kind of object to create from a catalog entry
548are *fixed* objects,
549for which a constant right ascension and declination are specified.
550These include stars, nebulae, global clusters, and galaxies.
551One example is Polaris, the North Star,
552which lies at the end of Ursa Minor's tail:
553
554>>> polaris = ephem.readdb("Polaris,f|M|F7,2:31:48.704,89:15:50.72,2.02,2000")
555>>> print(polaris.dec)
556Traceback (most recent call last):
557 ...
558RuntimeError: field dec undefined until first compute()
559
560We are able to create the object successfully —
561why should asking its position raise a runtime error?
562The reason is that fixed objects, like planets,
563have an undefined position and magnitude
564until you call their ``compute()`` method
565to determine their position for a particular date or ``Observer``:
566
567>>> polaris.compute()    # uses the current time by default
568>>> print(polaris.a_dec)
56989:15:50.7
570>>> print(ephem.degrees(ephem.degrees('90') - polaris.a_dec))
5710:44:09.3
572
573Much better; we see that the `North Star` lies
574less than forty-five arc minutes from the pole.
575But why should we have to call ``compute()``
576for something fixed —
577something whose position is considered permanent,
578and which should not move between one date and another?
579
580The reason is that, while `fixed` stars and nebulae
581are indeed nearly motionless over the span of human civilization,
582the coordinate system by which we designate their positions
583changes more rapidly.
584Right ascension and declination are based
585upon the orientation of the earth's pole —
586but it turns out that the pole slowly revolves
587(around the axis of the ecliptic plane)
588like the axis of a whirling top,
589completing each revolution in roughly 25,800 years.
590This motion is called *precession*.
591Because this makes the entire coordinate system shift slightly every year,
592is not sufficient to state that Polaris lies at
5932h31m right ascension and 89:15° declination;
594you have to say in *which year*.
595
596That is why the Polaris entry above ends with ``2000``;
597this gives the year for which the coordinates are correct,
598called the *epoch* of the coordinates.
599Because the year 2000 is currently a very popular epoch
600for quoting positions and orbital parameters,
601``compute()`` uses it by default;
602but we can provide an ``epoch=`` keyword parameter
603to have the coordinates translated into those for another year:
604
605>>> polaris.compute(epoch='2100')
606>>> print(polaris.a_dec)
60789:32:26.1
608
609Thus we see that in another hundred years Polaris
610will actually lie closer to the pole that it does today.
611(The ``'2100'`` is the same year/month/day format you have seen already,
612missing both its month and day
613because we are not bothering to be that specific.)
614If you enter subsequent years you will find
615that 2100 is very nearly the closest approach of the pole to Polaris,
616and that soon afterwards they move apart.
617For much of the twenty-five thousand year journey the pole makes,
618there are no stars very near;
619we may have been lucky to have held the Age of Exploration
620as the pole was approaching as convenient a star as Polaris.
621
622Today a dim star in Draco named Thuban
623lies more than twenty degrees from the pole:
624
625>>> thuban = ephem.readdb("Thuban,f|V|A0,14:4:23.3,64:22:33,3.65,2000")
626>>> thuban.compute()
627>>> print(thuban.a_dec)
62864:22:33.0
629
630But in 2801 BC, as the Egyptians built the pyramids,
631Thuban served as their pole star,
632while Polaris lay further from their pole than Thuban lies from ours today:
633
634>>> thuban.compute(epoch='-2800')
635>>> print(thuban.a_dec)
63689:54:35.0
637>>> polaris.compute(epoch='-2800')
638>>> print(polaris.a_dec)
63963:33:17.6
640
641Realize that in these examples I have been lazy
642by giving ``compute()`` an epoch without an actual date,
643which requests the *current* position of each star
644in the coordinates of another epoch.
645This makes no difference for these fixed objects,
646since their positions never change;
647but when dealing with moving objects
648one must always keep in mind the difference
649between the date for which you want their position computed,
650and the epoch in which you want those coordinates expressed.
651Here are some example ``compute()`` calls,
652beginning with one like the above but for a moving object:
653
654``halley.compute(epoch='1066')``
655 This is probably useless:
656 it computes the current position of ``halley``,
657 but returns coordinates relative
658 to the direction the earth's axis was pointing in the year 1066.
659 Unless you use a Conquest-era star atlas, this is not useful.
660
661``halley.compute('1066', epoch='1066')``
662 This is slightly more promising:
663 it computes the position of ``halley`` in 1066
664 and returns coordinates for the orientation of the earth in that year.
665 This might help you visualize
666 how the object was positioned above contemporary observers,
667 who considered it an ill omen in the imminent conflict
668 between King Harold of England and William the Bastard.
669 But to plot this position against a background of stars,
670 you would first have to recompute each star's position in 1066 coordinates.
671
672``halley.compute('1066')``
673 This is what you will probably use most often;
674 you get the position of ``halley`` in the year 1066
675 but expressed in the 2000 coordinates that your star atlas probably uses.
676
677When planning to observe with an equatorial telescope,
678you may want to use the current date as your epoch,
679because the rotation of the sky above your telescope
680is determined by where the pole points today,
681not where it pointed in 2000 or some other convenient epoch.
682Computing positions in the epoch of their date
683is accomplished by simply providing the same argument for both date and epoch:
684
685>>> j = ephem.Jupiter()
686>>> j.compute(epoch=ephem.now())   # so both date and epoch are now
687>>> print("%s %s" % (j.a_ra, j.a_dec))  # doctest: +SKIP
6888:44:29.49 19:00:10.23
689>>> j.compute('2003/3/25', epoch='2003/3/25')
690>>> print("%s %s" % (j.a_ra, j.a_dec))
6918:43:32.82 19:03:32.5
692
693Be careful when computing distances;
694comparing two positions in the coordinates of their own epochs
695will give slightly different results
696than if the two were based on the same epoch:
697
698>>> j1, j2 = ephem.Jupiter(), ephem.Jupiter()
699>>> j1.compute('2003/3/1')
700>>> j2.compute('2003/4/1')
701>>> print(ephem.separation(
702...     (j1.a_ra, j1.a_dec),
703...     (j2.a_ra, j2.a_dec)))   # coordinates are both epoch 2000
7041:46:35.9
705>>> j1.compute('2003/3/1', '2003/3/1')
706>>> j2.compute('2003/4/1', '2003/4/1')
707>>> print(ephem.separation(
708...     (j1.a_ra, j1.a_dec),
709...     (j2.a_ra, j2.a_dec)))   # coordinates are both epoch-of-date
7101:46:31.6
711
712Comparing coordinates of the same epoch, as in the first call above,
713measures motion against the background of stars;
714comparing coordinates from different epochs, as in the second call,
715measures motion against the slowly shifting coordinate system of the earth.
716Users are most often interested in the first kind of measurement,
717and stick with a single epoch the whole way through a computation.
718
719It was for the sake of simplicity
720that all of the examples in this section
721simply provided dates as arguments to the ``compute()`` function.
722If you are instead using an ``Observer`` argument,
723then you specify the epoch through the observer's ``epoch`` variable,
724not through the ``epoch=`` argument.
725Observers use epoch 2000 by default.
726
727Finally,
728make sure you understand
729that your choice of epoch only affects absolute position —
730the right ascension and declination returned for objects —
731*not* the azimuth and altitude of an object above an observer.
732This is because the sun will hang in the same position over Atlanta
733whether the star atlas with which you plot its position
734has epoch 2000, or 1950, or even 1066 coordinates;
735the epoch only affects how you name locations in the sky,
736not how they are positioned with respect to you.
737