1# The core functionalty of PyEphem lives in the C-language _libastro 2# module, which packages the astronomy routines from XEphem as 3# convenient Python types. 4 5import ephem._libastro as _libastro 6import re 7from datetime import datetime as _datetime 8from datetime import timedelta as _timedelta 9from datetime import tzinfo as _tzinfo 10from math import acos, cos, pi, sin 11from time import localtime as _localtime 12 13__version__ = '4.1.3' 14 15# As a favor, compile a regular expression that our C library would 16# really rather not compile for itself. 17 18_libastro._scansexa_split = re.compile(r''' 19 \s*:\s* # A colon optionally surrounded by whitespace, 20 | # or, 21 (?<!^)\s+(?!$) # whitespace not at the start or end of the string. 22''', re.X).split 23 24# Various constants. 25 26tau = 6.283185307179586476925287 27twopi = pi * 2. 28halfpi = pi / 2. 29quarterpi = pi / 4. 30eighthpi = pi / 8. 31 32degree = pi / 180. 33arcminute = degree / 60. 34arcsecond = arcminute / 60. 35half_arcsecond = arcsecond / 2. 36tiny = arcsecond / 360. 37 38c = 299792458. # exact speed of light in meters/second 39meters_per_au = _libastro.meters_per_au 40earth_radius = _libastro.earth_radius 41moon_radius = _libastro.moon_radius 42sun_radius = _libastro.sun_radius 43 44B1900 = 2415020.3135 - _libastro.MJD0 45B1950 = 2433282.4235 - _libastro.MJD0 46J2000 = _libastro.J2000 47 48_slightly_less_than_zero = -1e-15 49_slightly_more_than_pi = pi + 1e-15 50 51# We make available several basic types from _libastro. 52 53Angle = _libastro.Angle 54degrees = _libastro.degrees 55hours = _libastro.hours 56 57Date = _libastro.Date 58hour = 1. / 24. 59minute = hour / 60. 60second = minute / 60. 61 62default_newton_precision = second / 10. 63 64delta_t = _libastro.delta_t 65julian_date = _libastro.julian_date 66 67Body = _libastro.Body 68Planet = _libastro.Planet 69PlanetMoon = _libastro.PlanetMoon 70FixedBody = _libastro.FixedBody 71EllipticalBody = _libastro.EllipticalBody 72ParabolicBody = _libastro.ParabolicBody 73HyperbolicBody = _libastro.HyperbolicBody 74EarthSatellite = _libastro.EarthSatellite 75 76readdb = _libastro.readdb 77readtle = _libastro.readtle 78constellation = _libastro.constellation 79separation = _libastro.separation 80unrefract = _libastro.unrefract 81now = _libastro.now 82 83millennium_atlas = _libastro.millennium_atlas 84uranometria = _libastro.uranometria 85uranometria2000 = _libastro.uranometria2000 86 87# We also create a Python class ("Mercury", "Venus", etcetera) for 88# each planet and moon for which _libastro offers specific algorithms. 89 90for index, classname, name in _libastro.builtin_planets(): 91 exec(''' 92class %(name)s(_libastro.%(classname)s): 93 "Create a Body instance representing %(name)s" 94 __planet__ = %(index)r 95''' % dict(name=name, classname=classname, index=index)) 96 97del index, classname, name 98 99# We now replace two of the classes we have just created, because 100# _libastro actually provides separate types for two of the bodies. 101 102Jupiter = _libastro.Jupiter 103Saturn = _libastro.Saturn 104Moon = _libastro.Moon 105 106# Angles. 107 108def _plusminus_pi(angle): 109 return (angle - pi) % tau - pi 110 111# Newton's method. 112 113def newton(f, x0, x1, precision=default_newton_precision): 114 """Return an x-value at which the given function reaches zero. 115 116 Stops and declares victory once the x-value is within ``precision`` 117 of the solution, which defaults to a half-second of clock time. 118 119 """ 120 f0, f1 = f(x0), f(x1) 121 while f1 and abs(x1 - x0) > precision and f1 != f0: 122 x0, x1 = x1, x1 + (x1 - x0) / (f0/f1 - 1) 123 f0, f1 = f1, f(x1) 124 return x1 125 126# Find equinoxes and solstices. 127 128_sun = Sun() # used for computing equinoxes 129 130def holiday(d0, motion, offset): 131 """Function that assists the finding of equinoxes and solstices.""" 132 133 def f(d): 134 _sun.compute(d) 135 return (_sun.ra + eighthpi) % quarterpi - eighthpi 136 d0 = Date(d0) 137 _sun.compute(d0) 138 angle_to_cover = motion - (_sun.ra + offset) % motion 139 if abs(angle_to_cover) < tiny: 140 angle_to_cover = motion 141 d = d0 + 365.25 * angle_to_cover / twopi 142 return date(newton(f, d, d + hour)) 143 144def previous_vernal_equinox(date): 145 """Return the date of the previous vernal equinox.""" 146 return holiday(date, -twopi, 0) 147 148def next_vernal_equinox(date): 149 """Return the date of the next vernal equinox.""" 150 return holiday(date, twopi, 0) 151 152def previous_summer_solstice(date): 153 """Return the date of the previous summer solstice.""" 154 return holiday(date, -twopi, pi + halfpi) 155 156def next_summer_solstice(date): 157 """Return the date of the next summer solstice.""" 158 return holiday(date, twopi, pi + halfpi) 159 160def previous_autumnal_equinox(date): 161 """Return the date of the previous autumnal equinox.""" 162 return holiday(date, -twopi, pi) 163 164def next_autumnal_equinox(date): 165 """Return the date of the next autumnal equinox.""" 166 return holiday(date, twopi, pi) 167 168def previous_winter_solstice(date): 169 """Return the date of the previous winter solstice.""" 170 return holiday(date, -twopi, halfpi) 171 172def next_winter_solstice(date): 173 """Return the date of the next winter solstice.""" 174 return holiday(date, twopi, halfpi) 175 176# Common synonyms. 177 178next_spring_equinox = next_vernal_equinox 179previous_spring_equinox = previous_vernal_equinox 180 181next_fall_equinox = next_autumn_equinox = next_autumnal_equinox 182previous_fall_equinox = previous_autumn_equinox = previous_autumnal_equinox 183 184# More-general functions that find any equinox or solstice. 185 186def previous_equinox(date): 187 """Return the date of the previous equinox.""" 188 return holiday(date, -pi, 0) 189 190def next_equinox(date): 191 """Return the date of the next equinox.""" 192 return holiday(date, pi, 0) 193 194def previous_solstice(date): 195 """Return the date of the previous solstice.""" 196 return holiday(date, -pi, halfpi) 197 198def next_solstice(date): 199 """Return the date of the next solstice.""" 200 return holiday(date, pi, halfpi) 201 202# Find phases of the Moon. 203 204_moon = Moon() # used for computing Moon phases 205 206def _find_moon_phase(d0, motion, target): 207 """Function that assists the finding of moon phases.""" 208 209 def f(d): 210 _sun.compute(d) 211 _moon.compute(d) 212 slon = _libastro.eq_ecl(d, _sun.g_ra, _sun.g_dec)[0] 213 mlon = _libastro.eq_ecl(d, _moon.g_ra, _moon.g_dec)[0] 214 return (mlon - slon - antitarget) % twopi - pi 215 antitarget = target + pi 216 d0 = Date(d0) 217 f0 = f(d0) 218 angle_to_cover = (- f0) % motion 219 if abs(angle_to_cover) < tiny: 220 angle_to_cover = motion 221 d = d0 + 29.53 * angle_to_cover / twopi 222 return date(newton(f, d, d + hour)) 223 224def previous_new_moon(date): 225 """Return the date of the previous New Moon.""" 226 return _find_moon_phase(date, -twopi, 0) 227 228def next_new_moon(date): 229 """Return the date of the next New Moon.""" 230 return _find_moon_phase(date, twopi, 0) 231 232def previous_first_quarter_moon(date): 233 """Return the date of the previous First Quarter Moon.""" 234 return _find_moon_phase(date, -twopi, halfpi) 235 236def next_first_quarter_moon(date): 237 """Return the date of the next First Quarter Moon.""" 238 return _find_moon_phase(date, twopi, halfpi) 239 240def previous_full_moon(date): 241 """Return the date of the previous Full Moon.""" 242 return _find_moon_phase(date, -twopi, pi) 243 244def next_full_moon(date): 245 """Return the date of the next Full Moon.""" 246 return _find_moon_phase(date, twopi, pi) 247 248def previous_last_quarter_moon(date): 249 """Return the date of the previous Last Quarter Moon.""" 250 return _find_moon_phase(date, -twopi, pi + halfpi) 251 252def next_last_quarter_moon(date): 253 """Return the date of the next Last Quarter Moon.""" 254 return _find_moon_phase(date, twopi, pi + halfpi) 255 256# We provide a Python extension to our _libastro "Observer" class that 257# can search for circumstances like transits. 258 259class CircumpolarError(ValueError): pass 260class NeverUpError(CircumpolarError): pass 261class AlwaysUpError(CircumpolarError): pass 262 263def describe_riset_search(method): 264 if method.__doc__ is None: 265 return method 266 267 method.__doc__ += """, returning its date. 268 269 The search starts at the `date` of this `Observer` and is limited to 270 the single circuit of the sky, from antitransit to antitransit, that 271 the `body` was in the middle of describing at that date and time. 272 If the body did not, in fact, cross the horizon in the direction you 273 are asking about during that particular circuit, then the search 274 must raise a `CircumpolarError` exception like `NeverUpError` or 275 `AlwaysUpError` instead of returning a date. 276 277 """ 278 return method 279 280class Observer(_libastro.Observer): 281 """A location on earth for which positions are to be computed. 282 283 An `Observer` instance allows you to compute the positions of 284 celestial bodies as seen from a particular latitude and longitude on 285 the Earth's surface. The constructor takes no parameters; instead, 286 set its attributes once you have created it. Defaults: 287 288 `date` - the moment the `Observer` is created 289 `lat` - zero latitude 290 `lon` - zero longitude 291 `elevation` - 0 meters above sea level 292 `horizon` - 0 degrees 293 `epoch` - J2000 294 `temp` - 15 degrees Celsius 295 `pressure` - 1010 mBar 296 297 """ 298 __slots__ = [ 'name' ] 299 elev = _libastro.Observer.elevation 300 301 def copy(self): 302 o = self.__class__() 303 o.date = self.date 304 o.lat = self.lat 305 o.lon = self.lon 306 o.elev = self.elev 307 o.horizon = self.horizon 308 o.epoch = self.epoch 309 o.temp = self.temp 310 o.pressure = self.pressure 311 return o 312 313 __copy__ = copy 314 315 def __repr__(self): 316 """Return a useful textual representation of this Observer.""" 317 return ('<ephem.Observer date=%r epoch=%r' 318 " lon='%s' lat='%s' elevation=%sm" 319 ' horizon=%s temp=%sC pressure=%smBar>' 320 % (str(self.date), str(self.epoch), 321 self.lon, self.lat, self.elevation, 322 self.horizon, self.temp, self.pressure)) 323 324 def compute_pressure(self): 325 """Set the atmospheric pressure for the current elevation.""" 326 # Formula from the ISA Standard Atmosphere 327 self.pressure = (1013.25 * (1 - 0.0065 * self.elevation / 288.15) 328 ** 5.2558761132785179) 329 330 def _compute_transit(self, body, start, sign, offset): 331 """Internal function used to compute transits.""" 332 333 if isinstance(body, EarthSatellite): 334 raise TypeError( 335 'the next and previous transit methods do not' 336 ' support earth satellites because of their speed;' 337 ' please use the higher-resolution next_pass() method' 338 ) 339 340 def f(d): 341 self.date = d 342 body.compute(self) 343 return degrees(offset - sidereal_time() + body.g_ra).znorm 344 345 if start is not None: 346 self.date = start 347 sidereal_time = self.sidereal_time 348 body.compute(self) 349 ha = sidereal_time() - body.g_ra 350 ha_to_move = (offset - ha) % (sign * twopi) 351 if abs(ha_to_move) < tiny: 352 ha_to_move = sign * twopi 353 d = self.date + ha_to_move / twopi 354 result = Date(newton(f, d, d + minute)) 355 return result 356 357 def _previous_transit(self, body, start=None): 358 """Find the previous passage of a body across the meridian.""" 359 360 return self._compute_transit(body, start, -1., 0.) 361 362 def _next_transit(self, body, start=None): 363 """Find the next passage of a body across the meridian.""" 364 365 return self._compute_transit(body, start, +1., 0.) 366 367 def _previous_antitransit(self, body, start=None): 368 """Find the previous passage of a body across the anti-meridian.""" 369 370 return self._compute_transit(body, start, -1., pi) 371 372 def _next_antitransit(self, body, start=None): 373 """Find the next passage of a body across the anti-meridian.""" 374 375 return self._compute_transit(body, start, +1., pi) 376 377 def previous_transit(self, body, start=None): 378 """Find the previous passage of a body across the meridian.""" 379 380 original_date = self.date 381 d = self._previous_transit(body, start) 382 self.date = original_date 383 return d 384 385 def next_transit(self, body, start=None): 386 """Find the next passage of a body across the meridian.""" 387 388 original_date = self.date 389 d = self._next_transit(body, start) 390 self.date = original_date 391 return d 392 393 def previous_antitransit(self, body, start=None): 394 """Find the previous passage of a body across the anti-meridian.""" 395 396 original_date = self.date 397 d = self._previous_antitransit(body, start) 398 self.date = original_date 399 return d 400 401 def next_antitransit(self, body, start=None): 402 """Find the next passage of a body across the anti-meridian.""" 403 404 original_date = self.date 405 d = self._next_antitransit(body, start) 406 self.date = original_date 407 return d 408 409 def disallow_circumpolar(self, declination): 410 """Raise an exception if the given declination is circumpolar. 411 412 Raises NeverUpError if an object at the given declination is 413 always below this Observer's horizon, or AlwaysUpError if such 414 an object would always be above the horizon. 415 416 """ 417 if abs(self.lat - declination) >= halfpi: 418 raise NeverUpError('The declination %s never rises' 419 ' above the horizon at latitude %s' 420 % (declination, self.lat)) 421 if abs(self.lat + declination) >= halfpi: 422 raise AlwaysUpError('The declination %s is always' 423 ' above the horizon at latitude %s' 424 % (declination, self.lat)) 425 426 @describe_riset_search 427 def previous_rising(self, body, start=None, use_center=False): 428 """Search for the given body's previous rising""" 429 return self._find_rise_or_set(body, start, use_center, -1, True) 430 431 @describe_riset_search 432 def previous_setting(self, body, start=None, use_center=False): 433 """Search for the given body's previous setting""" 434 return self._find_rise_or_set(body, start, use_center, -1, False) 435 436 @describe_riset_search 437 def next_rising(self, body, start=None, use_center=False): 438 """Search for the given body's next rising""" 439 return self._find_rise_or_set(body, start, use_center, +1, True) 440 441 @describe_riset_search 442 def next_setting(self, body, start=None, use_center=False): 443 """Search for the given body's next setting""" 444 return self._find_rise_or_set(body, start, use_center, +1, False) 445 446 def _find_rise_or_set(self, body, start, use_center, direction, do_rising): 447 if isinstance(body, EarthSatellite): 448 raise TypeError( 449 'the rising and settings methods do not' 450 ' support earth satellites because of their speed;' 451 ' please use the higher-resolution next_pass() method' 452 ) 453 454 original_pressure = self.pressure 455 original_date = self.date 456 try: 457 self.pressure = 0.0 # otherwise geometry doesn't work 458 if start is not None: 459 self.date = start 460 prev_ha = None 461 while True: 462 body.compute(self) 463 horizon = self.horizon 464 if not use_center: 465 horizon -= body.radius 466 if original_pressure: 467 horizon = unrefract(original_pressure, self.temp, horizon) 468 abs_target_ha = self._target_hour_angle(body, horizon) 469 if do_rising: 470 target_ha = - abs_target_ha # rises in east (az 0-180) 471 else: 472 target_ha = abs_target_ha # sets in west (az 180-360) 473 ha = body.ha 474 difference = target_ha - ha 475 if prev_ha is None: 476 difference %= tau # force angle to be positive 477 if direction < 0: 478 difference -= tau 479 bump = difference / tau 480 if abs(bump) < default_newton_precision: 481 # Already at target event: move forward to next one. 482 bump += direction 483 else: 484 difference = _plusminus_pi(difference) 485 bump = difference / tau 486 if abs(bump) < default_newton_precision: 487 break 488 self.date += bump 489 prev_ha = ha 490 491 if abs_target_ha == _slightly_more_than_pi: 492 raise AlwaysUpError('%r is above the horizon at %s' 493 % (body.name, self.date)) 494 495 if abs_target_ha == _slightly_less_than_zero: 496 raise NeverUpError('%r is below the horizon at %s' 497 % (body.name, self.date)) 498 499 return self.date 500 finally: 501 if self.pressure != original_pressure: 502 self.pressure = original_pressure 503 body.compute(self) 504 self.date = original_date 505 506 def _target_hour_angle(self, body, alt): 507 lat = self.lat 508 dec = body.dec 509 arg = (sin(alt) - sin(lat) * sin(dec)) / (cos(lat) * cos(dec)) 510 511 if arg < -1.0: 512 return _slightly_more_than_pi 513 elif arg > 1.0: 514 return _slightly_less_than_zero 515 516 return acos(arg) 517 518 def next_pass(self, body, singlepass=True): 519 """Return the next rising, culmination, and setting of a satellite. 520 521 If singlepass is True, return next consecutive set of 522 ``(rising, culmination, setting)``. 523 524 If singlepass is False, return 525 ``(next_rising, next_culmination, next_setting)``. 526 527 """ 528 if not isinstance(body, EarthSatellite): 529 raise TypeError( 530 'the next_pass() method is only for use with' 531 ' EarthSatellite objects because of their high speed' 532 ) 533 534 result = _libastro._next_pass(self, body) 535 # _libastro behavior is singlepass=False 536 if ((not singlepass) 537 or (None in result) 538 or (result[4] >= result[0])): 539 return result 540 # retry starting just before next_rising 541 obscopy = self.copy() 542 # Almost always 1 minute before next_rising except 543 # in pathological case where set came immediately before rise 544 obscopy.date = result[0] - min(1.0/1440, 545 (result[0] - result[4])/2) 546 result = _libastro._next_pass(obscopy, body) 547 if result[0] <= result[2] <= result[4]: 548 return result 549 raise ValueError("this software is having trouble with those satellite parameters") 550 551 552del describe_riset_search 553 554# Time conversion. 555 556def _convert_to_seconds_and_microseconds(date): 557 """Converts a PyEphem date into seconds""" 558 microseconds = int(round(24 * 60 * 60 * 1000000 * date)) 559 seconds, microseconds = divmod(microseconds, 1000000) 560 seconds -= 2209032000 # difference between epoch 1900 and epoch 1970 561 return seconds, microseconds 562 563 564def localtime(date): 565 """Convert a PyEphem date into naive local time, returning a Python datetime.""" 566 seconds, microseconds = _convert_to_seconds_and_microseconds(date) 567 y, m, d, H, M, S, wday, yday, isdst = _localtime(seconds) 568 return _datetime(y, m, d, H, M, S, microseconds) 569 570 571class _UTC(_tzinfo): 572 ZERO = _timedelta(0) 573 def utcoffset(self, dt): 574 return self.ZERO 575 def dst(self, dt): 576 return self.ZERO 577 def __repr__(self): 578 return "<ephem.UTC>" 579 580 581UTC = _UTC() 582 583 584def to_timezone(date, tzinfo): 585 """"Convert a PyEphem date into a timezone aware Python datetime representation.""" 586 seconds, microseconds = _convert_to_seconds_and_microseconds(date) 587 date = _datetime.fromtimestamp(seconds, tzinfo) 588 date = date.replace(microsecond=microseconds) 589 return date 590 591# Coordinate transformations. 592 593class Coordinate(object): 594 def __init__(self, *args, **kw): 595 596 # Accept an optional "epoch" keyword argument. 597 598 epoch = kw.pop('epoch', None) 599 if epoch is not None: 600 self.epoch = epoch = Date(epoch) 601 if kw: 602 raise TypeError('"epoch" is the only keyword argument' 603 ' you can use during %s instantiation' 604 % (type(self).__name__)) 605 606 # Interpret a single-argument initialization. 607 608 if len(args) == 1: 609 a = args[0] 610 611 if isinstance(a, Body): 612 a = Equatorial(a.a_ra, a.a_dec, epoch = a.a_epoch) 613 614 for cls in (Equatorial, Ecliptic, Galactic): 615 if isinstance(a, cls): 616 617 # If the user omitted an "epoch" keyword, then 618 # use the epoch of the other object. 619 620 if epoch is None: 621 self.epoch = epoch = a.epoch 622 623 # If we are initialized from another of the same 624 # kind of coordinate and epoch, simply copy the 625 # coordinates and epoch into this new object. 626 627 if isinstance(self, cls) and epoch == a.epoch: 628 self.set(*a.get()) 629 return 630 631 # Otherwise, convert. 632 633 ra, dec = a.to_radec() 634 if epoch != a.epoch: 635 ra, dec = _libastro.precess( 636 a.epoch, epoch, ra, dec 637 ) 638 self.from_radec(ra, dec) 639 return 640 641 raise TypeError( 642 'a single argument used to initialize %s() must be either' 643 ' a coordinate or a Body, not an %r' % (type(a).__name__,) 644 ) 645 646 # Two arguments are interpreted as (ra, dec) or (lon, lat). 647 648 elif len(args) == 2: 649 self.set(*args) 650 if epoch is None: 651 self.epoch = epoch = Date(J2000) 652 653 else: 654 raise TypeError( 655 'to initialize %s you must pass either a Body,' 656 ' another coordinate, or two coordinate values,' 657 ' but not: %r' % (type(self).__name__, args,) 658 ) 659 660class Equatorial(Coordinate): 661 """An equatorial sky coordinate in right ascension and declination.""" 662 663 def get(self): 664 return self.ra, self.dec 665 666 def set(self, ra, dec): 667 self.ra, self.dec = hours(ra), degrees(dec) 668 669 to_radec = get 670 from_radec = set 671 672class LonLatCoordinate(Coordinate): 673 """A coordinate that is measured with a longitude and latitude.""" 674 675 def set(self, lon, lat): 676 self.lon, self.lat = degrees(lon), degrees(lat) 677 678 def get(self): 679 return self.lon, self.lat 680 681 @property 682 def long(self): 683 return self.lon 684 685 @long.setter 686 def long(self, value): 687 self.lon = value 688 689class Ecliptic(LonLatCoordinate): 690 """An ecliptic latitude and longitude.""" 691 692 def to_radec(self): 693 return _libastro.ecl_eq(self.epoch, self.lon, self.lat) 694 695 def from_radec(self, ra, dec): 696 self.lon, self.lat = _libastro.eq_ecl(self.epoch, ra, dec) 697 698class Galactic(LonLatCoordinate): 699 """A galactic latitude and longitude.""" 700 701 def to_radec(self): 702 return _libastro.gal_eq(self.epoch, self.lon, self.lat) 703 704 def from_radec(self, ra, dec): 705 self.lon, self.lat = _libastro.eq_gal(self.epoch, ra, dec) 706 707# For backwards compatibility, provide lower-case names for our Date 708# and Angle classes, and also allow "Lon" to be spelled "Long". 709 710date = Date 711angle = Angle 712LongLatCoordinate = LonLatCoordinate 713 714# Catalog boostraps. Each of these functions imports a catalog 715# module, then replaces itself with the function of the same name that 716# lives inside of the catalog. 717 718def star(name, *args, **kwargs): 719 """Load the stars database and return a star.""" 720 global star 721 import ephem.stars 722 star = ephem.stars.star 723 return star(name, *args, **kwargs) 724 725def city(name): 726 """Load the cities database and return a city.""" 727 global city 728 import ephem.cities 729 city = ephem.cities.city 730 return city(name) 731