1# -*- coding: utf-8 -*- 2"""Classes representing different kinds of astronomical position.""" 3 4from numpy import array, einsum, full, reshape, nan, nan_to_num 5from . import framelib 6from .constants import ANGVEL, AU_M, C, ERAD, DAY_S, RAD2DEG, pi, tau 7from .descriptorlib import reify 8from .earthlib import compute_limb_angle 9from .functions import ( 10 _T, _to_array, _to_spherical_and_rates, angle_between, from_spherical, 11 length_of, mxm, mxv, rot_z, to_spherical, 12) 13from .geometry import intersect_line_and_sphere 14from .relativity import add_aberration, add_deflection 15from .timelib import Time 16from .units import Angle, AngleRate, Distance, Velocity, _interpret_angle 17 18_GIGAPARSEC_AU = 206264806247096.38 # 1e9 * 360 * 3600 / tau 19 20def build_position(position_au, velocity_au_per_d=None, t=None, 21 center=None, target=None): 22 if center == 0: 23 cls = Barycentric 24 elif center == 399: 25 cls = Geocentric 26 elif hasattr(center, 'rotation_at'): # and thus deserves an altaz() method 27 cls = Geometric 28 else: 29 cls = ICRF 30 return cls(position_au, velocity_au_per_d, t, center, target) 31 32def position_of_radec(ra_hours, dec_degrees, distance_au=_GIGAPARSEC_AU, 33 epoch=None, t=None, center=None, target=None): 34 """Build a position object from a right ascension and declination. 35 36 If a specific ``distance_au`` is not provided, Skyfield returns a 37 position vector a gigaparsec in length. This puts the position at a 38 great enough distance that it will stand at the same right ascension 39 and declination from any viewing position in the Solar System, to 40 very high precision (within a few hundredths of a microarcsecond). 41 42 If an ``epoch`` is specified, the input coordinates are understood 43 to be in the dynamical system of that particular date. Otherwise, 44 they will be assumed to be ICRS (the modern replacement for J2000). 45 46 .. versionadded:: 1.21 47 This replaces a deprecated function ``position_from_radec()`` 48 whose ``distance`` argument was not as well designed. 49 50 """ 51 theta = _to_array(dec_degrees) / 360.0 * tau 52 phi = _to_array(ra_hours) / 24.0 * tau 53 position_au = from_spherical(distance_au, theta, phi) 54 if epoch is not None: 55 position_au = mxv(epoch.MT, position_au) 56 return build_position(position_au, None, t, center, target) 57 58def position_from_radec(ra_hours, dec_degrees, distance=1.0, epoch=None, 59 t=None, center=None, target=None): 60 """DEPRECATED version of ``position_of_radec()``. 61 62 Problems: 63 64 * The ``distance`` parameter specifies no unit, contrary to Skyfield 65 best practices. I have no idea what I was thinking. 66 67 * The default ``distance`` is far too small, since most objects for 68 which users specify an RA and declination are out on the celestial 69 sphere. The hope was that users would see the length 1.0 and 70 think, “ah, yes, that’s obviously a fake placeholder value.” But 71 it’s more likely that users will not even check the distance, or 72 maybe not even realize that a distance is involved. 73 74 """ 75 return position_of_radec(ra_hours, dec_degrees, distance, epoch, 76 t, center, target) 77 78class ICRF(object): 79 """An |xyz| position and velocity oriented to the ICRF axes. 80 81 The International Coordinate Reference Frame (ICRF) is a permanent 82 reference frame that is the replacement for J2000. Their axes agree 83 to within 0.02 arcseconds. It also supersedes older equinox-based 84 systems like B1900 and B1950. 85 86 Each instance of this class provides a ``.position`` vector and a 87 ``.velocity`` vector that specify |xyz| coordinates along the axes 88 of the ICRF. A specific time ``.t`` might be specified or might be 89 ``None``. 90 91 """ 92 center_barycentric = None 93 _observer_gcrs_au = None 94 _default_center = None 95 _ephemeris = None # cached so we can compute how light is deflected 96 97 def __init__(self, position_au, velocity_au_per_d=None, t=None, 98 center=None, target=None): 99 self.t = t 100 self.position = self.xyz = Distance(position_au) 101 if velocity_au_per_d is None: 102 velocity_au_per_d = full(self.position.au.shape, nan) 103 self.velocity = Velocity(velocity_au_per_d) 104 self.center = self._default_center if center is None else center 105 self.target = target 106 if center == 0: 107 self.center_barycentric = self 108 109 @classmethod 110 def from_radec(cls, ra_hours, dec_degrees, 111 distance_au=_GIGAPARSEC_AU, epoch=None): 112 theta = _to_array(dec_degrees) / 360.0 * tau 113 phi = _to_array(ra_hours) / 24.0 * tau 114 position_au = from_spherical(distance_au, theta, phi) 115 if epoch is not None: 116 position_au = mxv(epoch.MT, position_au) 117 return cls(position_au) 118 119 @classmethod 120 def from_time_and_frame_vectors(cls, t, frame, distance, velocity): 121 """Constructor: build a position from two vectors in a reference frame. 122 123 * ``t`` — The :class:`~skyfield.timelib.Time` of the position. 124 * ``frame`` — A reference frame listed at `reference_frames`. 125 * ``distance`` — A `Distance` |xyz| vector in the given frame. 126 * ``velocity`` — A `Velocity` ẋ,ẏ,ż vector in the given frame. 127 128 """ 129 r, v = distance.au, velocity.au_per_d 130 at = getattr(frame, '_dRdt_times_RT_at', None) 131 if at is not None: 132 V = at(t) 133 v = v - mxv(V, r) # subtract instead of transposing 134 RT = _T(frame.rotation_at(t)) 135 r = mxv(RT, r) 136 v = mxv(RT, v) 137 return cls(r, v, t) # TODO: args for center and target? 138 139 def __repr__(self): 140 name = self.__class__.__name__ 141 center = self.center 142 if name == 'Barycentric' and center == 0: 143 suffix = ' BCRS' 144 elif name == 'Apparent' and center == 399: 145 suffix = ' GCRS' 146 elif name != 'ICRF': 147 suffix = ' ICRS' 148 else: 149 suffix = '' 150 151 center = self.center 152 target = self.target 153 154 center_name = getattr(center, 'target_name', None) 155 if center_name is None: 156 center_name = str(center) 157 158 target_name = getattr(target, 'target_name', None) 159 if target_name is None: 160 target_name = str(target) 161 162 return '<{0}{1} position{2}{3}{4}{5}>'.format( 163 name, 164 suffix, 165 '' if (self.velocity is None) else ' and velocity', 166 '' if self.t is None else ' at date t', 167 '' if self.center is None else ' center={0}'.format(center_name), 168 '' if self.target is None else ' target={0}'.format(target_name), 169 ) 170 171 def __sub__(self, body): 172 """Subtract two ICRF vectors to produce a third.""" 173 p = self.position.au - body.position.au 174 if self.velocity is None or body.velocity is None: 175 v = None 176 else: 177 v = self.velocity.au_per_d - body.velocity.au_per_d 178 return build_position(p, v, self.t, body.target, self.target) 179 180 def __getitem__(self, i): 181 return type(self)( 182 self.position.au[:,i], 183 self.velocity.au_per_d[:,i], 184 self.t[i], 185 self.center, 186 self.target, 187 ) 188 189 def __neg__(self): 190 return type(self)( 191 -self.position.au, 192 -self.velocity.au_per_d, 193 self.t, 194 self.target, 195 self.center, 196 ) 197 198 def distance(self): 199 """Compute the distance from the origin to this position. 200 201 The return value is a :class:`~skyfield.units.Distance` that 202 prints itself out in astronomical units (au) but that also 203 offers attributes ``au``, ``km``, and ``m`` if you want to 204 access its magnitude as a number. 205 206 >>> v = ICRF([1, 1, 0]) 207 >>> print(v.distance()) 208 1.41421 au 209 210 """ 211 return Distance(length_of(self.position.au)) 212 213 def speed(self): 214 """Compute the magnitude of the velocity vector. 215 216 >>> v = ICRF([0, 0, 0], [1, 2, 3]) 217 >>> print(v.speed()) 218 3.74166 au/day 219 220 """ 221 return Velocity(length_of(self.velocity.au_per_d)) 222 223 @reify 224 def light_time(self): 225 """Length of this vector in days of light travel time.""" 226 # Note that this property is almost never called, since 227 # .light_time is set manually on the one kind of position 228 # (astrometric) that is ever likely to need it. Alas: back in 229 # 2015, I didn't think to either hide this attribute, or have it 230 # include its units in its name! 231 return self.distance().m / C / DAY_S 232 233 def radec(self, epoch=None): 234 r"""Compute equatorial RA, declination, and distance. 235 236 When called without a parameter, this returns standard ICRF 237 right ascension and declination: 238 239 >>> from skyfield.api import load 240 >>> ts = load.timescale() 241 >>> t = ts.utc(2020, 5, 13, 10, 32) 242 >>> eph = load('de421.bsp') 243 >>> astrometric = eph['earth'].at(t).observe(eph['sun']) 244 >>> ra, dec, distance = astrometric.radec() 245 >>> print(ra, dec, sep='\n') 246 03h 21m 47.67s 247 +18deg 28' 55.3" 248 249 If you instead want the coordinates referenced to the dynamical 250 system defined by the Earth's true equator and equinox, provide 251 a specific epoch time. 252 253 >>> ra, dec, distance = astrometric.apparent().radec(epoch='date') 254 >>> print(ra, dec, sep='\n') 255 03h 22m 54.73s 256 +18deg 33' 04.5" 257 258 To get J2000.0 coordinates, simply pass ``ts.J2000``. 259 260 """ 261 position_au = self.position.au 262 if epoch is not None: 263 if isinstance(epoch, Time): 264 pass 265 elif isinstance(epoch, float): 266 epoch = Time(None, tt=epoch) 267 elif epoch == 'date': 268 epoch = self.t 269 else: 270 raise ValueError('the epoch= must be a Time object,' 271 ' a floating point Terrestrial Time (TT),' 272 ' or the string "date" for epoch-of-date') 273 position_au = mxv(epoch.M, position_au) 274 r_au, dec, ra = to_spherical(position_au) 275 return (Angle(radians=ra, preference='hours'), 276 Angle(radians=dec, signed=True), 277 Distance(r_au)) 278 279 def hadec(self): 280 """Compute hour angle, declination, and distance. 281 282 Returns a tuple of two :class:`~skyfield.units.Angle` objects 283 plus the :class:`~skyfield.units.Distance` to the target. The 284 angles are the hour angle (±12 hours) east or west of your 285 meridian along the ITRS celestial equator, and the declination 286 (±90 degees) above or below it. This only works for positions 287 whose center is a geographic location; otherwise, there is no 288 local meridian from which to measure the hour angle. 289 290 Because this declination is measured from the plane of the 291 Earth’s physical geographic equator, it will be slightly 292 different than the declination returned by ``radec()`` if you 293 have loaded a :ref:`polar motion` file. 294 295 The coordinates are not adjusted for atmospheric refraction near 296 the horizon. 297 298 """ 299 R = framelib.itrs.rotation_at(self.t) 300 r = mxv(R, self.position.au) 301 au, dec, sublongtiude = to_spherical(r) 302 ha = self.center.longitude.radians - sublongtiude 303 ha += pi 304 ha %= tau 305 ha -= pi 306 return (Angle(radians=ha, preference='hours', signed=True), 307 Angle(radians=dec, signed=True), 308 Distance(au)) 309 310 def separation_from(self, another_icrf): 311 """Return the angle between this position and another. 312 313 >>> from skyfield.api import load 314 >>> ts = load.timescale() 315 >>> t = ts.utc(2020, 4, 18) 316 >>> eph = load('de421.bsp') 317 >>> sun, venus, earth = eph['sun'], eph['venus'], eph['earth'] 318 >>> e = earth.at(t) 319 >>> s = e.observe(sun) 320 >>> v = e.observe(venus) 321 >>> print(s.separation_from(v)) 322 43deg 23' 23.1" 323 324 You can also compute separations across an array of positions. 325 326 >>> t = ts.utc(2020, 4, [18, 19, 20]) 327 >>> e = earth.at(t) 328 >>> print(e.observe(sun).separation_from(e.observe(venus))) 329 3 values from 43deg 23' 23.1" to 42deg 49' 46.6" 330 331 """ 332 u = self.position.au 333 v = another_icrf.position.au 334 335 # Allow an array of positions to be compared with a single other 336 # position. 337 difference = len(u.shape) - len(v.shape) 338 if difference: 339 if difference > 0: 340 v = reshape(v, v.shape + (1,) * difference) 341 else: 342 u = reshape(u, u.shape + (1,) * -difference) 343 344 return Angle(radians=angle_between(u, v)) 345 346 # TODO: build a reference frame for the following two methods. 347 348 def cirs_xyz(self, epoch): 349 """Compute cartesian CIRS coordinates at a given epoch |xyz|. 350 351 Calculate coordinates in the Celestial Intermediate Reference System 352 (CIRS), a dynamical coordinate system referenced to the Celestial 353 Intermediate Origin (CIO). As this is a dynamical system it must be 354 calculated at a specific epoch. 355 """ 356 if isinstance(epoch, Time): 357 pass 358 elif isinstance(epoch, float): 359 epoch = Time(None, tt=epoch) 360 elif epoch == 'date': 361 epoch = self.t 362 else: 363 raise ValueError('the epoch= must be a Time object,' 364 ' a floating point Terrestrial Time (TT),' 365 ' or the string "date" for epoch-of-date') 366 367 vector = mxv(epoch.C, self.position.au) 368 return Distance(vector) 369 370 def cirs_radec(self, epoch): 371 """Get spherical CIRS coordinates at a given epoch (ra, dec, distance). 372 373 Calculate coordinates in the Celestial Intermediate Reference System 374 (CIRS), a dynamical coordinate system referenced to the Celestial 375 Intermediate Origin (CIO). As this is a dynamical system it must be 376 calculated at a specific epoch. 377 """ 378 r_au, dec, ra = to_spherical(self.cirs_xyz(epoch).au) 379 380 return (Angle(radians=ra, preference='hours'), 381 Angle(radians=dec, signed=True), 382 Distance(r_au)) 383 384 # Deprecated methods, that have been replaced by `framelib.py` plus 385 # the "frame" methods in the next section. 386 387 def ecliptic_xyz(self, epoch=None): 388 if epoch is None: 389 return self.frame_xyz(framelib.ecliptic_J2000_frame) 390 return _Fake(self, epoch).frame_xyz(framelib.ecliptic_frame) 391 def ecliptic_velocity(self): 392 return self.frame_xyz_and_velocity(framelib.ecliptic_J2000_frame)[1] 393 def ecliptic_latlon(self, epoch=None): 394 if epoch is None: 395 return self.frame_latlon(framelib.ecliptic_J2000_frame) 396 return _Fake(self, epoch).frame_latlon(framelib.ecliptic_frame) 397 398 def galactic_xyz(self): return self.frame_xyz(framelib.galactic_frame) 399 def galactic_latlon(self): return self.frame_latlon(framelib.galactic_frame) 400 ecliptic_position = ecliptic_xyz # old alias 401 galactic_position = galactic_xyz # old alias 402 403 # New methods for converting to and from `framelib.py` reference frames. 404 405 def frame_xyz(self, frame): 406 """Return this position as an |xyz| vector in a reference frame. 407 408 Returns a :class:`~skyfield.units.Distance` object giving the 409 |xyz| of this position in the given ``frame``. See 410 `reference_frames`. 411 412 """ 413 return Distance(mxv(frame.rotation_at(self.t), self.position.au)) 414 415 def frame_xyz_and_velocity(self, frame): 416 """Return |xyz| position and velocity vectors in a reference frame. 417 418 Returns two vectors in the given coordinate ``frame``: a 419 :class:`~skyfield.units.Distance` providing an |xyz| position 420 and a :class:`~skyfield.units.Velocity` giving (xdot,ydot,zdot) 421 velocity. See `reference_frames`. 422 423 """ 424 R = frame.rotation_at(self.t) 425 r, v = self.position.au, self.velocity.au_per_d 426 r = mxv(R, r) 427 v = mxv(R, v) 428 at = getattr(frame, '_dRdt_times_RT_at', None) 429 if at is not None: 430 V = at(self.t) 431 v += mxv(V, r) 432 return Distance(r), Velocity(v) 433 434 def frame_latlon(self, frame): 435 """Return longitude, latitude, and distance in the given frame. 436 437 Returns a 3-element tuple giving the latitude and longitude as 438 :class:`~skyfield.units.Angle` objects and the range to the 439 target as a :class:`~skyfield.units.Distance`. See 440 `reference_frames`. 441 442 """ 443 vector = mxv(frame.rotation_at(self.t), self.position.au) 444 d, lat, lon = to_spherical(vector) 445 return (Angle(radians=lat, signed=True), 446 Angle(radians=lon), 447 Distance(d)) 448 449 def frame_latlon_and_rates(self, frame): 450 """Return a reference frame longitude, latitude, range, and rates. 451 452 Return a 6-element tuple of 3 coordinates and 3 rates-of-change 453 for this position in the given reference ``frame``: 454 455 * Latitude :class:`~skyfield.units.Angle` from +90° north to −90° south 456 * Longitude :class:`~skyfield.units.Angle` 0°–360° east 457 * Radial :class:`~skyfield.units.Distance` 458 * Latitude :class:`~skyfield.units.AngleRate` 459 * Longitude :class:`~skyfield.units.AngleRate` 460 * Radial :class:`~skyfield.units.Velocity` 461 462 If the reference frame is the ICRS, or is J2000, or otherwise 463 involves the celestial equator and pole, then the latitude and 464 longitude returned will measure what are more commonly called 465 “declination” and “right ascension”. Note that right ascension 466 is usually expressed as hours (24 in a circle), rather than in 467 the degrees that this routine will return. 468 469 """ 470 r, v = self.frame_xyz_and_velocity(frame) 471 r = r.au 472 v = v.au_per_d 473 d, lat, lon, d_rate, lat_rate, lon_rate = _to_spherical_and_rates(r, v) 474 return (Angle(radians=lat, signed=True), 475 Angle(radians=lon), 476 Distance(d), 477 AngleRate._from_radians_per_day(lat_rate), 478 AngleRate._from_radians_per_day(lon_rate), 479 Velocity(d_rate)) 480 481 def to_skycoord(self, unit=None): 482 """Convert this distance to an AstroPy ``SkyCoord`` object.""" 483 from astropy.coordinates import SkyCoord 484 from astropy.units import au 485 x, y, z = self.position.au 486 return SkyCoord(representation_type='cartesian', x=x, y=y, z=z, unit=au) 487 488 def is_sunlit(self, ephemeris): 489 """Return whether a position in Earth orbit is in sunlight. 490 491 Returns ``True`` or ``False``, or an array of such values, to 492 indicate whether this position is in sunlight or is blocked by 493 the Earth’s shadow. It should work with positions produced 494 either by calling ``at()`` on a satellite object, or by calling 495 ``at()`` on the relative position ``sat - topos`` of a satellite 496 with respect to an Earth observer’s position. See 497 :ref:`satellite-is-sunlit`. 498 499 """ 500 if self.center == 399: 501 earth_m = - self.position.m 502 else: 503 gcrs_position = self._observer_gcrs_au 504 if gcrs_position is None: 505 raise ValueError('cannot tell whether this position is sunlit') 506 earth_m = - self.position.m - gcrs_position * AU_M 507 508 sun_m = (ephemeris['sun'] - ephemeris['earth']).at(self.t).position.m 509 near, far = intersect_line_and_sphere(sun_m + earth_m, earth_m, ERAD) 510 return nan_to_num(far) <= 0 511 512 def is_behind_earth(self): 513 """Return whether the Earth blocks the view of this object. 514 515 For a position centered on an Earth-orbiting satellite, return 516 whether the target is in eclipse behind the disc of the Earth. 517 See :ref:`is-behind-earth`. 518 519 """ 520 observer_gcrs_au = self._observer_gcrs_au 521 if observer_gcrs_au is None: 522 raise ValueError('can only compute Earth occultation for' 523 ' positions observed from an Earth satellite') 524 earth_m = - observer_gcrs_au * AU_M 525 vector_m = self.position.m 526 near, far = intersect_line_and_sphere(vector_m, earth_m, ERAD) 527 return nan_to_num(far) > 0 528 529 @reify 530 def _altaz_rotation(self): 531 # Return and cache (with @reify) the orientation of this 532 # observer, in case a single observer.at() position is used in 533 # several subsequent .observe().apparent().altaz() calls. 534 rotation_at = getattr(self.target, 'rotation_at', None) 535 if rotation_at is None: 536 raise ValueError(_altaz_message) 537 return rotation_at(self.t) 538 539 def from_altaz(self, alt=None, az=None, alt_degrees=None, az_degrees=None, 540 distance=Distance(au=0.1)): 541 """Generate an Apparent position from an altitude and azimuth. 542 543 The altitude and azimuth can each be provided as an `Angle` 544 object, or else as a number of degrees provided as either a 545 float or a tuple of degrees, arcminutes, and arcseconds:: 546 547 alt=Angle(...), az=Angle(...) 548 alt_degrees=23.2289, az_degrees=142.1161 549 alt_degrees=(23, 13, 44.1), az_degrees=(142, 6, 58.1) 550 551 The distance should be a :class:`~skyfield.units.Distance` 552 object, if provided; otherwise a default of 0.1 au is used. 553 554 """ 555 # TODO: should this method live on another class? 556 557 rotation_at = getattr(self.target, 'rotation_at', None) 558 if rotation_at is None: 559 raise ValueError(_altaz_message) 560 R = rotation_at(self.t) 561 562 alt = _interpret_angle('alt', alt, alt_degrees) 563 az = _interpret_angle('az', az, az_degrees) 564 r = distance.au 565 p = from_spherical(r, alt, az) 566 p = einsum('ji...,j...->i...', R, p) 567 return Apparent(p) 568 569# For compatibility with my original name for the class. Not an 570# important enough change to warrant a deprecation error for users, so: 571ICRS = ICRF 572 573 574class Geometric(ICRF): 575 """An |xyz| vector between two instantaneous position. 576 577 A geometric position is the difference between the Solar System 578 positions of two bodies at exactly the same instant. It is *not* 579 corrected for the fact that, in real physics, it will take time for 580 light to travel from one position to the other. 581 582 Both the ``.position`` and ``.velocity`` are |xyz| vectors 583 oriented along the axes of the International Celestial Reference 584 System (ICRS), the modern replacement for J2000 coordinates. 585 586 """ 587 def altaz(self, temperature_C=None, pressure_mbar='standard'): 588 """Compute (alt, az, distance) relative to the observer's horizon 589 590 The altitude returned is an :class:`~skyfield.units.Angle` 591 measured in degrees above the horizon, while the azimuth 592 :class:`~skyfield.units.Angle` measures east along the horizon 593 from geographic north (so 0 degrees means north, 90 is east, 180 594 is south, and 270 is west). 595 596 By default, Skyfield does not adjust the altitude for 597 atmospheric refraction. If you want Skyfield to estimate how 598 high the atmosphere might lift the body's image, give the 599 argument ``temperature_C`` either the temperature in degrees 600 centigrade, or the string ``'standard'`` (in which case 10°C is 601 used). 602 603 When calculating refraction, Skyfield uses the observer’s 604 elevation above sea level to estimate the atmospheric pressure. 605 If you want to override that value, simply provide a number 606 through the ``pressure_mbar`` parameter. 607 608 """ 609 return _to_altaz(self, temperature_C, pressure_mbar) 610 611 612class Barycentric(ICRF): 613 """An |xyz| position measured from the Solar System barycenter. 614 615 Skyfield generates a `Barycentric` position measured from the 616 gravitational center of the Solar System whenever you ask a body for 617 its location at a particular time: 618 619 >>> t = ts.utc(2003, 8, 29) 620 >>> mars.at(t) 621 <Barycentric BCRS position and velocity at date t center=0 target=499> 622 623 This class’s ``.position`` and ``.velocity`` are |xyz| vectors in 624 the Barycentric Celestial Reference System (BCRS), the modern 625 replacement for J2000 coordinates measured from the Solar System 626 Barycenter. 627 628 """ 629 _default_center = 0 630 631 def observe(self, body): 632 """Compute the `Astrometric` position of a body from this location. 633 634 To compute the body's astrometric position, it is first asked 635 for its position at the time `t` of this position itself. The 636 distance to the body is then divided by the speed of light to 637 find how long it takes its light to arrive. Finally, the light 638 travel time is subtracted from `t` and the body is asked for a 639 series of increasingly exact positions to learn where it was 640 when it emitted the light that is now reaching this position. 641 642 >>> earth.at(t).observe(mars) 643 <Astrometric ICRS position and velocity at date t center=399 target=499> 644 645 """ 646 p, v, t, light_time = body._observe_from_bcrs(self) 647 astrometric = Astrometric(p, v, t, self.target, body.target) 648 astrometric._ephemeris = self._ephemeris 649 astrometric.center_barycentric = self 650 astrometric.light_time = light_time 651 return astrometric 652 653# TODO: pre-create a Barycentric object representing the SSB, and make 654# it possible for it to observe() a planet. 655 656class Astrometric(ICRF): 657 """An astrometric |xyz| position relative to a particular observer. 658 659 The astrometric position of a body is its position relative to an 660 observer, adjusted for light-time delay. It is the position of the 661 body back when it emitted (or reflected) the light that is now 662 reaching the observer's eye or telescope. Astrometric positions are 663 usually generated in Skyfield by calling the `Barycentric` method 664 `observe()`, which performs the light-time correction. 665 666 Both the ``.position`` and ``.velocity`` are |xyz| vectors 667 oriented along the axes of the ICRF, the modern replacement for the 668 J2000 reference frame. 669 670 It is common to either call ``.radec()`` (with no argument) on an 671 astrometric position to generate an *astrometric place* right 672 ascension and declination with respect to the ICRF axes, or else to 673 call ``.apparent()`` to generate an :class:`Apparent` position. 674 675 """ 676 def apparent(self): 677 """Compute an :class:`Apparent` position for this body. 678 679 This applies two effects to the position that arise from 680 relativity and shift slightly where the other body will appear 681 in the sky: the deflection that the image will experience if its 682 light passes close to large masses in the Solar System, and the 683 aberration of light caused by the observer's own velocity. 684 685 >>> earth.at(t).observe(mars).apparent() 686 <Apparent GCRS position and velocity at date t center=399 target=499> 687 688 These transforms convert the position from the BCRS reference 689 frame of the Solar System barycenter and to the reference frame 690 of the observer. In the specific case of an Earth observer, the 691 output reference frame is the GCRS. 692 693 """ 694 t = self.t 695 target_au = self.position.au.copy() 696 697 cb = self.center_barycentric 698 bcrs_position = cb.position.au 699 bcrs_velocity = cb.velocity.au_per_d 700 observer_gcrs_au = cb._observer_gcrs_au 701 702 # If a single observer position (3,) is observing an array of 703 # targets (3,n), then deflection and aberration will complain 704 # that "operands could not be broadcast together" unless we give 705 # the observer another dimension too. 706 if len(bcrs_position.shape) < len(target_au.shape): 707 shape = bcrs_position.shape + (1,) 708 bcrs_position = bcrs_position.reshape(shape) 709 bcrs_velocity = bcrs_velocity.reshape(shape) 710 if observer_gcrs_au is not None: 711 observer_gcrs_au = observer_gcrs_au.reshape(shape) 712 713 if observer_gcrs_au is None: 714 include_earth_deflection = array((False,)) 715 else: 716 limb_angle, nadir_angle = compute_limb_angle( 717 target_au, observer_gcrs_au) 718 include_earth_deflection = nadir_angle >= 0.8 719 720 add_deflection(target_au, bcrs_position, 721 self._ephemeris, t, include_earth_deflection) 722 723 add_aberration(target_au, bcrs_velocity, self.light_time) 724 725 v = self.velocity.au_per_d 726 if v is not None: 727 pass # TODO: how to apply aberration and deflection to velocity? 728 729 apparent = Apparent(target_au, v, t, self.center, self.target) 730 apparent.center_barycentric = self.center_barycentric 731 apparent._observer_gcrs_au = observer_gcrs_au 732 return apparent 733 734class Apparent(ICRF): 735 """An apparent |xyz| position relative to a particular observer. 736 737 This class’s vectors provide the position and velocity of a body 738 relative to an observer, adjusted to predict where the body’s image 739 will really appear (hence "apparent") in the sky: 740 741 * Light-time delay, as already present in an `Astrometric` position. 742 743 * Deflection: gravity bends light, and thus the image of a distant 744 object, as the light passes massive objects like Jupiter, Saturn, 745 and the Sun. For an observer on the Earth’s surface or in Earth 746 orbit, the slight deflection by the gravity of the Earth itself is 747 also included. 748 749 * Aberration: incoming light arrives slanted because of the 750 observer's motion through space. 751 752 These positions are usually produced in Skyfield by calling the 753 `apparent()` method of an `Astrometric` object. 754 755 Both the ``.position`` and ``.velocity`` are |xyz| vectors 756 oriented along the axes of the ICRF, the modern replacement for the 757 J2000 reference frame. If the observer is at the geocenter, they 758 are more specifically GCRS coordinates. Two common coordinates that 759 this vector can generate are: 760 761 * *Proper place:* call ``.radec()`` without arguments to compute 762 right ascension and declination with respect to the fixed axes of 763 the ICRF. 764 765 * *Apparent place,* the most popular option: call ``.radec('date')`` 766 to generate right ascension and declination with respect to the 767 equator and equinox of date. 768 769 """ 770 def altaz(self, temperature_C=None, pressure_mbar='standard'): 771 """Compute (alt, az, distance) relative to the observer's horizon 772 773 The altitude returned is an :class:`~skyfield.units.Angle` 774 measured in degrees above the horizon, while the azimuth 775 :class:`~skyfield.units.Angle` measures east along the horizon 776 from geographic north (so 0 degrees means north, 90 is east, 180 777 is south, and 270 is west). 778 779 By default, Skyfield does not adjust the altitude for 780 atmospheric refraction. If you want Skyfield to estimate how 781 high the atmosphere might lift the body's image, give the 782 argument ``temperature_C`` either the temperature in degrees 783 centigrade, or the string ``'standard'`` (in which case 10°C is 784 used). 785 786 When calculating refraction, Skyfield uses the observer’s 787 elevation above sea level to estimate the atmospheric pressure. 788 If you want to override that value, simply provide a number 789 through the ``pressure_mbar`` parameter. 790 791 """ 792 return _to_altaz(self, temperature_C, pressure_mbar) 793 794 795class Geocentric(ICRF): 796 """An |xyz| position measured from the center of the Earth. 797 798 A geocentric position is the difference between the position of the 799 Earth at a given instant and the position of a target body at the 800 same instant, without accounting for light-travel time or the effect 801 of relativity on the light itself. 802 803 Its ``.position`` and ``.velocity`` vectors have |xyz| axes that 804 are those of the Geocentric Celestial Reference System (GCRS), an 805 inertial system that is an update to J2000 and that does not rotate 806 with the Earth itself. 807 808 """ 809 _default_center = 399 810 811 def itrf_xyz(self): 812 """Deprecated; instead, call ``.frame_xyz(itrs)``. \ 813 See `reference_frames`.""" 814 return self.frame_xyz(framelib.itrs) 815 816 def subpoint(self): 817 """Deprecated; instead, call either ``iers2010.subpoint(pos)`` or \ 818 ``wgs84.subpoint(pos)``.""" 819 from .toposlib import iers2010 820 return iers2010.subpoint(self) 821 822def _to_altaz(position, temperature_C, pressure_mbar): 823 """Compute (alt, az, distance) relative to the observer's horizon.""" 824 cb = position.center_barycentric 825 if cb is not None: 826 R = cb._altaz_rotation 827 else: 828 rotation_at = getattr(position.center, 'rotation_at') 829 if rotation_at is not None: 830 R = rotation_at(position.t) 831 else: 832 raise ValueError(_altaz_message) 833 834 position_au = mxv(R, position.position.au) 835 r_au, alt, az = to_spherical(position_au) 836 837 if temperature_C is None: 838 alt = Angle(radians=alt) 839 else: 840 refract = getattr(position.center, 'refract', None) 841 if refract is None: 842 raise ValueError(_altaz_message) 843 alt = position.center.refract( 844 alt * RAD2DEG, temperature_C, pressure_mbar, 845 ) 846 847 return alt, Angle(radians=az), Distance(r_au) 848 849_altaz_message = ( 850 'to compute an altazimuth position, you must observe from a' 851 ' specific Earth location or from a position on another body' 852 ' loaded from a set of planetary constants' 853) 854 855class _Fake(ICRF): # support for deprecated frame rotation methods above 856 def __init__(self, original, epoch): 857 self.position = original.position 858 if isinstance(epoch, Time): 859 self.t = epoch 860 elif isinstance(epoch, float): 861 self.t = Time(None, tt=epoch) 862 elif epoch == 'date': 863 self.t = original.t 864 else: 865 raise ValueError('the epoch= must be a Time object,' 866 ' a floating point Terrestrial Time (TT),' 867 ' or the string "date" for epoch-of-date') 868 869def ITRF_to_GCRS(t, rITRF): # Deprecated; for compatibility with old versions. 870 return mxv(_T(framelib.itrs.rotation_at(t)), rITRF) 871 872def ITRF_to_GCRS2(t, rITRF, vITRF, _high_accuracy=False): 873 position = array(rITRF) 874 velocity = array(vITRF) 875 876 # TODO: This is expensive, and should be extensively trimmed to only 877 # include the most important terms underlying GAST. But it improves 878 # the precision by something like 1e5 times when compared to using 879 # the round number skyfield.constants.ANGVEL! 880 # 881 # See the test `test_velocity_in_ITRF_to_GCRS2()`. 882 # 883 if _high_accuracy: 884 _one_second = 1.0 / DAY_S 885 t_later = t.ts.tt_jd(t.whole, t.tt_fraction + _one_second) 886 angvel = (t_later.gast - t.gast) / 24.0 * tau 887 else: 888 angvel = ANGVEL 889 890 spin = rot_z(t.gast / 24.0 * tau) 891 R = mxm(t.MT, spin) 892 893 z = 0.0 * angvel 894 895 V = array(( 896 (z,-DAY_S * angvel,z), 897 (DAY_S * angvel,z,z), 898 (z,z,z), 899 )) 900 901 velocity = velocity + mxv(V, position) 902 position = mxv(R, position) 903 velocity = mxv(R, velocity) 904 905 return position, velocity 906