Skip to content

Commit d486ee8

Browse files
leahkinerleahkiner
committed
[#1332] Add module unit test
1 parent ecdfef2 commit d486ee8

1 file changed

Lines changed: 322 additions & 0 deletions

File tree

‎src/simulation/dynamics/facetedSRPEffector/_UnitTest/test_facetedSRPEffector.py‎

Lines changed: 322 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,322 @@
1+
2+
# ISC License
3+
#
4+
# Copyright (c) 2026, Autonomous Vehicle Systems Lab, University of Colorado at Boulder
5+
#
6+
# Permission to use, copy, modify, and/or distribute this software for any
7+
# purpose with or without fee is hereby granted, provided that the above
8+
# copyright notice and this permission notice appear in all copies.
9+
#
10+
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11+
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12+
# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13+
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14+
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15+
# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16+
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17+
18+
#
19+
# Unit Test Script
20+
# Module Name: facetedSRPEffector
21+
# Author: Leah Kiner
22+
#
23+
24+
import numpy as np
25+
import pytest
26+
from Basilisk.utilities import SimulationBaseClass
27+
from Basilisk.utilities import simIncludeGravBody
28+
from Basilisk.utilities import macros
29+
from Basilisk.utilities import RigidBodyKinematics as rbk
30+
from Basilisk.simulation import facetedSpacecraftProjectedArea
31+
from Basilisk.simulation import facetedSRPEffector
32+
from Basilisk.simulation import spacecraft
33+
from Basilisk.architecture import messaging
34+
35+
# Required constants
36+
speed_light = 299792458.0 # [m/s] Speed of light
37+
ast_u = 149597870700.0 # [m] Astronomical unit
38+
solar_rad_flux = 1368.0 # [W/m^2] Solar radiation flux at 1 AU
39+
40+
@pytest.mark.parametrize(
41+
("r_BN_N_init", "r_SN_N_init"), # ([m], [m])
42+
[
43+
(np.array([-4020338.690396649, 7490566.741852513, 5248299.211589362]), np.array([0.0, 0.0, 0.0])), # normal
44+
(np.array([7490566.741852513, 5248299.211589362, -4020338.690396649]), np.array([0.0, 0.0, 0.0])), # normal
45+
(np.array([-4020338.690396649, 7490566.741852513, 5248299.211589362]), np.array([1.0e11, 2.0e10, -3.0e10])), # nonzero Sun position
46+
(np.array([2.0e11, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])) # farther separation
47+
]
48+
)
49+
@pytest.mark.parametrize("sigma_BN_init", # [-]
50+
[
51+
np.array([0.0, 0.0, 0.0]),
52+
np.array([0.41421356, 0.0, 0.0]), # +90 deg about x
53+
np.array([0.0, 0.41421356, 0.0]), # +90 deg about y
54+
np.array([0.1, -0.2, 0.05]), # arbitrary small rotation
55+
]
56+
)
57+
@pytest.mark.parametrize("omega_BN_B_init", # [rad/s]
58+
[
59+
np.array([0.0, 0.0, 0.0]),
60+
np.array([0.05, 0.0, 0.0]),
61+
np.array([0.0, 0.05, 0.0]),
62+
np.array([0.0, 0.0, 0.05]),
63+
np.array([0.05, -0.05, -0.05]),
64+
]
65+
)
66+
def test_facetedSRPEffector(show_plots,
67+
r_BN_N_init,
68+
r_SN_N_init,
69+
sigma_BN_init,
70+
omega_BN_B_init):
71+
r"""
72+
**Verification Test Description**
73+
74+
This unit test verifies that the faceted spacecraft solar radiation pressure (SRP) module correctly computes
75+
the aggregate force and torque acting on the spacecraft due to impinging photons from the Sun.
76+
The inertial Sun state information must be subscribed to the module ``SpicePlanetStateMsgPayload`` input message.
77+
The module assumes the spacecraft is modeled as a collection of facets, where the facet geometry information is
78+
passed as a vector of ``FacetElementBodyMsgPayload`` input messages to the module. The facet geometry information
79+
is required to be provided in the spacecraft body frame. This SRP module does not make any assumptions regarding
80+
whether the facets are rigid or articulate. The ``facetedSpacecraftModel`` module can be connected upstream and
81+
used to transform the facet geometry data from the facet frames to the spacecraft body frame. This upstream module
82+
can also be used to configure articulating facets. Finally, this SRP module also requires the sunlit area of all
83+
facets to be pre-computed and connected to the module vector of ``ProjectedAreaMsgPayload`` input messages.
84+
This information can be passed to the SRP module by connecting the facet geometry information to the
85+
``facetedSpacecraftProjectedArea`` module upstream.
86+
87+
This test sets up a simulation with a faceted spacecraft modeled as a cubic hub with two attached circular solar
88+
arrays. Six square facets represent the cubic hub and four circular facets represent the two solar arrays.
89+
The test varies the initial state information of both the spacecraft and the Sun. The test checks that the computed
90+
SRP forces and torques at each timestep match the values output from the module.
91+
92+
**Test Parameters**
93+
94+
Args:
95+
r_BN_N_init (np.ndarray): [m] shape (3,) Initial spacecraft inertial position
96+
r_SN_N_init (np.ndarray): [m] shape (3,) Initial Sun inertial position
97+
sigma_BN_init (np.ndarray): [-] shape (3,) Initial spacecraft inertial attitude
98+
omega_BN_B_init (np.ndarray): [rad/s] shape (3,) Initial spacecraft inertial angular velocity
99+
100+
**Description of Variables Being Tested**
101+
102+
The test checks that the module correctly computes the aggregate srp force and torque acting on the spacecraft
103+
body frame at each time step. The values logged from the module are checked to match the computed truth values.
104+
"""
105+
106+
task_name = "testTask"
107+
process_name = "testProcess"
108+
test_sim = SimulationBaseClass.SimBaseClass()
109+
test_time_step_sec = 0.001 # [s]
110+
test_process_rate = macros.sec2nano(test_time_step_sec)
111+
test_process = test_sim.CreateNewProcess(process_name)
112+
test_process.addTask(test_sim.CreateNewTask(task_name, test_process_rate))
113+
114+
# Create the Sun spice inertial state message
115+
sun_state_message_data = messaging.SpicePlanetStateMsgPayload()
116+
sun_state_message_data.PositionVector = r_SN_N_init # [m]
117+
sun_state_message_data.VelocityVector = [0.0, 0.0, 0.0] # [m/s]
118+
sun_state_message = messaging.SpicePlanetStateMsg().write(sun_state_message_data)
119+
120+
# Create the Sun
121+
grav_factory = simIncludeGravBody.gravBodyFactory()
122+
sun = grav_factory.createSun()
123+
sun.isCentralBody = True
124+
grav_factory.gravBodies['sun'].planetBodyInMsg.subscribeTo(sun_state_message)
125+
126+
# Create the spacecraft hub
127+
sc_object = spacecraft.Spacecraft()
128+
sc_object.ModelTag = "scObject"
129+
sc_object.hub.mHub = 750.0 # [kg]
130+
sc_object.hub.r_BcB_B = [[0.0], [0.0], [1.0]] # [m]
131+
sc_object.hub.IHubPntBc_B = [[900.0, 0.0, 0.0], [0.0, 800.0, 0.0], [0.0, 0.0, 600.0]] # [kg m^2]
132+
sc_object.hub.r_CN_NInit = r_BN_N_init # [m]
133+
sc_object.hub.v_CN_NInit = [[-5199.77710904224], [-3436.681645356935], [1041.576797498721]] # [m/s]
134+
sc_object.hub.sigma_BNInit = sigma_BN_init # [-]
135+
sc_object.hub.omega_BN_BInit = omega_BN_B_init # [rad/s]
136+
137+
# Facet geometry information
138+
num_facets = 10
139+
area_1 = 1.5 * 1.5 # [m^2]
140+
area_2 = np.pi * (0.5 * 7.5) * (0.5 * 7.5) # [m^2]
141+
facet_area_list = [area_1, area_1, area_1, area_1, area_1, area_1, area_2, area_2, area_2, area_2] # [m^2]
142+
facet_r_CopB_B_list = [np.array([0.75, 0.0, 0.0]),
143+
np.array([0.0, 0.75, 0.0]),
144+
np.array([-0.75, 0.0, 0.0]),
145+
np.array([0.0, -0.75, 0.0]),
146+
np.array([0.0, 0.0, 0.75]),
147+
np.array([0.0, 0.0, -0.75]),
148+
np.array([4.5, 0.0, 0.75]),
149+
np.array([4.5, 0.0, 0.75]),
150+
np.array([-4.5, 0.0, 0.75]),
151+
np.array([-4.5, 0.0, 0.75])] # [m]
152+
facet_nHat_B_list = [np.array([1.0, 0.0, 0.0]),
153+
np.array([0.0, 1.0, 0.0]),
154+
np.array([-1.0, 0.0, 0.0]),
155+
np.array([0.0, -1.0, 0.0]),
156+
np.array([0.0, 0.0, 1.0]),
157+
np.array([0.0, 0.0, -1.0]),
158+
np.array([0.0, 1.0, 0.0]),
159+
np.array([0.0, -1.0, 0.0]),
160+
np.array([0.0, 1.0, 0.0]),
161+
np.array([0.0, -1.0, 0.0])] # [-]
162+
facet_rotHat_B_list = [np.array([0.0, 0.0, 0.0]),
163+
np.array([0.0, 0.0, 0.0]),
164+
np.array([0.0, 0.0, 0.0]),
165+
np.array([0.0, 0.0, 0.0]),
166+
np.array([0.0, 0.0, 0.0]),
167+
np.array([0.0, 0.0, 0.0]),
168+
np.array([0.0, 0.0, 0.0]),
169+
np.array([0.0, 0.0, 0.0]),
170+
np.array([0.0, 0.0, 0.0]),
171+
np.array([0.0, 0.0, 0.0])] # [-]
172+
facet_diffuse_coeff_list = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
173+
facet_specular_coeff_list = np.array([0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9])
174+
175+
# Create the facet element input messages
176+
facet_element_message_list = list()
177+
for idx in range(num_facets):
178+
facet_element_message_data = messaging.FacetElementBodyMsgPayload(
179+
area = facet_area_list[idx],
180+
r_CopB_B = facet_r_CopB_B_list[idx],
181+
nHat_B = facet_nHat_B_list[idx],
182+
rotHat_B = facet_rotHat_B_list[idx],
183+
c_diffuse = facet_diffuse_coeff_list[idx],
184+
c_specular = facet_specular_coeff_list[idx],
185+
)
186+
facet_element_message = messaging.FacetElementBodyMsg().write(facet_element_message_data)
187+
facet_element_message_list.append(facet_element_message)
188+
189+
# Create the faceted spacecraft projected area module
190+
faceted_sc_projected_area = facetedSpacecraftProjectedArea.FacetedSpacecraftProjectedArea()
191+
faceted_sc_projected_area.ModelTag = "facetedSpacecraftProjectedArea"
192+
faceted_sc_projected_area.setNumFacets(num_facets)
193+
for idx in range(num_facets):
194+
faceted_sc_projected_area.facetElementBodyInMsgs[idx].subscribeTo(facet_element_message_list[idx])
195+
faceted_sc_projected_area.sunStateInMsg.subscribeTo(sun_state_message)
196+
faceted_sc_projected_area.spacecraftStateInMsg.subscribeTo(sc_object.scStateOutMsg)
197+
test_sim.AddModelToTask(task_name, faceted_sc_projected_area)
198+
test_sim.AddModelToTask(task_name, sc_object)
199+
200+
# Create the faceted SRP effector module
201+
faceted_srp_effector = facetedSRPEffector.FacetedSRPEffector()
202+
faceted_srp_effector.ModelTag = "facetedSRPEffector"
203+
faceted_srp_effector.setNumFacets(num_facets)
204+
for idx in range(num_facets):
205+
faceted_srp_effector.facetElementBodyInMsgs[idx].subscribeTo(facet_element_message_list[idx])
206+
faceted_srp_effector.facetProjectedAreaInMsgs[idx].subscribeTo(faceted_sc_projected_area.facetProjectedAreaOutMsgs[idx])
207+
faceted_srp_effector.sunStateInMsg.subscribeTo(sun_state_message)
208+
sc_object.addDynamicEffector(faceted_srp_effector)
209+
test_sim.AddModelToTask(task_name, faceted_srp_effector)
210+
211+
# Set up data logging
212+
spacecraft_state_data_log = sc_object.scStateOutMsg.recorder()
213+
sun_state_data_log = grav_factory.gravBodies['sun'].planetBodyInMsg.recorder()
214+
srp_data_log = faceted_srp_effector.logger(["forceExternal_B", "torqueExternalPntB_B"], test_process_rate)
215+
facet_element_projected_area_data_log = []
216+
for outMsg in faceted_sc_projected_area.facetProjectedAreaOutMsgs:
217+
facet_element_projected_area_data_log.append(outMsg.recorder())
218+
test_sim.AddModelToTask(task_name, facet_element_projected_area_data_log[-1])
219+
test_sim.AddModelToTask(task_name, spacecraft_state_data_log)
220+
test_sim.AddModelToTask(task_name, sun_state_data_log)
221+
test_sim.AddModelToTask(task_name, srp_data_log)
222+
223+
# Execute the simulation
224+
test_sim.InitializeSimulation()
225+
sim_time = macros.sec2nano(0.01) # [ns]
226+
test_sim.ConfigureStopTime(sim_time)
227+
test_sim.ExecuteSimulation()
228+
229+
# Retrieve the logged data
230+
timespan = spacecraft_state_data_log.times() * macros.NANO2SEC # [s]
231+
r_BN_N = spacecraft_state_data_log.r_BN_N # [m]
232+
sigma_BN = spacecraft_state_data_log.sigma_BN # [-]
233+
r_SN_N = sun_state_data_log.PositionVector # [m]
234+
srp_force_B_list_sim = srp_data_log.forceExternal_B # [N]
235+
srp_torque_B_list_sim = srp_data_log.torqueExternalPntB_B # [Nm]
236+
facet_element_projected_area_list_sim = []
237+
for data in facet_element_projected_area_data_log:
238+
facet_element_projected_area_list_sim.append(data.area)
239+
240+
# Compute truth data
241+
srp_force_B_list_truth, srp_torque_B_list_truth = compute_srp_force_torque(num_facets,
242+
facet_area_list,
243+
facet_r_CopB_B_list,
244+
facet_nHat_B_list,
245+
facet_diffuse_coeff_list,
246+
facet_specular_coeff_list,
247+
timespan,
248+
sigma_BN,
249+
r_BN_N,
250+
r_SN_N,
251+
facet_element_projected_area_list_sim)
252+
253+
# Check the simulated srp force and torque values match the computed truth values
254+
for idx in range(len(timespan)):
255+
np.testing.assert_allclose(srp_force_B_list_sim[idx],
256+
srp_force_B_list_truth[idx],
257+
atol=1e-7,
258+
verbose=True)
259+
np.testing.assert_allclose(srp_torque_B_list_sim[idx],
260+
srp_torque_B_list_truth[idx],
261+
atol=1e-7,
262+
verbose=True)
263+
264+
def compute_srp_force_torque(num_facets,
265+
facet_area_list,
266+
facet_r_CopB_B_list,
267+
facet_nHat_B_list,
268+
facet_diffuse_coeff_list,
269+
facet_specular_coeff_list,
270+
timespan,
271+
sigma_BN,
272+
r_BN_N,
273+
r_SN_N,
274+
facet_element_projected_area_list_sim):
275+
276+
srp_force_B_list_truth = [] # [N]
277+
srp_torque_B_list_truth = [] # [Nm]
278+
279+
for time_idx in range(len(timespan)):
280+
281+
# Determine unit direction vector pointing from sc to the Sun
282+
sigma_bn = sigma_BN[time_idx] # [-]
283+
dcm_BN = rbk.MRP2C(sigma_bn) # [-]
284+
r_bn_n = r_BN_N[time_idx] # [m]
285+
r_BN_B = np.matmul(dcm_BN, r_bn_n) # [m]
286+
r_sn_n = r_SN_N[time_idx] # [m]
287+
r_SN_B = np.matmul(dcm_BN, r_sn_n) # [m]
288+
r_SB_B = r_SN_B - r_BN_B # [m]
289+
s_hat = r_SB_B / np.linalg.norm(r_SB_B) # [-]
290+
291+
# Determine the SRP pressure at the current sc location
292+
num_au = ast_u / np.linalg.norm(r_SB_B)
293+
srp_pressure = (solar_rad_flux / speed_light) * num_au * num_au # [Pa]
294+
295+
srp_force_B = np.zeros([3,]) # [N]
296+
srp_torque_B = np.zeros([3,]) # [Nm]
297+
for facet_idx in range(num_facets):
298+
299+
projected_area = facet_element_projected_area_list_sim[facet_idx][time_idx] # [m^2]
300+
301+
if projected_area > 0:
302+
cos_theta = projected_area / facet_area_list[facet_idx]
303+
facet_force = -srp_pressure * projected_area * ((1 - facet_specular_coeff_list[facet_idx]) * s_hat
304+
+ 2 * ( (facet_diffuse_coeff_list[facet_idx] / 3)
305+
+ facet_specular_coeff_list[facet_idx] * cos_theta) * facet_nHat_B_list[facet_idx]) # [N]
306+
srp_force_B += facet_force # [N]
307+
srp_torque_B += np.cross(facet_r_CopB_B_list[facet_idx], facet_force) # [Nm]
308+
309+
srp_force_B_list_truth.append(srp_force_B)
310+
srp_torque_B_list_truth.append(srp_torque_B)
311+
312+
return np.array(srp_force_B_list_truth), np.array(srp_torque_B_list_truth)
313+
314+
315+
if __name__=="__main__":
316+
test_facetedSRPEffector(
317+
False, # show plots
318+
np.array([-4020338.690396649, 7490566.741852513, 5248299.211589362]), # [m] r_BN_N_init
319+
np.array([0.0, 0.0, 0.0]), # [m] r_SN_N_init
320+
np.array([0.0, 0.0, 0.0]), # [-] sigma_BN_init
321+
np.array([0.0, 0.0, 0.0]) # [rad/s] omega_BN_B_init
322+
)

0 commit comments

Comments
 (0)