Updates and correction for v0.3

setup.py

@@ -29,7 +29,7 @@
         'scipy>=1.7.1',
         'requests',
         'tqdm>=4.64',
-        'shapely>=1.8.2',
+        'shapely>=2.0.1',
     ],
     python_requires=">=3.7, <4",
     classifiers=[

sora/body/core.py

@@ -169,7 +169,7 @@ def __from_sbdb(self, name):
         pp = sbdb['phys_par']  # get the physical parameters (pp) of the sbdb
 
         if 'extent' in pp:
-            extent = np.array(pp['extent'].split('x'), dtype=np.float)/2
+            extent = np.array(pp['extent'].split('x'), dtype=float)/2
             self.shape = extent
         self.albedo = PhysicalData('Albedo', pp.get('albedo'), pp.get('albedo_sig'), pp.get('albedo_ref'), pp.get('albedo_note'))
         self.H = PhysicalData('Absolute Magnitude', pp.get('H'), pp.get('H_sig'), pp.get('H_ref'), pp.get('H_note'), unit=u.mag)
@@ -186,8 +186,9 @@ def __from_sbdb(self, name):
         self.UB = PhysicalData('U-B color', pp.get('UB'), pp.get('UB_sig'), pp.get('UB_ref'), pp.get('UB_note'))
         if 'pole' in pp:
             self.pole = SkyCoord(pp['pole'].replace('/', ' '), unit=('deg', 'deg'))
-            self.pole.ra.uncertainty = Longitude(pp['pole_sig'].split('/')[0], unit=u.deg)
-            self.pole.dec.uncertainty = Latitude(pp['pole_sig'].split('/')[1], unit=u.deg)
+            pole_err = pp['pole_sig'].split('/')
+            self.pole.ra.uncertainty = Longitude(pole_err[0], unit=u.deg)
+            self.pole.dec.uncertainty = Latitude(pole_err[0] if len(pole_err) == 1 else pole_err[1], unit=u.deg)
             self.pole.reference = pp['pole_ref'] or ""
             self.pole.notes = pp['pole_note'] or ""
         else:

sora/body/frame/iau_report.py

@@ -24,11 +24,11 @@
     '10': ["286.13 63.87", 0, 0, 84.176, 14.1844, True],  # Sun
     # planets
     '199': ["281.0103 61.4155", 0.0328, 0.0049, 329.5988, 6.1385108, False],  # Mercury
-    '299': ["272.76 67.16", 0, 0, 160.2, -1.4813688, False],  # Venus
-    '499': ["317.269202 54.432516", -0.10927547, -0.05827105, 176.049863, 350.891982443297, False],  # Marte
+    '299': ["272.76 67.16", 0, 0, 160.2, -1.4813688, True],  # Venus
+    '499': ["317.269202 54.432516", -0.10927547, -0.05827105, 176.049863, 350.891982443297, False],  # Mars
     '599': ["268.056595 64.495303", -0.006499, 0.002413, 284.95, 870.536, False],  # Jupiter
     '699': ["40.589 83.537", -0.036, -0.004, 38.9, 810.7939024, False],  # Saturn
-    '799': ["257.311 -15.175", 0, 0, 203.81, -501.1600928, False],  # Uranus
+    '799': ["257.311 -15.175", 0, 0, 203.81, -501.1600928, True],  # Uranus
     '899': ["299.36 46.46", 0, 0, 249.978, 541.1397757, False],  # Neptune
     # satellites
     '401': ["317.67071657 52.88627266", -0.10844326, -0.06134706, 34.9964842535, 1128.84475928, False],  # Phobos
@@ -53,12 +53,12 @@
     '611': ["40.58 83.52", -0.036, -0.004, 293.87, 518.4907239, False],  # Epimetheus
     '612': ["40.85 83.34", -0.036, -0.004, 245.12, 131.6174056, False],  # Helene
     '613': ["50.51 84.06", -0.036, -0.004, 56.88, 190.6979332, False],  # Telesto
-    '614': ["36.41 85.04", -0.036, -0.004, 153.51, 190.6142373, False],  # Calypso
+    '614': ["36.41 85.04", -0.036, -0.004, 153.51, 190.6742373, False],  # Calypso
     '615': ["40.58 83.53", -0.036, -0.004, 137.88, 598.306, False],  # Atlas
     '616': ["40.58 83.53", -0.036, -0.004, 296.14, 587.289, False],  # Prometheus
     '617': ["40.58 83.53", -0.036, -0.004, 162.92, 572.7891, False],  # Pandora
     '618': ["40.6 83.5", 0.036, 0.004, 48.8, 626.044, False],  # Pan
-    '701': ["257.43 -15.10", 0, 0, 156.22, -142.8956681, True],  # Ariel
+    '701': ["257.43 -15.10", 0, 0, 156.22, -142.8356681, True],  # Ariel
     '702': ["257.43 -15.10", 0, 0, 108.05, -86.8688923, True],  # Umbriel
     '703': ["257.43 -15.10", 0, 0, 77.74, -41.3514316, True],  # Titania
     '704': ["257.43 -15.10", 0, 0, 6.77, -26.7394932, True],  # Oberon
@@ -72,28 +72,39 @@
     '712': ["257.31 -15.18", 0, 0, 25.03, -701.4865870, True],  # Portia
     '713': ["257.31 -15.18", 0, 0, 314.90, -644.6311260, True],  # Rosalind
     '714': ["257.31 -15.18", 0, 0, 297.46, -577.3628170, True],  # Belinda
-    '715': ["257.31 -15.18", 0, 0, 91.24, -472.540690, True],  # Puck
-    '801': ["299.36 41.17", 0, 0, 296.57, -61.2572637, False],  # Triton
+    '715': ["257.31 -15.18", 0, 0, 91.24, -472.5450690, True],  # Puck
+    '801': ["299.36 41.17", 0, 0, 296.57, -61.2572637, True],  # Triton
     '803': ["299.36 43.36", 0, 0, 254.06, 1222.8441209, False],  # Naiad
     '804': ["299.36 43.45", 0, 0, 102.06, 1155.7555612, False],  # Thalassa
     '805': ["299.36 43.45", 0, 0, 306.51, 1075.7341562, False],  # Despina
     '806': ["299.36 43.43", 0, 0, 258.09, 839.6597686, False],  # Galatea
     '807': ["299.36 43.41", 0, 0, 179.41, 649.0534470, False],  # Larissa
     '808': ["299.27 42.91", 0, 0, 93.38, 320.7654228, False],  # Proteus
     # asteroids
-    '2000001': ["291.418 66.764", 0, 0, 170.65, 952.1532, False],  # Ceres
-    '2000002': ["33 -3", 0, 0, 38, 1105.8036, False],  # Pallas
-    '2000004': ["309.031 42.235", 0, 0, 285.39, 1617.3329428, False],  # Vesta
-    '2000021': ["52 12", 0, 0, 94, 1057.7515, False],  # Lutetia
-    '2000052': ["257 12", 0, 0, 55, 1534.6472187, False],  # Europa
-    '2000243': ["168.76 -87.12", 0, 0, 274.05, 1864.628007, False],  # Ida
-    '2000433': ["11.35 17.22", 0, 0, 326.07, 1639.38864745, False],  # Eros
-    '2000511': ["297 5", 0, 0, 268.1, 1684.4193549, False],  # Davida
-    '2000951': ["9.47 26.7", 0, 0, 83.67, 1226.911485, False],  # Gaspra
-    '2002867': ["91 -62", 0, 0, 321.76, 1428.09917, False],  # Steins
-    '2025143': ["90.53 -66.3", 0, 0, 0, 712.143, False],  # Itokawa
-    '901': ["132.993 -6.163", 0, 0, 122.695, 56.3625225, False],  # Charon
-    '999': ["132.993 -6.163", 0, 0, 302.695, 56.3625225, False],  # Pluto
+    '2000001': ["291.418 66.764", 0, 0, 170.65, 952.1532, True],  # Ceres
+    '20000001': ["291.418 66.764", 0, 0, 170.65, 952.1532, True],  # Ceres
+    '2000002': ["33 -3", 0, 0, 38, 1105.8036, True],  # Pallas
+    '20000002': ["33 -3", 0, 0, 38, 1105.8036, True],  # Pallas
+    '2000004': ["309.031 42.235", 0, 0, 285.39, 1617.3329428, True],  # Vesta
+    '20000004': ["309.031 42.235", 0, 0, 285.39, 1617.3329428, True],  # Vesta
+    '2000021': ["52 12", 0, 0, 94, 1057.7515, True],  # Lutetia
+    '20000021': ["52 12", 0, 0, 94, 1057.7515, True],  # Lutetia
+    '2000052': ["257 12", 0, 0, 55, 1534.6472187, True],  # Europa
+    '20000052': ["257 12", 0, 0, 55, 1534.6472187, True],  # Europa
+    '2000243': ["168.76 -87.12", 0, 0, 274.05, 1864.628007, True],  # Ida, error in delta0
+    '20000243': ["168.76 -87.12", 0, 0, 274.05, 1864.628007, True],  # Ida
+    '2000433': ["11.35 17.22", 0, 0, 326.07, 1639.38864745, True],  # Eros
+    '20000433': ["11.35 17.22", 0, 0, 326.07, 1639.38864745, True],  # Eros
+    '2000511': ["297 5", 0, 0, 268.1, 1684.4193549, True],  # Davida
+    '20000511': ["297 5", 0, 0, 268.1, 1684.4193549, True],  # Davida
+    '2000951': ["9.47 26.7", 0, 0, 83.67, 1226.911485, True],  # Gaspra
+    '20000951': ["9.47 26.7", 0, 0, 83.67, 1226.911485, True],  # Gaspra
+    '2002867': ["91 -62", 0, 0, 321.76, 1428.09917, True],  # Steins
+    '20002867': ["91 -62", 0, 0, 321.76, 1428.09917, True],  # Steins
+    '2025143': ["90.53 -66.3", 0, 0, 0, 712.143, False],  # Itokawa, error in thermal inertia
+    '20025143': ["90.53 -66.3", 0, 0, 0, 712.143, False],  # Itokawa
+    '901': ["132.993 -6.163", 0, 0, 122.695, 56.3625225, True],  # Charon
+    '999': ["132.993 -6.163", 0, 0, 302.695, 56.3625225, True],  # Pluto
 }
 
 # Include precession terms for sinusoidal functions
