From 52cebe9b994ac10fa55044b2969b6d4b2689f9e2 Mon Sep 17 00:00:00 2001 From: "Acciarini, Giacomo (PG/R - Maths & Physics)" Date: Wed, 16 Apr 2025 13:51:29 +0200 Subject: [PATCH 1/4] refactored newton method, removal of unused functions and some other updates --- doc/api.rst | 2 - .../gradient_based_optimization.ipynb | 23 +- dsgp4/newton_method.py | 322 ++++++------------ dsgp4/util.py | 177 ++++------ tests/test_differentiability.py | 19 +- tests/test_newton.py | 6 +- tests/test_utils.py | 57 +++- 7 files changed, 266 insertions(+), 340 deletions(-) diff --git a/doc/api.rst b/doc/api.rst index a10fd81..045a751 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -42,9 +42,7 @@ $\partial$SGP4 API util.from_datetime_to_fractional_day util.from_datetime_to_mjd util.from_datetime_to_jd - util.from_cartesian_to_tle_elements util.from_cartesian_to_keplerian - util.from_cartesian_to_keplerian_torch sgp4 sgp4_batched sgp4init.sgp4init diff --git a/doc/notebooks/gradient_based_optimization.ipynb b/doc/notebooks/gradient_based_optimization.ipynb index 207c070..821069f 100644 --- a/doc/notebooks/gradient_based_optimization.ipynb +++ b/doc/notebooks/gradient_based_optimization.ipynb @@ -48,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 1, "id": "95627789", "metadata": {}, "outputs": [], @@ -59,7 +59,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 2, "id": "91d2d6c5", "metadata": {}, "outputs": [], @@ -72,7 +72,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 8, "id": "c376c651", "metadata": {}, "outputs": [ @@ -80,8 +80,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "F(y): 1.0751821836135118e-13\n", - "Solution found, at iter: 6\n", + "Solution F(y) = 2.4867357777951112e-14, converged in 5 iterations\n", "TLE(\n", "0 COSMOS 2251 DEB\n", "1 34454U 93036SX 22068.91971155 .00000319 00000-0 11812-3 0 9996\n", @@ -95,14 +94,14 @@ ], "source": [ "#when I propagate to zero, I expect the returned TLE to be identical to my_tle:\n", - "found_tle, y=dsgp4.newton_method(tle_0=my_tle,time_mjd=my_tle.date_mjd,new_tol=1e-12,max_iter=100)\n", + "found_tle, y=dsgp4.newton_method(tle0=my_tle,time_mjd=my_tle.date_mjd,new_tol=1e-12,max_iter=10,verbose=True)\n", "\n", "print(my_tle,found_tle)" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 7, "id": "e0f85268", "metadata": {}, "outputs": [ @@ -110,8 +109,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "F(y): 4.547475677268469e-13\n", - "Solution found, at iter: 7\n", + "Solution F(y) = 4.547517012268725e-13, converged in 5 iterations\n", "TLE(\n", "0 COSMOS 2251 DEB\n", "1 34454U 93036SX 22068.91971155 .00000319 00000-0 11812-3 0 9996\n", @@ -125,10 +123,11 @@ ], "source": [ "#I now propagate until 1000 minutes after\n", - "found_tle, y=dsgp4.newton_method(tle_0=my_tle,\n", + "found_tle, y=dsgp4.newton_method(tle0=my_tle,\n", " time_mjd=my_tle.date_mjd+800.,\n", " new_tol=1e-12,\n", - " max_iter=20)\n", + " max_iter=10,\n", + " verbose=True)\n", "#Newton still converges, and the TLE is of course now different:\n", "print(my_tle, found_tle)" ] @@ -136,7 +135,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "dsgp4", "language": "python", "name": "python3" }, diff --git a/dsgp4/newton_method.py b/dsgp4/newton_method.py index c301438..6b30ad1 100644 --- a/dsgp4/newton_method.py +++ b/dsgp4/newton_method.py @@ -5,232 +5,130 @@ from . import util from .tle import TLE +def update_TLE(old_tle, y0): + xpdotp = 1440.0 / (2.0 * np.pi) + mean_motion = float(y0[4]) * xpdotp * (np.pi / 43200.0) + + tle_elements = { + 'b_star': old_tle._bstar, + 'raan': float(y0[5]), + 'eccentricity': float(y0[0]), + 'argument_of_perigee': float(y0[1]), + 'inclination': float(y0[2]), + 'mean_anomaly': float(y0[3]), + 'mean_motion': mean_motion, + 'mean_motion_first_derivative': old_tle.mean_motion_first_derivative, + 'mean_motion_second_derivative': old_tle.mean_motion_second_derivative, + 'epoch_days': old_tle.epoch_days, + 'epoch_year': old_tle.epoch_year, + 'classification': old_tle.classification, + 'satellite_catalog_number': old_tle.satellite_catalog_number, + 'ephemeris_type': old_tle.ephemeris_type, + 'international_designator': old_tle.international_designator, + 'revolution_number_at_epoch': old_tle.revolution_number_at_epoch, + 'element_number': old_tle.element_number, + } + + return TLE(tle_elements) + +def initial_guess_tle(time_mjd, tle_object, gravity_constant_name="wgs-84"): + #first, let's decompose the time into -> epoch of the year and days + datetime_obj=util.from_mjd_to_datetime(time_mjd) + epoch_days=util.from_datetime_to_fractional_day(datetime_obj) + #then we need to propagate the state, and extract the keplerian elements: + util.initialize_tle(tle_object) + tsince=(time_mjd-util.from_datetime_to_mjd(tle_object._epoch))*1440. + target_state=util.propagate(tle_object, tsince).detach().numpy()*1e3 + _,mu_earth,_,_,_,_,_,_=util.get_gravity_constants(gravity_constant_name) + mu_earth=float(mu_earth)*1e9 + kepl_el=util.from_cartesian_to_keplerian(target_state[0],target_state[1],mu_earth) + #we need to convert the keplerian elements to TLE elements: + data = dict( + satellite_catalog_number=tle_object.satellite_catalog_number, + classification=tle_object.classification, + international_designator=tle_object.international_designator, + epoch_year=datetime_obj.year, + epoch_days=epoch_days, + ephemeris_type=tle_object.ephemeris_type, + element_number=tle_object.element_number, + revolution_number_at_epoch=tle_object.revolution_number_at_epoch, + mean_motion=np.sqrt(mu_earth/((kepl_el[0])**(3.0))), + mean_motion_first_derivative=tle_object.mean_motion_first_derivative, + mean_motion_second_derivative=tle_object.mean_motion_second_derivative, + eccentricity=kepl_el[1], + inclination=kepl_el[2], + argument_of_perigee=kepl_el[4], + raan=kepl_el[3], + mean_anomaly=kepl_el[5], + b_star=tle_object.b_star) + return TLE(data) + def _propagate(x, tle_sat, tsince, gravity_constant_name="wgs-84"): - from .sgp4init import sgp4init - from .sgp4 import sgp4 whichconst=util.get_gravity_constants(gravity_constant_name) sgp4init(whichconst=whichconst, opsmode='i', satn=tle_sat.satellite_catalog_number, epoch=(tle_sat._jdsatepoch+tle_sat._jdsatepochF)-2433281.5, - xbstar=x[0], - xndot=x[1], - xnddot=x[2], - xecco=x[3], - xargpo=x[4], - xinclo=x[5], - xmo=x[6], - xno_kozai=x[7], - xnodeo=x[8], + xbstar=tle_sat._bstar, + xndot=tle_sat._ndot, + xnddot=tle_sat._nddot, + xecco=x[0], + xargpo=x[1], + xinclo=x[2], + xmo=x[3], + xno_kozai=x[4], + xnodeo=x[5], satellite=tle_sat) state=sgp4(tle_sat, tsince*torch.ones(1,1)) return state -def initial_guess(tle_0, time_mjd, target_state=None): - """ - This method takes an initial TLE and the time at which we want to propagate it, and returns - a set of parameter related to an initial guess useful to find the TLE observation that corresponds to the propagated state - - Parameters: - ---------------- - tle_0 (``dsgp4.tle.TLE``): starting TLE at time - new_date (``datetime.datetime``): new date of the TLE, as a datetime objec - target_state (``torch.tensor``): a 2x3 tensor for the position and velocity of the state. If None, then this is computed - by propagating `tle_0` at the `new_date`. - - Returns: - ---------------- - y0 (``torch.tensor``): initial guess for the TLE elements. In particular, y0 contains - the following elements (see SGP4 for a thorough description of these parameters): - * y0[0]: bstar (B* parameter) - * y0[1]: ndot (mean motion first derivative) - * y0[2]: nddot (mean motion second derivative) - * y0[3]: ecco (eccentricity) - * y0[4]: argpo (argument of perigee) - * y0[5]: inclo (inclination) - * y0[6]: mo (mean anomaly) - * y0[7]: no_kozai (mean anomaly in certain units) - * y0[8]: nodeo (right ascension of the ascending node) - - target_state (``torch.tensor``): this expresses the cartesian state as position & velocity in km & km/s, at the propagated time. - The objective of the Newton method is to find the TLE observation that corresponds to that, at the propagated - time. - - new_tle (``dsgp4.tle.TLE``): TLE constructed with `y0`, in order to find the TLE that corresponds to `target_state` - - tle_elements_0 (``dict``): dictionary used to construct `new_tle` - - """ - new_date=util.from_mjd_to_datetime(time_mjd) +def newton_method(tle0, time_mjd, max_iter=50, new_tol=1e-12, verbose=False, target_state=None): if target_state is None: - whichconst=util.get_gravity_constants("wgs-84") - sgp4init(whichconst=whichconst, - opsmode=tle_0._opsmode, - satn=tle_0.satellite_catalog_number, - epoch=(tle_0._jdsatepoch+tle_0._jdsatepochF)-2433281.5, - xbstar=tle_0._bstar, - xndot=tle_0._ndot, - xnddot=tle_0._nddot, - xecco=tle_0._ecco, - xargpo=tle_0._argpo, - xinclo=tle_0._inclo, - xmo=tle_0._mo, - xno_kozai=tle_0._no_kozai, - xnodeo=tle_0._nodeo, - satellite=tle_0) - tsince=(time_mjd-util.from_datetime_to_mjd(tle_0._epoch))*1440. - x=tsince*torch.ones(1,1,requires_grad=True)#torch.rand(1,1, requires_grad=True) - target_state=sgp4(tle_0, x) - #print(target_state.size(), target_state) - tle_elements_0=util.from_cartesian_to_tle_elements(target_state.detach().numpy()*1e3) - #you need to transform it to the appropriate units - xpdotp=1440.0 / (2.0 *np.pi) - tle_elements_0['epoch_year']=new_date.year - tle_elements_0['epoch_days']=util.from_datetime_to_fractional_day(new_date) - #tle_elements_0['no_kozai']=tle_elements_0['mean_motion']/(np.pi/43200.0)/xpdotp - #tle_elements_0['inclo']=tle_elements_0['inclination'] - #tle_elements_0['nodeo']=tle_elements_0['raan'] - #tle_elements_0['argpo']=tle_elements_0['argument_of_perigee'] - #tle_elements_0['mo']=tle_elements_0['mean_anomaly'] - #tle_elements_0['ecco']=tle_elements_0['eccentricity'] - tle_elements_0['mean_motion_first_derivative']=tle_0.mean_motion_first_derivative - tle_elements_0['mean_motion_second_derivative']=tle_0.mean_motion_second_derivative - tle_elements_0['classification']=tle_0.classification - tle_elements_0['international_designator']=tle_0.international_designator - tle_elements_0['satellite_catalog_number']=tle_0.satellite_catalog_number - tle_elements_0['ephemeris_type']=tle_0.ephemeris_type - tle_elements_0['element_number']=tle_0.element_number - tle_elements_0['b_star']=tle_0.b_star - tle_elements_0['revolution_number_at_epoch']=tle_0.revolution_number_at_epoch - new_tle=TLE(tle_elements_0) - y0=torch.tensor([new_tle._bstar, - new_tle._ndot, - new_tle._nddot, - new_tle._ecco, - new_tle._argpo, - new_tle._inclo, - new_tle._mo, - new_tle._no_kozai, - new_tle._nodeo, - 0.],requires_grad=True) - return y0, target_state, new_tle, tle_elements_0 - -def update_TLE(old_tle,y0): - """ - This method takes a TLE and an initial guess, and returns a TLE updated accordingly. It - is useful while doing Newton iterations. - - Parameters: - ---------------- - old_tle (``dsgp4.tle.TLE``): TLE corresponding to previous guess - y0 (``torch.tensor``): initial guess (see the docstrings for `initial_guess` to know the content of `y0`) - - Returns: - ---------------- - new_tle (``dsgp4.tle.TLE``): updated TLE - - """ - tle_elements={} - xpdotp=1440.0/(2.0 *np.pi) - #I need to convert no_kozai to consistent units for what is expected by TLE class: - no_kozai=float(y0[7])*xpdotp - tle_elements['mean_motion']=no_kozai*np.pi/43200.0 - - #the rest is consistent: - tle_elements['b_star']=float(y0[0]) - tle_elements['raan']=float(y0[8]) - tle_elements['eccentricity']=float(y0[3]) - tle_elements['argument_of_perigee']=float(y0[4]) - tle_elements['inclination']=float(y0[5]) - tle_elements['mean_anomaly']=float(y0[6]) - - #now the ones that stay the same: - tle_elements['mean_motion_first_derivative']=old_tle.mean_motion_first_derivative - tle_elements['mean_motion_second_derivative']=old_tle.mean_motion_second_derivative - tle_elements['epoch_days']=old_tle.epoch_days - tle_elements['epoch_year']=old_tle.epoch_year - tle_elements['classification']=old_tle.classification - tle_elements['satellite_catalog_number']=old_tle.satellite_catalog_number - tle_elements['ephemeris_type']=old_tle.ephemeris_type - tle_elements['international_designator']=old_tle.international_designator - tle_elements['revolution_number_at_epoch']=old_tle.revolution_number_at_epoch - tle_elements['element_number']=old_tle.element_number - return TLE(tle_elements) - -def newton_method(tle_0, time_mjd, target_state=None, new_tol=1e-12,max_iter=50, verbose=False): - """ - This method performs Newton method starting from an initial TLE and a given propagation time. The objective - is to find a TLE that accurately reconstructs the propagated state, at observation time. - - Parameters: - ---------------- - tle_0 (``dsgp4.tle.TLE``): starting TLE (i.e., TLE at a given initial time) - time_mjd (``float``): time at which we want to propagate the TLE (in Modified Julian Date) - target_state (``torch.tensor``): a 2x3 tensor for the position and velocity of the state. If None, then this is computed by propagating `tle_0` at the `new_date`. - new_tol (``float``): tolerance for the Newton method - max_iter (``int``): maximum number of iterations for the Newton method - - Returns: - ---------------- - next_tle (``dsgp4.tle.TLE``): TLE corresponding to the propagated state - y0 (``torch.tensor``): initial guess for the TLE elements. In particular, y0 contains the following elements (see SGP4 for a thorough description of these parameters): + util.initialize_tle(tle0) + target_state=util.propagate(tle0, (time_mjd-util.from_datetime_to_mjd(tle0._epoch))*1440.) - """ - i=0 - tol=1e9 - y0, state_target, next_tle, tle_elements_0=initial_guess(tle_0=tle_0,time_mjd=time_mjd,target_state=target_state) - #print(y0,states) - #newton iterations: + i,tol=0,1e9 + next_tle=initial_guess_tle(time_mjd, tle0) + y0 = torch.tensor([ + next_tle._ecco, + next_tle._argpo, + next_tle._inclo, + next_tle._mo, + next_tle._no_kozai, + next_tle._nodeo, + ], requires_grad=True) + def propagate_fn(x): + tsince=(time_mjd-util.from_datetime_to_mjd(next_tle._epoch))*1440. + return _propagate(x,next_tle,tsince) while inew_tol: - #print(f"y0: {y0}") - propagate=lambda x: _propagate(x,next_tle,(time_mjd-util.from_datetime_to_mjd(next_tle._epoch))*1440.) - y1=util.clone_w_grad(y0) - y2=util.clone_w_grad(y0) - y3=util.clone_w_grad(y0) - y4=util.clone_w_grad(y0) - y5=util.clone_w_grad(y0) - r_x=propagate(y0)[0][0] - r_x.backward() - gradient_rx=y0.grad - r_y=propagate(y1)[0][1] - r_y.backward() - gradient_ry=y1.grad - r_z=propagate(y2)[0][2] - r_z.backward() - gradient_rz=y2.grad - v_x=propagate(y3)[1][0] - v_x.backward() - gradient_vx=y3.grad - v_y=propagate(y4)[1][1] - v_y.backward() - gradient_vy=y4.grad - v_z=propagate(y5)[1][2] - v_z.backward() - gradient_vz=y5.grad - F=np.array([float(r_x)-float(state_target[0][0]),float(r_y)-float(state_target[0][1]), float(r_z)-float(state_target[0][2]), float(v_x)-float(state_target[1][0]), float(v_y)-float(state_target[1][1]), float(v_z)-float(state_target[1][2])]) + grads=[] + F=[] + for idx, (row,col) in enumerate([(0,0), (0,1), (0,2), (1,0), (1,1), (1,2)]): + y=util.clone_w_grad(y0) + val=propagate_fn(y)[row][col] + val.backward() + grads.append(y.grad) + F.append((val-target_state[row][col]).item()) tol=np.linalg.norm(F) - #print(f"|F|: {tol}") - DF=np.stack((gradient_rx, gradient_ry, gradient_rz, gradient_vx, gradient_vy, gradient_vz)) - #we remove the first three columns (they are all zeros): - DF=DF[:,3:] - DF=DF[:,:-1] - dY=-np.matmul(np.matmul(np.linalg.pinv(np.matmul(DF.T,DF)),DF.T),F) - #we make sure the eccentricity does not go to negative values: - if float(y0[3])+dY[0]<0: - dY[0]=-float(y0[3])*0.9999 - dY=torch.tensor([0.,0.,0.]+list(dY)+[0.], requires_grad=True) if tol1.: + dY[0]=1-1e-10 + dY=torch.tensor(dY, requires_grad=True) + #update the state: + y0 = torch.tensor([float(a) + float(b) for a, b in zip(y0, dY)], requires_grad=True) + next_tle = update_TLE(next_tle, y0) i+=1 - if verbose: + if verbose: print("Solution not found, returning best found so far") - print(f"F(y): {np.linalg.norm(F)}") - return next_tle, y0 + print(f"F(y): {tol:.3e}") + return next_tle, y0 \ No newline at end of file diff --git a/dsgp4/util.py b/dsgp4/util.py index 74e1442..f9f6888 100644 --- a/dsgp4/util.py +++ b/dsgp4/util.py @@ -416,120 +416,79 @@ def from_datetime_to_jd(datetime_obj): """ return sum(jday(year=datetime_obj.year, mon=datetime_obj.month, day=datetime_obj.day, hr=datetime_obj.hour, minute=datetime_obj.minute, sec=datetime_obj.second+float('0.'+str(datetime_obj.microsecond)))) -def from_cartesian_to_tle_elements(state, gravity_constant_name='wgs-72'): +def from_cartesian_to_keplerian(r_vec, v_vec, mu): """ - This function converts the provided state from Cartesian to TLE elements. - - Parameters: - ---------------- - state (``np.array``): state in Cartesian coordinates - gravity_constant_name (``str``): name of the gravity constant to be used (default: 'wgs-72') - - Returns: - ---------------- - ``dict``: dictionary of TLE elements - """ - _,mu_earth,_,_,_,_,_,_=get_gravity_constants(gravity_constant_name) - mu_earth=float(mu_earth)*1e9 - kepl_el = from_cartesian_to_keplerian(state, mu_earth) - tle_elements={} - tle_elements['mean_motion'] = np.sqrt(mu_earth/((kepl_el[0])**(3.0))) - tle_elements['eccentricity'] = kepl_el[1] - tle_elements['inclination'] = kepl_el[2] - tle_elements['raan'] = kepl_el[3] - tle_elements['argument_of_perigee'] = kepl_el[4] - mean_anomaly = kepl_el[5] - kepl_el[1]*np.sin(kepl_el[5]) - tle_elements['mean_anomaly'] = mean_anomaly%(2*np.pi) - return tle_elements - -def from_cartesian_to_keplerian(state, mu): - """ - This function takes the state in cartesian coordinates and the gravitational - parameter of the central body, and returns the state in Keplerian elements. - - Parameters: - ---------------- - state (``np.array``): numpy array of 2 rows and 3 columns, where - the first row represents position, and the second velocity. - mu (``float``): gravitational parameter of the central body + This function converts the provided state from Cartesian to Keplerian elements. - Returns: - ---------------- - ``np.array``: numpy array of the six keplerian elements: (a,e,i,omega,Omega,mean_anomaly) - (i.e., semi major axis, eccentricity, inclination, - right ascension of ascending node, argument of perigee, - mean anomaly). All the angles are in radiants, eccentricity is unitless - and semi major axis is in SI. - """ - h_bar = np.cross(np.array([state[0,0], state[0,1], state[0,2]]), np.array([state[1,0], state[1,1], state[1,2]])) - h = np.linalg.norm(h_bar) - r = np.linalg.norm(np.array([state[0,0], state[0,1], state[0,2]])) - v = np.linalg.norm(np.array([state[1,0], state[1,1], state[1,2]])) - E = 0.5*(v**2)-mu/r - a = -mu/(2*E) - e = np.sqrt(1-(h**2)/(a*mu)) - i = np.arccos(h_bar[2]/h) - Omega = np.arctan2(h_bar[0],-h_bar[1]) - - lat = np.arctan2(np.divide(state[0,2],(np.sin(i))), (state[0,0]*np.cos(Omega) + state[0,1]*np.sin(Omega))) - p = a*(1-e**2) - nu = np.arctan2(np.sqrt(p/mu)*np.dot(np.array([state[0,0], state[0,1], state[0,2]]),np.array([state[1,0], state[1,1], state[1,2]])), p-r) - omega = (lat-nu) - eccentric_anomaly = 2*np.arctan(np.sqrt((1-e)/(1+e))*np.tan(nu/2)) - n = np.sqrt(mu/(a**3)) - mean_anomaly=eccentric_anomaly-e*np.sin(eccentric_anomaly) - #I make sure they are always in 0,2pi - if mean_anomaly<0: - mean_anomaly = 2*np.pi-abs(mean_anomaly) - if omega<0: - omega=2*np.pi-abs(omega) - if Omega<0: - Omega=2*np.pi-abs(Omega) - if abs(mean_anomaly)>2*np.pi: - mean_anomaly=mean_anomaly%(2*np.pi) - if abs(omega)>2*np.pi: - omega=omega%(2*np.pi) - if abs(Omega)>2*np.pi: - Omega=Omega%(2*np.pi) - return np.array([a, e, i, Omega, omega, mean_anomaly]) - -def from_cartesian_to_keplerian_torch(state, mu): - """ - Same as from_cartesian_to_keplerian, but for torch tensors. - Parameters: ---------------- - state (``torch.tensor``): torch tensor of 2 rows and 3 columns, where - the first row represents position, and the second velocity. - mu (``float``): gravitational parameter of the central body + r_vec (``np.array``): position vector in Cartesian coordinates + v_vec (``np.array``): velocity vector in Cartesian coordinates + mu (``float``): gravitational parameter of the central body Returns: ---------------- - ``torch.tensor``: torch tensor of the six keplerian elements: (a,e,i,omega,Omega,mean_anomaly) - (i.e., semi major axis, eccentricity, inclination, + ``np.array``: array of Keplerian elements: (a, e, i, Omega, omega, M) + (i.e., semi-major axis, eccentricity, inclination, right ascension of ascending node, argument of perigee, - mean anomaly). All the angles are in radiants, eccentricity is unitless - and semi major axis is in SI. - """ - h_bar = torch.cross(state[0], state[1]) - h = h_bar.norm() - r = state[0].norm() - v = torch.norm(state[1]) - E = 0.5*(v**2)-mu/r - a = -mu/(2*E) - e = torch.sqrt(1-(h**2)/(a*mu)) - i = torch.arccos(h_bar[2]/h) - Omega = torch.arctan2(h_bar[0],-h_bar[1]) - lat = torch.arctan2(torch.divide(state[0,2],(torch.sin(i))), (state[0,0]*torch.cos(Omega) + state[0,1]*torch.sin(Omega))) - p = a*(1-e**2) - nu = torch.arctan2(torch.sqrt(p/mu)*torch.dot(state[0],state[1]), p-r) - omega = (lat-nu) - eccentric_anomaly = 2*torch.arctan(torch.sqrt((1-e)/(1+e))*torch.tan(nu/2)) - n = torch.sqrt(mu/(a**3)) - mean_anomaly=eccentric_anomaly-e*torch.sin(eccentric_anomaly) - #I make sure they are always in 0,2pi - mean_motion=torch.sqrt(mu/((a)**(3.0))) - xpdotp = 1440.0 / (2.0 *np.pi) - no_kozai_conversion_factor=xpdotp/43200.0* np.pi - no_kozai=mean_motion/no_kozai_conversion_factor - return [no_kozai, e, i, Omega, omega, mean_anomaly] + mean anomaly). All the angles are in radians, eccentricity is unitless + and semi-major axis is in SI. + """ + # Norms + r = np.linalg.norm(r_vec) + v = np.linalg.norm(v_vec) + + # Angular momentum vector and its magnitude + h_vec = np.cross(r_vec, v_vec) + h = np.linalg.norm(h_vec) + + # Inclination + i = np.arccos(h_vec[2] / h) if h != 0 else 0 + + # Node vector + K = np.array([0, 0, 1]) + n_vec = np.cross(K, h_vec) + n = np.linalg.norm(n_vec) + + # Eccentricity vector + e_vec = (1/mu) * ((v**2 - mu/r) * r_vec - np.dot(r_vec, v_vec) * v_vec) + e = np.linalg.norm(e_vec) + + # Semi-major axis (a) + energy = v**2 / 2 - mu / r + if abs(e - 1.0) > 1e-8: # Elliptical or hyperbolic orbit + a = -mu / (2 * energy) + else: # Parabolic orbit + a = np.inf + + # Right ascension of ascending node (RAAN) + Omega = 0 + if n != 0: + Omega = np.arccos(n_vec[0] / n) if n_vec[1] >= 0 else 2 * np.pi - np.arccos(n_vec[0] / n) + + # Argument of perigee (ω) + omega = 0 + if n != 0 and e != 0: + omega = np.arccos(np.dot(n_vec, e_vec) / (n * e)) if e_vec[2] >= 0 else 2 * np.pi - np.arccos(np.dot(n_vec, e_vec) / (n * e)) + + # True anomaly (ν) + nu = 0 + if e != 0: + nu = np.arccos(np.dot(e_vec, r_vec) / (e * r)) if np.dot(r_vec, v_vec) >= 0 else 2 * np.pi - np.arccos(np.dot(e_vec, r_vec) / (e * r)) + + # Eccentric anomaly (E) and Mean anomaly (M) + if e < 1.0: # Elliptical orbit + E = 2 * np.arctan(np.tan(nu / 2) * np.sqrt((1 - e) / (1 + e))) + M = E - e * np.sin(E) + elif e > 1.0: # Hyperbolic orbit + F = 2 * np.arctanh(np.tan(nu / 2) * np.sqrt((e - 1) / (e + 1))) + M = e * np.sinh(F) - F + else: # Parabolic orbit + M = np.nan # Mean anomaly is undefined for parabolic orbits + + # Normalize angles to [0, 2π) + Omega %= 2 * np.pi + omega %= 2 * np.pi + M %= 2 * np.pi + + return np.array([a, e, i, Omega, omega, M]) diff --git a/tests/test_differentiability.py b/tests/test_differentiability.py index 376cbfc..1869c33 100644 --- a/tests/test_differentiability.py +++ b/tests/test_differentiability.py @@ -5,7 +5,6 @@ import sgp4.io import torch import unittest -from dsgp4.newton_method import _propagate error_string="Error: deep space propagation not supported (yet). The provided satellite has \ an orbital period above 225 minutes. If you want to let us know you need it or you want to \ @@ -87,6 +86,24 @@ def test_input_gradients(self): state_teme.flatten()[i].backward(retain_graph=True) partial_derivatives[i,:] = tles_elements[idx].grad #now, let's compare with finite differences: + def _propagate(x, tle_sat, tsince, gravity_constant_name="wgs-84"): + whichconst=dsgp4.util.get_gravity_constants(gravity_constant_name) + dsgp4.sgp4init(whichconst=whichconst, + opsmode='i', + satn=tle_sat.satellite_catalog_number, + epoch=(tle_sat._jdsatepoch+tle_sat._jdsatepochF)-2433281.5, + xbstar=x[0], + xndot=x[1], + xnddot=x[2], + xecco=x[3], + xargpo=x[4], + xinclo=x[5], + xmo=x[6], + xno_kozai=x[7], + xnodeo=x[8], + satellite=tle_sat) + state=dsgp4.sgp4(tle_sat, tsince*torch.ones(1,1)) + return state fun=lambda xx: _propagate(xx,tle,time) #I create 6 copies of the inputs to be used for each function (three #position and three velocity coordinates) diff --git a/tests/test_newton.py b/tests/test_newton.py index 6af9e1d..49c65b1 100644 --- a/tests/test_newton.py +++ b/tests/test_newton.py @@ -11,7 +11,7 @@ def test_newton_method_1(self): lines.append('1 43437U 18100A 20143.90384230 .00041418 00000-0 10000-3 0 99968') lines.append('2 43437 97.8268 249.9127 0221000 123.9136 259.1144 15.12608579563539') my_tle=dsgp4.tle.TLE(lines) - whichconst=dsgp4.util.get_gravity_constants("wgs-72") + whichconst=dsgp4.util.get_gravity_constants("wgs-84") dsgp4.sgp4init(whichconst=whichconst, opsmode=my_tle._opsmode, satn=my_tle.satellite_catalog_number, @@ -27,7 +27,7 @@ def test_newton_method_1(self): xnodeo=my_tle._nodeo, satellite=my_tle) - found_tle,_=dsgp4.newton_method(tle_0=my_tle, time_mjd=dsgp4.util.from_datetime_to_mjd(my_tle._epoch)) + found_tle,_=dsgp4.newton_method(tle0=my_tle, time_mjd=dsgp4.util.from_datetime_to_mjd(my_tle._epoch)) self.assertEqual(my_tle.line1, found_tle.line1) self.assertEqual(my_tle.line2, found_tle.line2) #I now propagate for 1000 minutes and make sure that the found TLE w Newton @@ -37,6 +37,6 @@ def test_newton_method_1(self): target_state=dsgp4.sgp4(my_tle,tsince_2*torch.ones(1,1,requires_grad=True)) time_mjd_0=dsgp4.util.from_datetime_to_mjd(my_tle._epoch) time_mjd=time_mjd_0+tsince_2/24/60 - found_tle_2,_=dsgp4.newton_method(tle_0=my_tle, time_mjd=time_mjd, new_tol=new_tol, target_state=target_state) + found_tle_2,_=dsgp4.newton_method(tle0=my_tle, time_mjd=time_mjd, new_tol=new_tol, target_state=target_state) found_state=dsgp4.sgp4(found_tle_2,torch.tensor((time_mjd-dsgp4.util.from_datetime_to_mjd(found_tle_2._epoch))*1440.))#(dsgp4.util.from_datetime_to_mjd(found_tle_2._epoch)-dsgp4.util.from_datetime_to_mjd(new_date))*1440.*torch.ones(1,1,requires_grad=True)) self.assertTrue(torch.norm(found_state-target_state) Keplerian using our function + a, e, i, Omega, omega, M = dsgp4.util.from_cartesian_to_keplerian(r_vec, v_vec, mu) + + # # Use poliastro to compute the same transformations + # # Create an orbit from Cartesian state vectors + # orbit = Orbit.from_vectors(Earth, + # r_vec * u.m, + # v_vec * u.m / u.s) + # # Extract Keplerian elements from poliastro + # a_poliastro = orbit.a.to(u.m).value + # e_poliastro = orbit.ecc.value + # i_poliastro = orbit.inc.to(u.rad).value + # Omega_poliastro = orbit.raan.to(u.rad).value + # omega_poliastro = orbit.argp.to(u.rad).value + # nu_poliastro = orbit.nu.to(u.rad).value # True anomaly + # # Compute mean anomaly (M) from true anomaly (ν) + # if e < 1.0: # Elliptical orbit + # E = 2 * np.arctan(np.sqrt((1 - e) / (1 + e)) * np.tan(nu_poliastro / 2)) + # M_poliastro = E - e * np.sin(E) + # else: # Hyperbolic orbit + # F = 2 * np.arctanh(np.sqrt((e - 1) / (e + 1)) * np.tanh(nu_poliastro / 2)) + # M_poliastro = e * np.sinh(F) - F + # Omega_poliastro %= 2 * np.pi + # omega_poliastro %= 2 * np.pi + # M_poliastro %= 2 * np.pi + a_poliastro=7023679.5817881366237998 + e_poliastro=0.0041649630912143 + i_poliastro=1.2919744331609118 + Omega_poliastro=5.3551396410293757 + omega_poliastro=0.4409272281996022 + M_poliastro=5.8458093777349722 + + # Compare the results + self.assertAlmostEqual(a, a_poliastro, places=6) + self.assertAlmostEqual(e, e_poliastro, places=6) + self.assertAlmostEqual(i, i_poliastro, places=6) + self.assertAlmostEqual(Omega, Omega_poliastro, places=6) + self.assertAlmostEqual(omega, omega_poliastro, places=6) + self.assertAlmostEqual(M, M_poliastro, places=6) From 393569537206d5f7a6f26d27747fe2f65baa0b56 Mon Sep 17 00:00:00 2001 From: "Acciarini, Giacomo (PG/R - Maths & Physics)" Date: Wed, 16 Apr 2025 14:01:04 +0200 Subject: [PATCH 2/4] update docstrings --- dsgp4/newton_method.py | 56 ++++++++++++++++++++++++++++++++++++++---- tests/test_newton.py | 7 +++--- 2 files changed, 55 insertions(+), 8 deletions(-) diff --git a/dsgp4/newton_method.py b/dsgp4/newton_method.py index 6b30ad1..5690e18 100644 --- a/dsgp4/newton_method.py +++ b/dsgp4/newton_method.py @@ -6,6 +6,18 @@ from .tle import TLE def update_TLE(old_tle, y0): + """ + This function updates the TLE object with the new keplerian elements. + + Parameters: + ---------------- + old_tle (``dsgp4.TLE``): The old TLE object to be updated. + y0 (``torch.tensor``): The new keplerian elements. + + Returns: + ---------------- + TLE: The updated TLE object. + """ xpdotp = 1440.0 / (2.0 * np.pi) mean_motion = float(y0[4]) * xpdotp * (np.pi / 43200.0) @@ -32,11 +44,26 @@ def update_TLE(old_tle, y0): return TLE(tle_elements) def initial_guess_tle(time_mjd, tle_object, gravity_constant_name="wgs-84"): + """ + This function generates an initial guess for the TLE object at a given time. + It uses the current TLE object to propagate the state and extract the keplerian elements. + The function returns a new TLE object with the updated elements. + + Parameters: + ---------------- + time_mjd (``float``): The time in MJD format. + tle_object (``dsgp4.TLE``): The TLE object to be updated. + gravity_constant_name (``str``): The name of the gravity constant to be used. Default is "wgs-84". + + Returns: + ---------------- + TLE: The updated TLE object. + """ #first, let's decompose the time into -> epoch of the year and days datetime_obj=util.from_mjd_to_datetime(time_mjd) epoch_days=util.from_datetime_to_fractional_day(datetime_obj) #then we need to propagate the state, and extract the keplerian elements: - util.initialize_tle(tle_object) + util.initialize_tle(tle_object,gravity_constant_name=gravity_constant_name) tsince=(time_mjd-util.from_datetime_to_mjd(tle_object._epoch))*1440. target_state=util.propagate(tle_object, tsince).detach().numpy()*1e3 _,mu_earth,_,_,_,_,_,_=util.get_gravity_constants(gravity_constant_name) @@ -82,13 +109,32 @@ def _propagate(x, tle_sat, tsince, gravity_constant_name="wgs-84"): state=sgp4(tle_sat, tsince*torch.ones(1,1)) return state -def newton_method(tle0, time_mjd, max_iter=50, new_tol=1e-12, verbose=False, target_state=None): +def newton_method(tle0, time_mjd, max_iter=50, new_tol=1e-12, verbose=False, target_state=None, gravity_constant_name="wgs-84"): + """ + This function implements the Newton-Raphson method to find the TLE elements that match a given target state. + It uses the SGP4 propagator to propagate the TLE elements and compare them with the target state. + The function returns the updated TLE object and the final state vector. + Parameters: + ---------------- + tle0 (``dsgp4.TLE``): The initial TLE object to be updated. + time_mjd (``float``): The time in MJD format. + max_iter (``int``): The maximum number of iterations for the Newton-Raphson method. Default is 50. + new_tol (``float``): The tolerance for convergence. Default is 1e-12. + verbose (``bool``): If True, prints the convergence information. Default is False. + target_state (``torch.tensor``): The target state vector to be matched. If None, the function will propagate the initial TLE object to find the target state. + gravity_constant_name (``str``): The name of the gravity constant to be used. Default is "wgs-84". + + Returns: + ---------------- + TLE: The updated TLE object. + torch.tensor: The tensor of the updated elements. + """ if target_state is None: - util.initialize_tle(tle0) + util.initialize_tle(tle0,gravity_constant_name=gravity_constant_name) target_state=util.propagate(tle0, (time_mjd-util.from_datetime_to_mjd(tle0._epoch))*1440.) i,tol=0,1e9 - next_tle=initial_guess_tle(time_mjd, tle0) + next_tle=initial_guess_tle(time_mjd, tle0,gravity_constant_name=gravity_constant_name) y0 = torch.tensor([ next_tle._ecco, next_tle._argpo, @@ -99,7 +145,7 @@ def newton_method(tle0, time_mjd, max_iter=50, new_tol=1e-12, verbose=False, tar ], requires_grad=True) def propagate_fn(x): tsince=(time_mjd-util.from_datetime_to_mjd(next_tle._epoch))*1440. - return _propagate(x,next_tle,tsince) + return _propagate(x,next_tle,tsince,gravity_constant_name=gravity_constant_name) while inew_tol: grads=[] F=[] diff --git a/tests/test_newton.py b/tests/test_newton.py index 49c65b1..11773c8 100644 --- a/tests/test_newton.py +++ b/tests/test_newton.py @@ -11,7 +11,8 @@ def test_newton_method_1(self): lines.append('1 43437U 18100A 20143.90384230 .00041418 00000-0 10000-3 0 99968') lines.append('2 43437 97.8268 249.9127 0221000 123.9136 259.1144 15.12608579563539') my_tle=dsgp4.tle.TLE(lines) - whichconst=dsgp4.util.get_gravity_constants("wgs-84") + gravity_constant_name="wgs-84" + whichconst=dsgp4.util.get_gravity_constants(gravity_constant_name) dsgp4.sgp4init(whichconst=whichconst, opsmode=my_tle._opsmode, satn=my_tle.satellite_catalog_number, @@ -27,7 +28,7 @@ def test_newton_method_1(self): xnodeo=my_tle._nodeo, satellite=my_tle) - found_tle,_=dsgp4.newton_method(tle0=my_tle, time_mjd=dsgp4.util.from_datetime_to_mjd(my_tle._epoch)) + found_tle,_=dsgp4.newton_method(tle0=my_tle, time_mjd=dsgp4.util.from_datetime_to_mjd(my_tle._epoch),gravity_constant_name=gravity_constant_name) self.assertEqual(my_tle.line1, found_tle.line1) self.assertEqual(my_tle.line2, found_tle.line2) #I now propagate for 1000 minutes and make sure that the found TLE w Newton @@ -37,6 +38,6 @@ def test_newton_method_1(self): target_state=dsgp4.sgp4(my_tle,tsince_2*torch.ones(1,1,requires_grad=True)) time_mjd_0=dsgp4.util.from_datetime_to_mjd(my_tle._epoch) time_mjd=time_mjd_0+tsince_2/24/60 - found_tle_2,_=dsgp4.newton_method(tle0=my_tle, time_mjd=time_mjd, new_tol=new_tol, target_state=target_state) + found_tle_2,_=dsgp4.newton_method(tle0=my_tle, time_mjd=time_mjd, new_tol=new_tol, target_state=target_state,gravity_constant_name=gravity_constant_name) found_state=dsgp4.sgp4(found_tle_2,torch.tensor((time_mjd-dsgp4.util.from_datetime_to_mjd(found_tle_2._epoch))*1440.))#(dsgp4.util.from_datetime_to_mjd(found_tle_2._epoch)-dsgp4.util.from_datetime_to_mjd(new_date))*1440.*torch.ones(1,1,requires_grad=True)) self.assertTrue(torch.norm(found_state-target_state) Date: Wed, 16 Apr 2025 14:02:55 +0200 Subject: [PATCH 3/4] fix linkcheck for esa website --- doc/conf.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/doc/conf.py b/doc/conf.py index e29fa08..89d27d3 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -85,6 +85,13 @@ html_logo = "_static/logo_dsgp4.png" + +linkcheck_ignore = [ + r'https://www\.esa\.int/gsp/ACT/team/giacomo_acciarini/', + r'https://www.esa.int/gsp/ACT/team/dario_izzo/', + ] + + html_theme_options = { "repository_url": "https://github.com/esa/dSGP4/", "repository_branch": "master", From 34142aa6878d06b9aa7c04f8de892d6d09579ae0 Mon Sep 17 00:00:00 2001 From: "Acciarini, Giacomo (PG/R - Maths & Physics)" Date: Wed, 16 Apr 2025 14:04:48 +0200 Subject: [PATCH 4/4] update version --- dsgp4/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsgp4/__init__.py b/dsgp4/__init__.py index 815ad2e..1b8aa91 100644 --- a/dsgp4/__init__.py +++ b/dsgp4/__init__.py @@ -1,4 +1,4 @@ -__version__ = '1.1.4' +__version__ = '1.1.5' import torch torch.set_default_dtype(torch.float64)