LRO Solar Radiation Pressure

This example demonstrates how to compute solar radiation pressure (SRP) accelerations using the object-oriented interface of the pyRTX library.

This approach computes SRP accelerations on-the-fly using the spacecraft geometry and solar position at each time step.

lro_srp.py
  1### ------------------------------------------------------------------------------------------------------- ###
  2
  3# Example purpose:
  4# Show the object-oriented interface of the pyRTX library
  5#
  6# Example case:
  7# Compute the srp acceleration for LRO spacecraft, using the SPICE trajectory and frames
  8
  9### ------------------------------------------------------------------------------------------------------- ###
 10### IMPORTS
 11
 12import spiceypy as sp
 13import xarray as xr
 14import matplotlib.pyplot as plt
 15import logging, timeit
 16
 17from pyRTX.classes.Spacecraft import Spacecraft
 18from pyRTX.classes.Planet import Planet
 19from pyRTX.classes.PixelPlane import PixelPlane
 20from pyRTX.classes.RayTracer import RayTracer
 21from pyRTX.classes.SRP import SunShadow, SolarPressure 
 22from pyRTX.classes.Precompute import Precompute
 23from pyRTX.core.analysis_utils import epochRange2
 24import logging
 25
 26from numpy import floor, mod
 27
 28import warnings
 29warnings.filterwarnings('ignore')
 30
 31### ------------------------------------------------------------------------------------------------------- ###
 32### INPUTS
 33
 34ref_epc 	=  "2010 may 10 09:25:00"
 35duration    =  10000  									  # seconds
 36timestep    =  100
 37spacing	    =  0.01
 38METAKR      =  '../example_data/LRO/metakernel_lro.tm'     # metakernel
 39obj_path    =  '../example_data/LRO/'				       # folder with shape .obj files
 40base_flux   =  1361.5
 41ref_radius  =  1737.4
 42n_cores	    =  10
 43
 44# The spacecraft mass can be a float, int or a xarray with times and values [kg]
 45sc_mass = xr.open_dataset('mass/lro_mass.nc')
 46sc_mass.load()
 47sc_mass.close()
 48
 49### ------------------------------------------------------------------------------------------------------- ###
 50### OBJECTS DEFINITION
 51
 52# Time initialization
 53tic = timeit.default_timer()
 54
 55# Load the metakernel containing references to the necessary SPICE frames
 56sp.furnsh(METAKR)
 57
 58# Define a basic epoch
 59epc_et0 = sp.str2et( ref_epc ) 
 60epc_et1 = epc_et0 + duration
 61epochs  = epochRange2(startEpoch = epc_et0, endEpoch = epc_et1, step = timestep)
 62
 63# Define the Spacecraft Object (Refer to the class documentation for further details)
 64lro = Spacecraft( name = 'LRO',
 65                 
 66				  base_frame = 'LRO_SC_BUS', 					     # Name of the spacecraft body-fixed frame
 67      
 68                  mass = sc_mass,
 69      
 70				  spacecraft_model = {						         # Define a spacecraft model
 71                          
 72					'LRO_BUS': { 
 73							 'file' : obj_path + 'bus_rotated.obj',	 # .obj file of the spacecraft component
 74							 'frame_type': 'Spice',				     # type of frame (can be 'Spice' or 'UD'
 75							 'frame_name': 'LRO_SC_BUS',			 # Name of the frame
 76							 'center': [0.0,0.0,0.0],			     # Origin of the component
 77							 'diffuse': 0.1,				         # Diffuse reflect. coefficient
 78							 'specular': 0.3,				         # Specular reflect. coefficient
 79							 },
 80
 81					'LRO_SA': {	
 82							'file': obj_path + 'SA_recentred.obj',
 83							'frame_type': 'Spice',
 84							'frame_name': 'LRO_SA',
 85							'center': [-1,-1.1, -0.1],
 86							'diffuse': 0,
 87							'specular': 0.3,
 88							},
 89
 90
 91					'LRO_HGA': { 	
 92							'file': obj_path + 'HGA_recentred.obj',
 93							'frame_type': 'Spice',
 94							'frame_name': 'LRO_HGA',
 95							'center':[-0.99,    -0.3,  -3.1],
 96							'diffuse': 0.2,
 97							'specular': 0.1,
 98							},
 99					}
100				)
101
102
103# Define the Moon object
104moon = Planet(  fromFile      = None,
105                radius        = ref_radius,
106                name          = 'Moon',
107                bodyFrame     = 'MOON_PA',
108                sunFixedFrame = 'GSE_MOON',
109                units         = 'km',
110                subdivs       = 5,
111                )
112
113
114# Define the Sun rays object
115rays = PixelPlane( spacecraft  = lro,         # Spacecraft object 
116			       mode        = 'Dynamic',   # Mode: can be 'Dynamic' ( The sun orientation is computed from the kernels), or 'Fixed'
117			       distance    = 100,	      # Distance of the ray origin from the spacecraft
118			       source      = 'Sun',       # Source body (used to compute the orientation of the rays wrt. spacecraft)
119			       width       = 10,	      # Width of the pixel plane
120			       height      = 10,          # Height of the pixel plane
121			       ray_spacing = spacing,     # Ray spacing (in m)
122				   )
123
124
125# Define the ray tracer
126rtx = RayTracer( lro,                    # Spacecraft object
127                 rays,                   # pixelPlane object
128                 kernel = 'Embree3',     # The RTX kernel to use
129                 bounces = 2,            # The number of bounces to account for
130                 diffusion = False,      # Account for secondary diffusion
131                ) 
132
133# Precomputation object. This object performs all the calls to spiceypy before 
134# calculating the acceleration. This is necessary when calculating the acceleration
135# with parallel cores.
136prec = Precompute(epochs = epochs,)
137prec.precomputeSolarPressure(lro, moon, correction='LT+S')
138prec.dump()
139
140# Define the shadow function object
141shadow = SunShadow( spacecraft     = lro,
142				    body           = 'Moon',
143				    bodyShape      = moon,
144				    limbDarkening  = 'Eddington',
145        			precomputation = prec,
146				    )
147
148# Define the solar pressure object
149srp = SolarPressure( lro, 
150				     rtx,
151				     baseflux       = base_flux,    # Here we use the None option to obtain the generalized geometry vector, used also for the computation of albedo and thermal infrared
152				     shadowObj      = shadow,
153					 precomputation = prec,
154				     )
155
156# Managing Error messages from trimesh
157# (when concatenating textures, in this case, withouth .mtl definition, trimesh returns a warning that
158#  would fill the stdout. Deactivate it for a clean output)
159log = logging.getLogger('trimesh')
160log.disabled = True
161
162### ------------------------------------------------------------------------------------------------------- ###
163### COMPUTATIONS
164
165# Compute the SRP acceleration 
166accel = srp.compute(epochs, n_cores = n_cores) * 1e3
167        
168log.disabled = False
169
170# Always unload the SPICE kernels
171sp.unload(METAKR)
172
173### ... Elapsed time
174toc = timeit.default_timer()
175time_min = int(floor((toc-tic)/60))
176time_sec = int(mod((toc-tic), 60))
177print("")
178print("\t Elapsed time: %d min, %d sec" %(time_min, time_sec))
179print("")
180
181### ------------------------------------------------------------------------------------------------------- ###
182### PLOT
183
184epochs  = [float( epc - epc_et0)/3600 for epc in epochs]
185
186fig, ax = plt.subplots(3, 1, figsize=(14,8), sharex = True)
187
188ax[0].plot(epochs, accel[:,0], linewidth = 2, color = "tab:blue")
189ax[0].set_ylabel('X [m/s^2]')
190ax[1].plot(epochs, accel[:,1], linewidth = 2, color = "tab:blue")
191ax[1].set_ylabel('Y [m/s^2]')
192ax[2].plot(epochs, accel[:,2], linewidth = 2, color = "tab:blue")
193ax[2].set_ylabel('Z [m/s^2]')
194ax[2].set_xlabel('Hours from t0')
195
196plt.tight_layout()
197plt.show()
198
199### ------------------------------------------------------------------------------------------------------- ###
200
201