@@ -238,8 +249,8 @@
     "601": {
         "extra_alpha": [[13.56, 177.40, -36505.5]],
         "extra_delta": [[-1.53, 177.40, -36505.5]],
-        "extra_w": [[13.48, 177.40, -36505.5],
-                    [44.85, 316.45, 506.2]]
+        "extra_w": [[-13.48, 177.40, -36505.5],
+                    [-44.85, 316.45, 506.2]]
     },
     "603": {
         "extra_alpha": [[9.66, 300.00, -7225.9]],
@@ -337,7 +348,7 @@
     },
     "713": {
         "extra_alpha": [[-0.29, 157.36, 16652.76]],
-        "extra_delta": [[0.28], 157.36, 16652.76],
+        "extra_delta": [[0.28, 157.36, 16652.76]],
         "extra_w": [[-0.08, 157.36, 16652.76]]
     },
     "714": {

sora/body/meta.py

@@ -89,6 +89,8 @@ def uncertainty(self):
     def uncertainty(self, value):
         unit = self.unit
         given_unit = unit
+        if '%' in str(value):
+            value = float(str(value).replace('%', ''))*self.value/100
         if isinstance(value, u.core.CompositeUnit):
             unit = u.dimensionless_unscaled
             for i in np.arange(len(value.bases)):
@@ -397,7 +399,7 @@ def shape(self, value):
         elif isinstance(value, str):
             self._shape = Shape3D(value)
         else:
-            value = np.array(value, ndmin=1, dtype=np.float)
+            value = np.array(value, ndmin=1, dtype=float)
             if len(value) <= 3:
                 self._shape = Ellipsoid(*value)
             else:

sora/body/shape/limb.py

@@ -8,11 +8,15 @@
 zero = geometry.Point(0, 0)
 
 
-class Limb(geometry.LineString):
+class Limb:
+    __doc__ = geometry.LineString.__doc__
+
+    def __init__(self, *args, **kwargs):
+        self._contour = geometry.LineString(*args, **kwargs)
 
     @property
     def xy(self):
-        return np.array(super(Limb, self).xy)
+        return np.array(self._contour.xy)
 
     def plot(self, center_f=0, center_g=0, scale=1, ax=None, **kwargs):
         """Plots the limb on the tangent plane
@@ -44,7 +48,7 @@ def plot(self, center_f=0, center_g=0, scale=1, ax=None, **kwargs):
     @property
     def maxdist(self):
         """Computes the maximum distance of the limb from the origin"""
-        return self.hausdorff_distance(zero)
+        return self._contour.hausdorff_distance(zero)
 
     def radial_residual_to(self, fg):
         """Calculates radial residuals from points.
@@ -70,10 +74,16 @@ def radial_residual_to(self, fg):
         """
         points = geometry.MultiPoint(fg)
         endlines = (fg.T * 1.1 * self.maxdist / np.linalg.norm(fg, axis=-1)).T
-        vals = [point.distance(self.intersection(geometry.LineString([zero, endline]))) for point, endline in
+        vals = [point.distance(self._contour.intersection(geometry.LineString([zero, endline]))) for point, endline in
                 zip(points.geoms, endlines)]
         return np.array(vals)
 
+    def __str__(self):
+        return self._contour.__str__()
+
+    def __repr__(self):
+        return self._contour.__repr__()
+
 
 def limb_radial_residual(limb, fg, center_f=0, center_g=0, scale=1, position_angle=0):
     """

sora/body/shape/utils.py

@@ -30,9 +30,9 @@ def read_obj_file(filename):
     faces = []
     for line in lines:
         if line.startswith('v '):
-            vertices.append(line.strip().split(' ')[1:4])
+            vertices.append(line.strip().split()[1:4])
         elif line.startswith('f '):
-            values = line.strip().split(' ')[1:]
+            values = line.strip().split()[1:]
             faces.append([i.split('/')[0] for i in values])
 
     vertices = np.array(vertices, dtype='float')

sora/body/utils.py

@@ -31,7 +31,7 @@ def search_sbdb(name):
     from . import values
 
     print('Obtaining data for {} from SBDB'.format(name))
-    sbdb = SBDB.query(name, full_precision=True, solution_epoch=True, validity=True, phys=True, discovery=True)
+    sbdb = SBDB.query(name, full_precision=True, solution_epoch=True, validity=True, phys=True, discovery=True, cache=False)
     if 'message' in sbdb:
         if sbdb['message'] == values.not_found_message:
             raise ValueError(values.not_found_message + " on SBDB")

sora/occultation/core.py

@@ -378,7 +378,7 @@ def check_velocities(self):
                       format(normal_vel))
 
     @deprecated_alias(log='verbose')  # remove this line in v1.0
