1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2""" 3This module contains convenience functions for retrieving solar system 4ephemerides from jplephem. 5""" 6 7from urllib.parse import urlparse 8import os.path 9 10import numpy as np 11import erfa 12 13from .sky_coordinate import SkyCoord 14from astropy.utils.data import download_file 15from astropy.utils.decorators import classproperty, deprecated 16from astropy.utils.state import ScienceState 17from astropy.utils import indent 18from astropy import units as u 19from astropy.constants import c as speed_of_light 20from .representation import CartesianRepresentation, CartesianDifferential 21from .builtin_frames import GCRS, ICRS, ITRS, TETE 22from .builtin_frames.utils import get_jd12 23 24__all__ = ["get_body", "get_moon", "get_body_barycentric", 25 "get_body_barycentric_posvel", "solar_system_ephemeris"] 26 27 28DEFAULT_JPL_EPHEMERIS = 'de430' 29 30"""List of kernel pairs needed to calculate positions of a given object.""" 31BODY_NAME_TO_KERNEL_SPEC = { 32 'sun': [(0, 10)], 33 'mercury': [(0, 1), (1, 199)], 34 'venus': [(0, 2), (2, 299)], 35 'earth-moon-barycenter': [(0, 3)], 36 'earth': [(0, 3), (3, 399)], 37 'moon': [(0, 3), (3, 301)], 38 'mars': [(0, 4)], 39 'jupiter': [(0, 5)], 40 'saturn': [(0, 6)], 41 'uranus': [(0, 7)], 42 'neptune': [(0, 8)], 43 'pluto': [(0, 9)], 44} 45 46"""Indices to the plan94 routine for the given object.""" 47PLAN94_BODY_NAME_TO_PLANET_INDEX = { 48 'mercury': 1, 49 'venus': 2, 50 'earth-moon-barycenter': 3, 51 'mars': 4, 52 'jupiter': 5, 53 'saturn': 6, 54 'uranus': 7, 55 'neptune': 8, 56} 57 58_EPHEMERIS_NOTE = """ 59You can either give an explicit ephemeris or use a default, which is normally 60a built-in ephemeris that does not require ephemeris files. To change 61the default to be the JPL ephemeris:: 62 63 >>> from astropy.coordinates import solar_system_ephemeris 64 >>> solar_system_ephemeris.set('jpl') # doctest: +SKIP 65 66Use of any JPL ephemeris requires the jplephem package 67(https://pypi.org/project/jplephem/). 68If needed, the ephemeris file will be downloaded (and cached). 69 70One can check which bodies are covered by a given ephemeris using:: 71 72 >>> solar_system_ephemeris.bodies 73 ('earth', 'sun', 'moon', 'mercury', 'venus', 'earth-moon-barycenter', 'mars', 'jupiter', 'saturn', 'uranus', 'neptune') 74"""[1:-1] 75 76 77class solar_system_ephemeris(ScienceState): 78 """Default ephemerides for calculating positions of Solar-System bodies. 79 80 This can be one of the following:: 81 82 - 'builtin': polynomial approximations to the orbital elements. 83 - 'de430', 'de432s', 'de440', 'de440s': short-cuts for recent JPL dynamical models. 84 - 'jpl': Alias for the default JPL ephemeris (currently, 'de430'). 85 - URL: (str) The url to a SPK ephemeris in SPICE binary (.bsp) format. 86 - PATH: (str) File path to a SPK ephemeris in SPICE binary (.bsp) format. 87 - `None`: Ensure an Exception is raised without an explicit ephemeris. 88 89 The default is 'builtin', which uses the ``epv00`` and ``plan94`` 90 routines from the ``erfa`` implementation of the Standards Of Fundamental 91 Astronomy library. 92 93 Notes 94 ----- 95 Any file required will be downloaded (and cached) when the state is set. 96 The default Satellite Planet Kernel (SPK) file from NASA JPL (de430) is 97 ~120MB, and covers years ~1550-2650 CE [1]_. The smaller de432s file is 98 ~10MB, and covers years 1950-2050 [2]_ (and similarly for the newer de440 99 and de440s). Older versions of the JPL ephemerides (such as the widely 100 used de200) can be used via their URL [3]_. 101 102 .. [1] https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/aareadme_de430-de431.txt 103 .. [2] https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/aareadme_de432s.txt 104 .. [3] https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/ 105 """ 106 _value = 'builtin' 107 _kernel = None 108 109 @classmethod 110 def validate(cls, value): 111 # make no changes if value is None 112 if value is None: 113 return cls._value 114 # Set up Kernel; if the file is not in cache, this will download it. 115 cls.get_kernel(value) 116 return value 117 118 @classmethod 119 def get_kernel(cls, value): 120 # ScienceState only ensures the `_value` attribute is up to date, 121 # so we need to be sure any kernel returned is consistent. 122 if cls._kernel is None or cls._kernel.origin != value: 123 if cls._kernel is not None: 124 cls._kernel.daf.file.close() 125 cls._kernel = None 126 kernel = _get_kernel(value) 127 if kernel is not None: 128 kernel.origin = value 129 cls._kernel = kernel 130 return cls._kernel 131 132 @classproperty 133 def kernel(cls): 134 return cls.get_kernel(cls._value) 135 136 @classproperty 137 def bodies(cls): 138 if cls._value is None: 139 return None 140 if cls._value.lower() == 'builtin': 141 return (('earth', 'sun', 'moon') + 142 tuple(PLAN94_BODY_NAME_TO_PLANET_INDEX.keys())) 143 else: 144 return tuple(BODY_NAME_TO_KERNEL_SPEC.keys()) 145 146 147def _get_kernel(value): 148 """ 149 Try importing jplephem, download/retrieve from cache the Satellite Planet 150 Kernel corresponding to the given ephemeris. 151 """ 152 if value is None or value.lower() == 'builtin': 153 return None 154 155 try: 156 from jplephem.spk import SPK 157 except ImportError: 158 raise ImportError("Solar system JPL ephemeris calculations require " 159 "the jplephem package " 160 "(https://pypi.org/project/jplephem/)") 161 162 if value.lower() == 'jpl': 163 value = DEFAULT_JPL_EPHEMERIS 164 165 if value.lower() in ('de430', 'de432s', 'de440', 'de440s'): 166 value = ('https://naif.jpl.nasa.gov/pub/naif/generic_kernels' 167 '/spk/planets/{:s}.bsp'.format(value.lower())) 168 169 elif os.path.isfile(value): 170 return SPK.open(value) 171 172 else: 173 try: 174 urlparse(value) 175 except Exception: 176 raise ValueError('{} was not one of the standard strings and ' 177 'could not be parsed as a file path or URL'.format(value)) 178 179 return SPK.open(download_file(value, cache=True)) 180 181 182def _get_body_barycentric_posvel(body, time, ephemeris=None, 183 get_velocity=True): 184 """Calculate the barycentric position (and velocity) of a solar system body. 185 186 Parameters 187 ---------- 188 body : str or other 189 The solar system body for which to calculate positions. Can also be a 190 kernel specifier (list of 2-tuples) if the ``ephemeris`` is a JPL 191 kernel. 192 time : `~astropy.time.Time` 193 Time of observation. 194 ephemeris : str, optional 195 Ephemeris to use. By default, use the one set with 196 ``astropy.coordinates.solar_system_ephemeris.set`` 197 get_velocity : bool, optional 198 Whether or not to calculate the velocity as well as the position. 199 200 Returns 201 ------- 202 position : `~astropy.coordinates.CartesianRepresentation` or tuple 203 Barycentric (ICRS) position or tuple of position and velocity. 204 205 Notes 206 ----- 207 Whether or not velocities are calculated makes little difference for the 208 built-in ephemerides, but for most JPL ephemeris files, the execution time 209 roughly doubles. 210 """ 211 # If the ephemeris is to be taken from solar_system_ephemeris, or the one 212 # it already contains, use the kernel there. Otherwise, open the ephemeris, 213 # possibly downloading it, but make sure the file is closed at the end. 214 default_kernel = ephemeris is None or ephemeris is solar_system_ephemeris._value 215 kernel = None 216 try: 217 if default_kernel: 218 if solar_system_ephemeris.get() is None: 219 raise ValueError(_EPHEMERIS_NOTE) 220 kernel = solar_system_ephemeris.kernel 221 else: 222 kernel = _get_kernel(ephemeris) 223 224 jd1, jd2 = get_jd12(time, 'tdb') 225 if kernel is None: 226 body = body.lower() 227 earth_pv_helio, earth_pv_bary = erfa.epv00(jd1, jd2) 228 if body == 'earth': 229 body_pv_bary = earth_pv_bary 230 231 elif body == 'moon': 232 # The moon98 documentation notes that it takes TT, but that TDB leads 233 # to errors smaller than the uncertainties in the algorithm. 234 # moon98 returns the astrometric position relative to the Earth. 235 moon_pv_geo = erfa.moon98(jd1, jd2) 236 body_pv_bary = erfa.pvppv(moon_pv_geo, earth_pv_bary) 237 else: 238 sun_pv_bary = erfa.pvmpv(earth_pv_bary, earth_pv_helio) 239 if body == 'sun': 240 body_pv_bary = sun_pv_bary 241 else: 242 try: 243 body_index = PLAN94_BODY_NAME_TO_PLANET_INDEX[body] 244 except KeyError: 245 raise KeyError("{}'s position and velocity cannot be " 246 "calculated with the '{}' ephemeris." 247 .format(body, ephemeris)) 248 body_pv_helio = erfa.plan94(jd1, jd2, body_index) 249 body_pv_bary = erfa.pvppv(body_pv_helio, sun_pv_bary) 250 251 body_pos_bary = CartesianRepresentation( 252 body_pv_bary['p'], unit=u.au, xyz_axis=-1, copy=False) 253 if get_velocity: 254 body_vel_bary = CartesianRepresentation( 255 body_pv_bary['v'], unit=u.au/u.day, xyz_axis=-1, 256 copy=False) 257 258 else: 259 if isinstance(body, str): 260 # Look up kernel chain for JPL ephemeris, based on name 261 try: 262 kernel_spec = BODY_NAME_TO_KERNEL_SPEC[body.lower()] 263 except KeyError: 264 raise KeyError("{}'s position cannot be calculated with " 265 "the {} ephemeris.".format(body, ephemeris)) 266 else: 267 # otherwise, assume the user knows what their doing and intentionally 268 # passed in a kernel chain 269 kernel_spec = body 270 271 # jplephem cannot handle multi-D arrays, so convert to 1D here. 272 jd1_shape = getattr(jd1, 'shape', ()) 273 if len(jd1_shape) > 1: 274 jd1, jd2 = jd1.ravel(), jd2.ravel() 275 # Note that we use the new jd1.shape here to create a 1D result array. 276 # It is reshaped below. 277 body_posvel_bary = np.zeros((2 if get_velocity else 1, 3) + 278 getattr(jd1, 'shape', ())) 279 for pair in kernel_spec: 280 spk = kernel[pair] 281 if spk.data_type == 3: 282 # Type 3 kernels contain both position and velocity. 283 posvel = spk.compute(jd1, jd2) 284 if get_velocity: 285 body_posvel_bary += posvel.reshape(body_posvel_bary.shape) 286 else: 287 body_posvel_bary[0] += posvel[:4] 288 else: 289 # spk.generate first yields the position and then the 290 # derivative. If no velocities are desired, body_posvel_bary 291 # has only one element and thus the loop ends after a single 292 # iteration, avoiding the velocity calculation. 293 for body_p_or_v, p_or_v in zip(body_posvel_bary, 294 spk.generate(jd1, jd2)): 295 body_p_or_v += p_or_v 296 297 body_posvel_bary.shape = body_posvel_bary.shape[:2] + jd1_shape 298 body_pos_bary = CartesianRepresentation(body_posvel_bary[0], 299 unit=u.km, copy=False) 300 if get_velocity: 301 body_vel_bary = CartesianRepresentation(body_posvel_bary[1], 302 unit=u.km/u.day, copy=False) 303 304 return (body_pos_bary, body_vel_bary) if get_velocity else body_pos_bary 305 306 finally: 307 if not default_kernel and kernel is not None: 308 kernel.daf.file.close() 309 310 311def get_body_barycentric_posvel(body, time, ephemeris=None): 312 """Calculate the barycentric position and velocity of a solar system body. 313 314 Parameters 315 ---------- 316 body : str or list of tuple 317 The solar system body for which to calculate positions. Can also be a 318 kernel specifier (list of 2-tuples) if the ``ephemeris`` is a JPL 319 kernel. 320 time : `~astropy.time.Time` 321 Time of observation. 322 ephemeris : str, optional 323 Ephemeris to use. By default, use the one set with 324 ``astropy.coordinates.solar_system_ephemeris.set`` 325 326 Returns 327 ------- 328 position, velocity : tuple of `~astropy.coordinates.CartesianRepresentation` 329 Tuple of barycentric (ICRS) position and velocity. 330 331 See also 332 -------- 333 get_body_barycentric : to calculate position only. 334 This is faster by about a factor two for JPL kernels, but has no 335 speed advantage for the built-in ephemeris. 336 337 Notes 338 ----- 339 {_EPHEMERIS_NOTE} 340 """ 341 return _get_body_barycentric_posvel(body, time, ephemeris) 342 343 344def get_body_barycentric(body, time, ephemeris=None): 345 """Calculate the barycentric position of a solar system body. 346 347 Parameters 348 ---------- 349 body : str or list of tuple 350 The solar system body for which to calculate positions. Can also be a 351 kernel specifier (list of 2-tuples) if the ``ephemeris`` is a JPL 352 kernel. 353 time : `~astropy.time.Time` 354 Time of observation. 355 ephemeris : str, optional 356 Ephemeris to use. By default, use the one set with 357 ``astropy.coordinates.solar_system_ephemeris.set`` 358 359 Returns 360 ------- 361 position : `~astropy.coordinates.CartesianRepresentation` 362 Barycentric (ICRS) position of the body in cartesian coordinates 363 364 See also 365 -------- 366 get_body_barycentric_posvel : to calculate both position and velocity. 367 368 Notes 369 ----- 370 {_EPHEMERIS_NOTE} 371 """ 372 return _get_body_barycentric_posvel(body, time, ephemeris, 373 get_velocity=False) 374 375 376def _get_apparent_body_position(body, time, ephemeris, obsgeoloc=None): 377 """Calculate the apparent position of body ``body`` relative to Earth. 378 379 This corrects for the light-travel time to the object. 380 381 Parameters 382 ---------- 383 body : str or other 384 The solar system body for which to calculate positions. Can also be a 385 kernel specifier (list of 2-tuples) if the ``ephemeris`` is a JPL 386 kernel. 387 time : `~astropy.time.Time` 388 Time of observation. 389 ephemeris : str, optional 390 Ephemeris to use. By default, use the one set with 391 ``~astropy.coordinates.solar_system_ephemeris.set`` 392 obsgeoloc : `~astropy.coordinates.CartesianRepresentation`, optional 393 The GCRS position of the observer 394 395 Returns 396 ------- 397 cartesian_position : `~astropy.coordinates.CartesianRepresentation` 398 Barycentric (ICRS) apparent position of the body in cartesian coordinates 399 400 Notes 401 ----- 402 {_EPHEMERIS_NOTE} 403 """ 404 if ephemeris is None: 405 ephemeris = solar_system_ephemeris.get() 406 407 # Calculate position given approximate light travel time. 408 delta_light_travel_time = 20. * u.s 409 emitted_time = time 410 light_travel_time = 0. * u.s 411 earth_loc = get_body_barycentric('earth', time, ephemeris) 412 if obsgeoloc is not None: 413 earth_loc += obsgeoloc 414 while np.any(np.fabs(delta_light_travel_time) > 1.0e-8*u.s): 415 body_loc = get_body_barycentric(body, emitted_time, ephemeris) 416 earth_distance = (body_loc - earth_loc).norm() 417 delta_light_travel_time = (light_travel_time - 418 earth_distance/speed_of_light) 419 light_travel_time = earth_distance/speed_of_light 420 emitted_time = time - light_travel_time 421 422 return get_body_barycentric(body, emitted_time, ephemeris) 423 424 425def get_body(body, time, location=None, ephemeris=None): 426 """ 427 Get a `~astropy.coordinates.SkyCoord` for a solar system body as observed 428 from a location on Earth in the `~astropy.coordinates.GCRS` reference 429 system. 430 431 Parameters 432 ---------- 433 body : str or list of tuple 434 The solar system body for which to calculate positions. Can also be a 435 kernel specifier (list of 2-tuples) if the ``ephemeris`` is a JPL 436 kernel. 437 time : `~astropy.time.Time` 438 Time of observation. 439 location : `~astropy.coordinates.EarthLocation`, optional 440 Location of observer on the Earth. If not given, will be taken from 441 ``time`` (if not present, a geocentric observer will be assumed). 442 ephemeris : str, optional 443 Ephemeris to use. If not given, use the one set with 444 ``astropy.coordinates.solar_system_ephemeris.set`` (which is 445 set to 'builtin' by default). 446 447 Returns 448 ------- 449 skycoord : `~astropy.coordinates.SkyCoord` 450 GCRS Coordinate for the body 451 452 Notes 453 ----- 454 The coordinate returned is the apparent position, which is the position of 455 the body at time *t* minus the light travel time from the *body* to the 456 observing *location*. 457 458 {_EPHEMERIS_NOTE} 459 """ 460 if location is None: 461 location = time.location 462 463 if location is not None: 464 obsgeoloc, obsgeovel = location.get_gcrs_posvel(time) 465 else: 466 obsgeoloc, obsgeovel = None, None 467 468 cartrep = _get_apparent_body_position(body, time, ephemeris, obsgeoloc) 469 icrs = ICRS(cartrep) 470 gcrs = icrs.transform_to(GCRS(obstime=time, 471 obsgeoloc=obsgeoloc, 472 obsgeovel=obsgeovel)) 473 474 return SkyCoord(gcrs) 475 476 477def get_moon(time, location=None, ephemeris=None): 478 """ 479 Get a `~astropy.coordinates.SkyCoord` for the Earth's Moon as observed 480 from a location on Earth in the `~astropy.coordinates.GCRS` reference 481 system. 482 483 Parameters 484 ---------- 485 time : `~astropy.time.Time` 486 Time of observation 487 location : `~astropy.coordinates.EarthLocation` 488 Location of observer on the Earth. If none is supplied, taken from 489 ``time`` (if not present, a geocentric observer will be assumed). 490 ephemeris : str, optional 491 Ephemeris to use. If not given, use the one set with 492 ``astropy.coordinates.solar_system_ephemeris.set`` (which is 493 set to 'builtin' by default). 494 495 Returns 496 ------- 497 skycoord : `~astropy.coordinates.SkyCoord` 498 GCRS Coordinate for the Moon 499 500 Notes 501 ----- 502 {_EPHEMERIS_NOTE} 503 """ 504 505 return get_body('moon', time, location=location, ephemeris=ephemeris) 506 507 508# Add note about the ephemeris choices to the docstrings of relevant functions. 509# Note: sadly, one cannot use f-strings for docstrings, so we format explicitly. 510for f in [f for f in locals().values() if callable(f) and f.__doc__ is not None 511 and '{_EPHEMERIS_NOTE}' in f.__doc__]: 512 f.__doc__ = f.__doc__.format(_EPHEMERIS_NOTE=indent(_EPHEMERIS_NOTE)[4:]) 513 514 515deprecation_msg = """ 516The use of _apparent_position_in_true_coordinates is deprecated because 517astropy now implements a True Equator True Equinox Frame (TETE), which 518should be used instead. 519""" 520 521 522@deprecated('4.2', deprecation_msg) 523def _apparent_position_in_true_coordinates(skycoord): 524 """ 525 Convert Skycoord in GCRS frame into one in which RA and Dec 526 are defined w.r.t to the true equinox and poles of the Earth 527 """ 528 location = getattr(skycoord, 'location', None) 529 if location is None: 530 gcrs_rep = skycoord.obsgeoloc.with_differentials( 531 {'s': CartesianDifferential.from_cartesian(skycoord.obsgeovel)}) 532 location = (GCRS(gcrs_rep, obstime=skycoord.obstime) 533 .transform_to(ITRS(obstime=skycoord.obstime)) 534 .earth_location) 535 tete_frame = TETE(obstime=skycoord.obstime, location=location) 536 return skycoord.transform_to(tete_frame) 537