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