-    def new_astrometric_position(self, time=None, offset=None, error=None, verbose=True, observer=None):
+    def new_astrometric_position(self, time=None, offset=None, error=None, verbose=True, observer=None, sun_ld=False):
         """Calculates the new astrometric position for the object given fitted parameters.
 
         Parameters
@@ -418,6 +418,10 @@ def new_astrometric_position(self, time=None, offset=None, error=None, verbose=T
         observer : `str`, `sora.Observer`, `sora.Spacecraft`
             IAU code of the observer (must be present in given list of kernels),
             a SORA observer object or a string: ['geocenter', 'barycenter']
+
+        sun_ld : `bool`
+            If true, it computes the differential light deflection caused by the Sun
+            in the starlight before being occulted.
         """
         from astropy.coordinates import SkyCoord, SkyOffsetFrame
 
@@ -481,6 +485,17 @@ def new_astrometric_position(self, time=None, offset=None, error=None, verbose=T
         if e_dist:
             e_off_ra = np.arctan2(e_off_ra, distance)
             e_off_dec = np.arctan2(e_off_dec, distance)
+
+        if sun_ld:
+            from .utils import calc_sun_dif_ld
+
+            pos_star = self.star.get_position(time=time, observer=observer)
+            ndra, nddec = calc_sun_dif_ld(body_coord=coord, star_coord=pos_star, time=time, observer=observer)
+            off_ra -= ndra
+            off_dec -= nddec
+            print('Differential Sun Light Deflection included (mas): DRA={:.3f}, DDEC={:.3f}'.format(
+                  ndra.to(u.mas).value, nddec.to(u.mas).value))
+
         new_pos = SkyCoord(lon=off_ra, lat=off_dec, frame=coord_frame)
         new_pos = new_pos.icrs
 

sora/occultation/utils.py

@@ -1,5 +1,10 @@
-import astropy.units as u
 import numpy as np
+import astropy.constants as const
+import astropy.units as u
+from astropy.coordinates import SkyOffsetFrame
+
+from sora.body import Body
+from sora.ephem import EphemHorizons
 
 __all__ = ['filter_negative_chord', 'positionv', 'add_arrow']
 
@@ -222,4 +227,65 @@ def add_arrow(line, position=None, direction='right', size=15, color=None):
                        xy=(xdata[end_ind], ydata[end_ind]),
                        arrowprops=dict(arrowstyle="->", color=color),
                        size=size
-                       )
+                       )
+
+
+sun = Body(name='Sun', spkid='10', ephem=EphemHorizons(name='Sun', spkid=10, id_type='majorbody'),
+           GM=const.G*(1*u.M_sun), database=None)
+
+
+def calc_sun_dif_ld(body_coord, star_coord, time, observer="geocenter"):
+    """Computes differential light deflection of star and body caused by the Sun.
+
+    Parameters
+    ----------
+    body_coord : `astropy.coordinates.SkyCoord`
+        The astrometric coordinate of the body at given time
+
+    star_coord : `astropy.coordinates.SkyCoord`
+        The astrometric coordinate of the star at given time
+
+    time : `astropy.time.Time`
+        The time which to compute the light deflection.
+        Use to compute the position of the Sun
+
+    observer : `sora.observer.Observer`
+        The observer to compute the position of the Sun.
+
+    Returns
+    -------
+    ld_ra, ld_dec : `astropy.coordinates.Angle`
+        The offsets in right ascension and declination of the object relative to the star,
+        considering the light deflection caused by the Sun.
+
+    """
+    import erfa
+    from astropy.coordinates import SkyCoord
+
+    pos_sun = sun.get_position(time=time, observer=observer)
+    bm = 1  # Sun Mass in Solar Masses
+
+    e = -(pos_sun.cartesian.xyz / pos_sun.cartesian.norm()).value  # unitary vector Sun -> observer
+    em = pos_sun.distance.to(u.AU).value  # distance Sun -> observer in AU
+
+    pbody = (body_coord.cartesian.xyz / body_coord.cartesian.norm()).value  # unitary vector observer -> body
+    pbs = (body_coord.cartesian - pos_sun.cartesian)  # vector Sun -> body
+    q = (pbs.xyz / pbs.norm()).value  # unitary vector Sun -> body
+
+    # computes the gravitational light deflection, return new normalized vector observer -> body
+    ds_body = erfa.ld(bm, pbody, q, e, em, 0)
+
+    os_body = SkyCoord(*ds_body * body_coord.cartesian.norm(), representation_type='cartesian')
+    dra_b, ddec_b = body_coord.spherical_offsets_to(os_body)  # computes offset between corrected and astrometric
+
+    pstar = (star_coord.cartesian.xyz / star_coord.cartesian.norm()).value  # unitary vector observer -> star
+    pbs = (star_coord.cartesian - pos_sun.cartesian)  # vector Sun -> star
+    q = (pbs.xyz / pbs.norm()).value  # unitary vector Sun -> star
+
+    # computes the gravitational light deflection, return new normalized vector observer -> body
+    ds_star = erfa.ld(bm, pstar, q, e, em, 0)
+
+    os_star = SkyCoord(*ds_star * star_coord.cartesian.norm(), representation_type='cartesian')
+    dra_s, ddec_s = star_coord.spherical_offsets_to(os_star)  # computes offset between corrected and astrometric
+
+    return dra_b - dra_s, ddec_b - ddec_s

