1"""Vector functions and their composition.""" 2 3from jplephem.names import target_names as _jpl_code_name_dict 4from numpy import max 5from .constants import C_AUDAY 6from .descriptorlib import reify 7from .errors import DeprecationError 8from .functions import length_of 9from .positionlib import build_position 10from .timelib import Time 11 12class VectorFunction(object): 13 """Given a time, computes a corresponding position.""" 14 15 ephemeris = None 16 17 @reify 18 def vector_name(self): 19 return type(self).__name__ 20 21 @reify 22 def center_name(self): 23 return _jpl_name(self.center) 24 25 @reify 26 def target_name(self): 27 return _jpl_name(self.target) 28 29 def __repr__(self): 30 return '<{0} {1}>'.format(type(self).__name__, str(self)) 31 32 def __str__(self): 33 if self.target is self: 34 return self.target_name 35 return self.arrow_str() 36 37 def arrow_str(self): 38 return '{0} {1} -> {2}'.format( 39 self.vector_name, self.center_name, self.target_name, 40 ) 41 42 def __add__(self, other): 43 if self.target != other.center: 44 if other.target == self.center: 45 self, other = other, self 46 else: 47 raise ValueError( 48 "you can only add two vectors" 49 " if the target where one of the vectors ends" 50 " is the center where the other vector starts" 51 ) 52 53 self_vfs = getattr(self, 'vector_functions', None) or (self,) 54 other_vfs = getattr(other, 'vector_functions', None) or (other,) 55 56 return VectorSum(self.center, other.target, self_vfs + other_vfs) 57 58 def __neg__(self): 59 return ReversedVector(self) 60 61 def __sub__(self, other): 62 if self.center != other.center: 63 raise ValueError( 64 "you can only subtract two vectors" 65 " if they both start at the same center" 66 ) 67 68 self_vfs = getattr(self, 'vector_functions', None) or (self,) 69 other_vfs = getattr(other, 'vector_functions', None) or (other,) 70 other_vfs = tuple(reversed([-vf for vf in other_vfs])) 71 72 return VectorSum(other.target, self.target, other_vfs + self_vfs) 73 74 def at(self, t): 75 """At time ``t``, compute the target's position relative to the center. 76 77 If ``t`` is an array of times, then the returned position object 78 will specify as many positions as there were times. The kind of 79 position returned depends on the value of the ``center`` 80 attribute: 81 82 * Solar System Barycenter: :class:`~skyfield.positionlib.Barycentric` 83 * Center of the Earth: :class:`~skyfield.positionlib.Geocentric` 84 * Difference: :class:`~skyfield.positionlib.Geometric` 85 * Anything else: :class:`~skyfield.positionlib.ICRF` 86 87 """ 88 if not isinstance(t, Time): 89 raise ValueError('please provide the at() method with a Time' 90 ' instance as its argument, instead of the' 91 ' value {0!r}'.format(t)) 92 p, v, gcrs_position, message = self._at(t) 93 center = self.center 94 position = build_position(p, v, t, center, self.target) 95 position._ephemeris = self.ephemeris 96 position._observer_gcrs_au = gcrs_position 97 position.message = message 98 return position 99 100 def _observe_from_bcrs(self, observer): 101 if self.center != 0: 102 raise ValueError('you can only observe() a body whose vector' 103 ' center is the Solar System Barycenter,' 104 ' but this vector has the center {0}' 105 .format(self.center_name)) 106 return _correct_for_light_travel_time(observer, self) 107 108 def geometry_of(self, other): 109 raise DeprecationError( 110"""the geometry_of() method has, alas, been deprecated 111 112This old method has been replaced by an improved interface. If you just 113need your software working again, install Skyfield 0.9.1 for a quick fix: 114 115 pip install skyfield==0.9.1 116 117Or, to update your old code, replace each operation that looks like: 118 119 position = boston.geometry_of(satellite).at(t) 120 121with the vector math that was previously hiding inside the old method: 122 123 position = (satellite - boston).at(t)""") 124 125 def topos(self, latitude=None, longitude=None, latitude_degrees=None, 126 longitude_degrees=None, elevation_m=0.0, x=0.0, y=0.0): 127 raise DeprecationError( 128"""the topos() method has, alas, been deprecated 129 130This old method has been replaced by an improved interface. If you just 131need your software working again, install Skyfield 0.9.1 for a quick fix: 132 133 pip install skyfield==0.9.1 134 135Or, to update your old code, replace each operation that looks like: 136 137 boston = earth.topos(...) 138 139with the vector math that was previously hiding inside the old method: 140 141 from skyfield.api import Topos 142 boston = earth + Topos(...)""") 143 144 def satellite(self, text): 145 raise DeprecationError( 146"""the satellite() method has, alas, been deprecated 147 148This old method has been replaced by an improved interface. If you just 149need your software working again, install Skyfield 0.9.1 for a quick fix: 150 151 pip install skyfield==0.9.1 152 153Or, to update your old code, replace each operation that looks like: 154 155 sat = earth.satellite(tle_text) 156 157with the vector math (and the little bit of text manipulation) that was 158previously hiding inside the old method: 159 160 from skyfield.api import EarthSatellite 161 line1, line2 = tle_text.splitlines()[-2:] 162 sat = earth + EarthSatellite(line1, line2)""") 163 164class ReversedVector(VectorFunction): 165 def __init__(self, vector_function): 166 self.center = vector_function.target 167 self.target = vector_function.center 168 self.vector_function = vector_function 169 170 @reify 171 def vector_name(self): 172 return 'Reversed ' + self.vector_function.vector_name 173 174 @reify 175 def center_name(self): 176 return self.vector_function.target_name 177 178 @reify 179 def target_name(self): 180 return self.vector_function.center_name 181 182 def __neg__(self): 183 return self.vector_function 184 185 def _at(self, t): 186 p, v, _, message = self.vector_function._at(t) 187 return -p, -v, None, message 188 189class VectorSum(VectorFunction): 190 def __init__(self, center, target, vector_functions): 191 self.center = center 192 self.target = target 193 self.vector_functions = vector_functions 194 195 # For now, just grab the first ephemeris we can find. 196 ephemerides = (segment.ephemeris for segment in vector_functions 197 if segment.ephemeris) 198 self.ephemeris = next(ephemerides, None) 199 200 def __str__(self): 201 vector_functions = self.vector_functions 202 lines = [' ' + segment.arrow_str() for segment in vector_functions] 203 return 'Sum of {0} vectors:\n{1}'.format( 204 len(vector_functions), 205 '\n'.join(lines), 206 ) 207 208 def __repr__(self): 209 return '<Vector{0}>'.format(self) 210 211 def _at(self, t): 212 p, v = 0.0, 0.0 213 gcrs_position = None 214 vfs = self.vector_functions 215 for vf in vfs: 216 p2, v2, _, message = vf._at(t) 217 if vf.center == 399: 218 gcrs_position = -p 219 p += p2 220 v += v2 221 if vfs[0].center == 0 and vf.center == 399: 222 gcrs_position = p2 223 return p, v, gcrs_position, message 224 225def _correct_for_light_travel_time(observer, target): 226 """Return a light-time corrected astrometric position and velocity. 227 228 Given an `observer` that is a `Barycentric` position somewhere in 229 the solar system, compute where in the sky they will see the body 230 `target`, by computing the light-time between them and figuring out 231 where `target` was back when the light was leaving it that is now 232 reaching the eyes or instruments of the `observer`. 233 234 """ 235 t = observer.t 236 ts = t.ts 237 whole = t.whole 238 tdb_fraction = t.tdb_fraction 239 240 cposition = observer.position.au 241 cvelocity = observer.velocity.au_per_d 242 243 tposition, tvelocity, gcrs_position, message = target._at(t) 244 245 distance = length_of(tposition - cposition) 246 light_time0 = 0.0 247 for i in range(10): 248 light_time = distance / C_AUDAY 249 delta = light_time - light_time0 250 if max(abs(delta)) < 1e-12: 251 break 252 253 # We assume a light travel time of at most a couple of days. A 254 # longer light travel time would best be split into a whole and 255 # fraction, for adding to the whole and fraction of TDB. 256 t2 = ts.tdb_jd(whole, tdb_fraction - light_time) 257 258 tposition, tvelocity, gcrs_position, message = target._at(t2) 259 distance = length_of(tposition - cposition) 260 light_time0 = light_time 261 else: 262 raise ValueError('light-travel time failed to converge') 263 return tposition - cposition, tvelocity - cvelocity, t, light_time 264 265def _jpl_name(target): 266 if not isinstance(target, int): 267 return type(target).__name__ 268 name = _jpl_code_name_dict.get(target) 269 if name is None: 270 return str(target) 271 return '{0} {1}'.format(target, name) 272