diff --git a/setup.py b/setup.py index 02ac8a2..b3256fc 100644 --- a/setup.py +++ b/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=[ diff --git a/sora/body/core.py b/sora/body/core.py index 90c2c34..6e03913 100644 --- a/sora/body/core.py +++ b/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: diff --git a/sora/body/frame/iau_report.py b/sora/body/frame/iau_report.py index 696f286..e291b29 100644 --- a/sora/body/frame/iau_report.py +++ b/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,8 +72,8 @@ '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 @@ -81,19 +81,30 @@ '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": { diff --git a/sora/body/meta.py b/sora/body/meta.py index 56cf98e..c518730 100644 --- a/sora/body/meta.py +++ b/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: diff --git a/sora/body/shape/limb.py b/sora/body/shape/limb.py index f191650..72c86ff 100644 --- a/sora/body/shape/limb.py +++ b/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): """ diff --git a/sora/body/shape/utils.py b/sora/body/shape/utils.py index 972d28d..7541df3 100644 --- a/sora/body/shape/utils.py +++ b/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') diff --git a/sora/body/utils.py b/sora/body/utils.py index 3872f1c..984283e 100644 --- a/sora/body/utils.py +++ b/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") diff --git a/sora/occultation/core.py b/sora/occultation/core.py index ac5cdcb..1e3da7e 100644 --- a/sora/occultation/core.py +++ b/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 diff --git a/sora/occultation/utils.py b/sora/occultation/utils.py index 2712abb..8506039 100644 --- a/sora/occultation/utils.py +++ b/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 - ) \ No newline at end of file + ) + + +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 diff --git a/sora/prediction/occmap.py b/sora/prediction/occmap.py index 305eb36..5f4d181 100644 --- a/sora/prediction/occmap.py +++ b/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)