sora/prediction/occmap.py

@@ -617,12 +617,18 @@ def plot_occ_map(name, radius, coord, time, ca, pa, vel, dist, mag=0, longi=0, *
         lon1, lat1 = xy2latlon(ax2.to(u.m).value, by2.to(u.m).value, centers.lon.value, centers.lat.value, datas1)
         j = np.where(lon1 < 1e+30)
         axf.plot(lon1[j], lat1[j], '--', transform=ccrs.Geodetic(), color=ercolor)
+        j = np.where(lon1 > 1e+30)
+        if 'centerproj' not in kwargs:
+            plt.plot(ax2[j].to(u.m).value, by2[j].to(u.m).value, '--', color=ercolor, clip_on=(not centert), zorder=-0.2)
 
         ax3 = ax + errd*np.sin(paplus) + radius*np.sin(paplus)
         by3 = by + errd*np.cos(paplus) + radius*np.cos(paplus)
         lon2, lat2 = xy2latlon(ax3.to(u.m).value, by3.to(u.m).value, centers.lon.value, centers.lat.value, datas1)
         j = np.where(lon2 < 1e+30)
         axf.plot(lon2[j], lat2[j], '--', transform=ccrs.Geodetic(), color=ercolor)
+        j = np.where(lon1 > 1e+30)
+        if 'centerproj' not in kwargs:
+            plt.plot(ax3[j].to(u.m).value, by3[j].to(u.m).value, '--', color=ercolor, clip_on=(not centert), zorder=-0.2)
 
     # plots ring
     if ring is not None:
@@ -632,12 +638,18 @@ def plot_occ_map(name, radius, coord, time, ca, pa, vel, dist, mag=0, longi=0, *
         lon1, lat1 = xy2latlon(ax2.to(u.m).value, by2.to(u.m).value, centers.lon.value, centers.lat.value, datas1)
         j = np.where(lon1 < 1e+30)
         axf.plot(lon1[j], lat1[j], '--', transform=ccrs.Geodetic(), color=rncolor)
+        j = np.where(lon1 > 1e+30)
+        if 'centerproj' not in kwargs:
+            plt.plot(ax2[j].to(u.m).value, by2[j].to(u.m).value, '--', color=rncolor, clip_on=(not centert), zorder=-0.2)
 
         ax3 = ax + rng*np.sin(paplus)
         by3 = by + rng*np.cos(paplus)
         lon2, lat2 = xy2latlon(ax3.to(u.m).value, by3.to(u.m).value, centers.lon.value, centers.lat.value, datas1)
         j = np.where(lon2 < 1e+30)
         axf.plot(lon2[j], lat2[j], '--', transform=ccrs.Geodetic(), color=rncolor)
+        j = np.where(lon1 > 1e+30)
+        if 'centerproj' not in kwargs:
+            plt.plot(ax3[j].to(u.m).value, by3[j].to(u.m).value, '--', color=rncolor, clip_on=(not centert), zorder=-0.2)
 
     # plots atm
     if atm is not None:
@@ -647,12 +659,18 @@ def plot_occ_map(name, radius, coord, time, ca, pa, vel, dist, mag=0, longi=0, *
         lon1, lat1 = xy2latlon(ax2.to(u.m).value, by2.to(u.m).value, centers.lon.value, centers.lat.value, datas1)
         j = np.where(lon1 < 1e+30)
         axf.plot(lon1[j], lat1[j], '--', transform=ccrs.Geodetic(), color=atcolor)
+        j = np.where(lon1 > 1e+30)
+        if 'centerproj' not in kwargs:
+            plt.plot(ax2[j].to(u.m).value, by2[j].to(u.m).value, '--', color=atcolor, clip_on=(not centert), zorder=-0.2)
 
         ax3 = ax + atmo*np.sin(paplus)
         by3 = by + atmo*np.cos(paplus)
         lon2, lat2 = xy2latlon(ax3.to(u.m).value, by3.to(u.m).value, centers.lon.value, centers.lat.value, datas1)
         j = np.where(lon2 < 1e+30)
         axf.plot(lon2[j], lat2[j], '--', transform=ccrs.Geodetic(), color=atcolor)
+        j = np.where(lon1 > 1e+30)
+        if 'centerproj' not in kwargs:
+            plt.plot(ax3[j].to(u.m).value, by3[j].to(u.m).value, '--', color=atcolor, clip_on=(not centert), zorder=-0.2)
 
     # plots center points
     vec = np.arange(0, int(8000/(np.absolute(occs['vel'].value))), cpoints)