Doing the Mars run with fusion propulsion at 1 G

We received a comment on LinkedIn about how fast the “Mars run” could be achieved with a sustained 1 G acceleration. The reader suggested this could be done in 40 hours. What engine parameters would be required to make that happen?

Using a simple constant-acceleration, straight-line analysis, you can indeed compute that the trip should take only a couple of days. Assuming a Mars conjunction, the straight distance is about 0.5 AU. At this speed you can ignore the gravitational effects of the sun and so the distance is a simple integral of the acceleration: d = 1/2 at2. The ship accelerates for half the time then decelerates, and the change in velocity is ΔV = at. Combining the two halves of the trip, at an acceleration of 9.8 m/s2, the trip takes about 2.1 days.

% straight line: distance s = 0.5*at^2
acc = 9.8;             % accel, m/s^2
aU  = Constant('au');  % km
dF  = 0.5*aU*1000;     % distance, m
t   = sqrt(4*dF/acc);  % time for dF, s
dV  = t*acc/1000; % km/s
fprintf('\nAccel: %g m/s^2\n',acc)
fprintf('Time: %g days\n',t/86400)
fprintf('Delta-V: %g km/s\n',dV)

Accel: 9.8 m/s^2
Time: 2.02232 days
Delta-V: 1712.34 km/s

Now, your ship mass includes your payload, your engine, your fuel tanks and your fuel. Assume we want to move a payload of 50,000 kg, somewhat larger than the NASA Deep Space Habitat. The engine mass is computed using a parameter called the specific power, in units of W/kg. The fuel tank mass is scaled from the fuel mass, typically adding another 10%. When we run the numbers, we find that the engine needs to have a specific power of about 1×108 W/kg, and an exhaust velocity of about 5000 km/s results in the maximum payload fraction. We can compute the fuel mass and trajectory using our MassFuelElectricConstantUE and StraightLineConstantAccel toolbox functions:

d = StraightLineDataStructure;
d.dF = 0.5*aU;   % 0.5 AU
d.tF = 2*86400;  % 2.1 days
d.uE = 5000;     % km/s
d.eta = 1;       % jet power sigma
d.sigma = 1e8;   % W/kg
d.mP = 50000;    % kg
dOut = MassFuelElectricConstantUE( d )
StraightLineConstantAccel(dOut)
fprintf('Acceleration: %g m/s^2\n',dOut.a*1e3)
fprintf('Exhaust velocity: %g km/s\n',d.uE)
fprintf('Initial power: %g GW\n',dOut.p(1)*1e-9)
fprintf('Engine mass: %g kg\n',dOut.mE)
fprintf('Payload fraction: %g\n',dOut.mP/dOut.m0)

Acceleration: 10.02 m/s^2
Exhaust velocity: 5000 km/s
Specific Power: 1e+08 W/kg
Initial power: 2832.61 GW
Engine mass: 28326.1 kg
Payload fraction: 0.442172

This produces the following plots:

The power needed is… over 2.8 terawatts! That’s about equal to the total power output of the entire Earth, which had an installed power capacity of 2.8 terawatts in 2020. And the engine would need to weigh less than 30 tons, about the size of a loaded tractor-trailer truck. For comparison, we estimate a Direct Fusion Drive would produce about 1 MW per ton, which is a specific power of 1×103 W/kg. So, this is why you see us trying to design an engine that can do the Mars transfer in 90 days and not 3 days!

Now, there is another consideration here. Namely, constant acceleration at 1 G is not the optimal solution by any means. The optimal solution for a fast, light transfer is actually a linear acceleration profile. This knowledge goes way back: 1961! Here’s a reference:

Leitmann, George. "Minimum Transfer Time for a Power-Limited Rocket." Journal of Applied Mechanics 28, no. 2 (June 1, 1961): 171-78. https://doi.org/10.1115/1.3641648.

This would mean that the engine changes its exhaust velocity during trip, passing through infinity at the switch point. We compute this in our “straight-line, power-limited” or SLPL function series. While this can’t be done physically, even an approximation of this with a variable impulse thruster will one day be more efficient than constant acceleration or thrust. How much better? The power needed is nearly 1/2 the constant acceleration solution, 1.5 TW, and the specific power needed is reduced by half, to 5.6×107 W/kg. However, those are still insane numbers!

mD = 80000;  % dry mass: engine, tanks, payload
m0 = 1.5*mD; % wet mass: with fuel
tF = 3*86400;
vF = 0;

[Pj,A,tau] = SLPLFindPower( aU, tF, vF, mD, m0 );

mTank = 0.05*(m0-mD); % tanks, scale with fuel
mLeft = mD-mTank;
mEngine = mLeft - mPayload;

disp('Straight-line Power-limited (linear accel)')
fprintf('Engine power is %g GW\n',Pj*1e-9);
fprintf('Engine mass is %g kg\n',mEngine);
fprintf('Payload mass is %g kg\n',mPayload);
fprintf('sigma is %g W/kg\n',Pj/mEngine);

SLPLTrajectory( A, tau, Pj, m0, tF )

Straight-line Power-limited (linear accel)
Engine power is 1573.26 GW
Engine mass is 28000 kg
Payload fraction is 0.416667
sigma is 5.6188e+07 W/kg

The trajectory and engine output are plotted below. The linear acceleration results in a curved velocity plot, while in the constant acceleration case, we saw a linear velocity plot. You can see the spike in exhaust velocity at the switch point, which occurs exactly at the halfway point.

After all, who needs 1G gravity when the trip only takes 2 days?

Foe even more fun though, we computed a planar trajectory to Mars using the parameters we found – just to confirm the straight-line analysis is in fact a good approximation. This figure shows the paths the optimization takes:

Earth to Mars Trajectory, 2.1 days, 0.5 AU traversed

It is in fact approximately a straight line!

In reality though, these power system numbers are not even remotely plausible with any technology we are aware of today. That’s why we are designing engines to reduce the Mars trip time to 90 days from 8 or 9 months – still a big improvement!

New Fusion Reactor Design Function

The Fusion Energy Toolbox for MATLAB is a toolbox for designing fusion reactors and for studying plasma physics. It includes a wide variety of physics and engineering tools. The latest addition to this toolbox is a new function for designing tokamaks, based on the paper in reference [1]. Tokamaks have been the leading magnetic confinement devices investigated in the pursuit of fusion net energy gain. Well-known tokamaks that either have ongoing experiments or are under development include JET, ITER, DIII-D, KSTAR, EAST, and Commonwealth Fusion Systems’ SPARC. The new capability of our toolboxes to conduct trade studies on tokamaks allows our customers to take part in this exciting field of fusion reactor design and development.

The Fusion Reactor Design function checks that the reactor satisfies key operational constraints for tokamaks. These operational constraints result from the plasma physics of the fusion reactor, where there are requirements for the plasma to remain stable (e.g., not crash into the walls) and to maintain enough electric current to help sustain itself. The tunable parameters include: the plasma minor radius ‘a’ (see figure below), the H-mode enhancement factor ‘H’, the maximum magnetic field at the coils ‘B_max’, the electric power output of the reactor ‘P_E’, and the neutron wall loading ‘P_W’, which are all essential variables to tokamak design and operation. H-mode is the high confinement mode used in many machines.

Illustration of the toroidal plasma of a tokamak. R is the major radius while a is the minor radius of the plasma. The red line represents a magnetic field line which helically winds along the torus. Image from [2].

This function captures all figure and table results in the original paper. We implemented a numerical solver which allows the user to choose a variable over which to perform a parameter sweep. A ‘mode’ option has been incorporated which allows one to select a desired parameter sweep variable (‘a’, ‘H’, ‘B_max’, ‘P_E’, or ‘P_W’) when calling the function. Some example outputs of the function are described below.

As an example, we will consider the case of tuning the maximum magnetic field at the coils ‘B_max’. The figure below plots the normalized operation constraint parameters for a tokamak as functions of B_max from 10 Tesla to 25 Tesla. The unshaded region, where the vertical axis is below the value of 1, is the region where operational constraints are met. We see that for magnetic fields below about 17.5 Tesla there is at least one operation constraint that is not met, while for higher magnetic fields all operation constraints are satisfied, thus meeting the conditions for successful operation. This high magnetic field approach is the design approach of Commonwealth Fusion Systems for the reactor they are developing [3].

Operational constraint curves as a function of B_max. Successful operation occurs if all of the curves are in the unshaded region. Note, f_B/f_NC, a ratio of the achievable to required bootstrap current, is set equal to 1. In this case P_E = 1000 MW, P_W = 4 MW/m2, and H = 1. For more details on the plotted parameters and how they function as operational plasma constraints, see reference [1].

Note, however, that there is a material cost associated with achieving higher magnetic fields, as described in reference [1]. This is illustrated in the figure below, which plots the cost parameter (the ratio of engineering components volume V_I to electric power output P_E) against B_max. There is a considerable increase in cost at high magnetic fields due to the need to add material volume that can structurally handle the higher current loads required.

Cost parameter (units of volume in cubic meters per megawatt of power, m3/MW) as a function of B_max.

In this post we illustrated the case of a tunable maximum magnetic field at the coils, though as mentioned earlier, there are other parameters you can tune. This function is part of release 2022.1 of the Fusion Energy Toolbox. Contact us at info@psatellite.com or call us at +01 609 276-9606 for more information.

Thank you to interns Emma Suh and Paige Cromley for their contributions to the development of this function.

[1] Freidberg, Mangiarotti, and Minervini, “Designing a tokamak fusion reactor–How does plasma physics fit in?”, Physics of Plasmas 22, 070901 (2015); https://doi.org/10.1063/1.4923266
[2] http://www-fusion-magnetique.cea.fr/gb/iter/iter02.htm
[3] https://news.mit.edu/2021/MIT-CFS-major-advance-toward-fusion-energy-0908

Hohmann Transfer Simulation with the Spacecraft Control Toolbox

Hohmann transfers are a well-known maneuver used to change the semi-major axis of an orbit. The Spacecraft Control Toolbox allows you to compute the required velocity changes, and integrate them into a full simulation.

In this demonstration, we create a 6U CubeSat that has 3 orthogonal reaction wheels and a single hydrazine thruster. The thruster is aligned with the body x-axis and must be aligned with the velocity vector to do the maneuver. An ideal Hohmann maneuver is done with impulsive burns at two points in the orbit. In reality, with a thruster, we have to do a finite burn.

The Hohmann transfer is computed with the following Spacecraft Control Toolbox code:

rI = [-7000;0;0];
vI = [0;-sqrt(mu/Mag(rI));0];
OrbMnvrHohmann(Mag(rI),rF);
[dV,tOF] = OrbMnvrHohmann(Mag(rI),rF);

The first time OrbMnvrHohmann is called, it generates the plot below of the planned Hohmann transfer. The function computes the delta-V and also the time of flight, which will be used to determine the start time of the second thruster burn.

We create a short script with numerical integration to implement the maneuver using a thruster. The burn durations are computed based on the thrust and the mass of the spacecraft. In this case, they are about three minutes long. The maneuver is quite small, so the mass change is not important. The attitude control system uses the PID3Axis function which is a general-purpose attitude control algorithm. The simulation is a for loop, shown below. The ECI vector for the burn is passed to the attitude control system, which updates every step of the simulation.

% Simulation loop
for k = 1:n

  % Update the controller
  dC.eci_vector = uBurn(:,kMnvr);
  [tRWA, dC]    = PID3Axis( x(7:10), dC );

  % Start the first burn
  inMnvr = false;
  if( t(k) > tStart(1) && t(k) < tEnd(1) )
    inMnvr = true;
  end

  % Switch orientation
  if( t(k) > tEnd(1) )
    kMnvr = 2;
  end

  % Start the second burn
  if( t(k) > tStart(2) && t(k) < tEnd(2) )
    inMnvr = true;
  end
  
  if( inMnvr )
    dRHS.force = thrustE*QTForm(x(7:10),dC.body_vector)*nToKN; % kN
  else
    dRHS.force = [0;0;0];
  end
  el = RV2El(x(1:3),x(4:6));
  xP(:,k) = [x;tRWA;Mag(dRHS.force)/nToKN;el(1);el(5)];

  % Right hand side
  dRHS.torqueRWA = -tRWA;
  x = RK4(@RHSRWAOrbit,x,dT,0,dRHS);
end

The maneuver logic just waits a quarter orbit then performs the first burn, by applying the thrust along the body vector. It then waits for the time of flight and then starts the next burn. The start and stop times are pre-computed. RK4 is Fourth Order Runge-Kutta, a popular numerical algorithm included with the toolbox.

At the final orbit radius an attitude maneuver is needed to reorient for the final burn.

The spacecraft body rates, in the body frame, during the maneuver are shown below.

The reaction wheel rates are shown below. The simulation does not model any particular wheel. Friction is not included in the simulation, although the right-hand-side function can include friction.

The wheel torques and rocket thrust are shown below. The thruster is a 0.2 lbf hydrazine thruster that is based on the Aerojet-Rocketdyne MR-103. The PID controller does not demand much torque.

The semi-major axis and eccentricity are shown below. The middle portion is during the transfer orbit.

The eccentricity is zero at the start and finish. Note the slope in both eccentricity and semi-major axis due to the finite acceleration. At the end of the simulation, we print the achieved orbital elements:

Final SMA        7099.72 km
  SMA error         0.28 km
Final e          1.3e-05

The result is very close to the ideal solution!

This post shows how you can easily integrate attitude and orbit control. Email us for more information! We’d be happy to share the script. We can also offer a 30 day demo to let you explore the software.

PSS Toolboxes 2021.1 Now Available!

Version 2021.1 of Princeton Satellite Systems toolboxes for MATLAB is now available! Over 50 new functions and scripts are included. Many other existing functions have been improved.

One new function is AtmNRLMSISE.m, an atmosphere function based on the NRL MSISE model. It is uses extensive flight data and includes sun effects. It computes the overall density and the number density of all atmospheric constituents. Our function has an easy to use interface that automatically incorporates the sun information and lets you input your spacecrafts ECI coordinates. You can also choose to use the original interface. Here is a comparison with the existing scale height model.

We provide a complete set of functions for planning lunar missions in the Missions module. The software includes landing control systems and trajectory optimizaton tools. You can use our Optical Navigation system for your cis-lunar missions and explore our cutting-edge neural network terminal descent software.

Here are two images from an optical navigation simulation for a solar sail.

Solar Sail and Earth paths in the heliocentric frame.
Navigation camera view.

The Spacecraft Control Toolbox provides you with a lot of ways to do things, so you can use your own creativity to perform analyses or design a mission.

Contact us to purchase or for a demo!

Updates for the 2019 Aircraft Control Toolbox

We’ve added some new tools to the Aircraft Control Toolbox for our upcoming 2019 release. The first is a new GUI for creating aircraft models. You import a Wavefront OBJ files and then you point and click to define leading edges, wing areas, engine locations and so forth. This makes it easier to import the geometric data. The GUI is shown below. It illuminates the view that you need to use for a given geometric element in red. The inertia matrix is generated from the mass and the surface geometry.

The new Model Creation GUI

A new simulation function was added to use the data from this GUI. It has a flat Earth aircraft model with a plugins architecture. You can add your own lift, drag and thrust models or use the simple built-in models. It is much simpler than AC.m which is designed to be a comprehensive high-fidelity simulation. We’ve added a new animation GUI to show you the results of your simulations.

We expect 2019.1 to be available in June. You can get a demo with previews of the new functions now.

Lunar Orbit Insertion Maneuver

New functions in the Lunar Cube module in 2016.1 allow you to easily plan lunar insertion and orbit change maneuvers. In the following pictures you can see a lunar orbit insertion from a hyperbolic orbit. In all figures the lunar terrain is exaggerated by a factor of 10.

The same maneuver looking down on the orbit plane. The green arrows are the force vectors.

The following figure shows a two maneuver sequence. The first puts the spacecraft into an elliptical orbit. The second circularizes the orbit.

Lunar Cube Module for 2016.1

We are adding the Lunar Cube Module in 2016.1 to our CubeSat Toolbox for MATLAB! It allows users to analyze and simulateCubeSats in lunar transfer and lunar orbit. It includes a new dynamical model for CubeSats that includes:

  • Earth, Moon and Sun gravity based on the JPL ephemerides
  • Spherical harmonic lunar gravity model
  • Reaction wheels
  • Thrusters
  • Power generation from solar panels
  • Battery energy storage
  • Variable mass due to fuel consumption
  • Solar pressure disturbances
  • Lunar topographic model
  • New graphics functions for lunar orbit operations
  • Lunar targeting function
  • Lunar mission control function for attitude control and orbit control

The module includes a script with a simulation of a 6U Cubesat leaving Earth orbit and reaching the moon. The following figure shows the Earth to Moon trajectory.

LunarTrajectory

This figure shows the transfer orbit near the moon. The lunar topography is exaggerated by a factor of 10 to make it visible. It is based on Clementine measurements.

LunarEncounter

Here are results from the new LunarTargeting function. It finds optimal transfers to lunar orbits. The first shows the transfer path to the Moon’s sphere of influence.

Test21

The next shows the lunar hyperbolic orbit. In this case the transfer is into a high inclination lunar orbit.

Test22

Contact us for more information!

Pluto Spacecraft

Here is a picture of our DFD transfer vehicle. You can see the lander on the front and two Deep Space Optical Communication System (DSOC) assemblies mounted on trusses. There are 2 DFD engines.

PlutoDFD

A picture of the Pluto Lander. The solar panels are illuminated by a laser from the orbiter. The lander has a dry mass of 150 kg.

LunarLander

Both were designed in the Spacecraft Control Toolbox v2015.2.

You can get more information about the Pluto orbital mission on Slideshare.

Fission Powered Lunar Lander

Settlements on the moon, for mining and scientific research, will require routine travel between lunar orbit and the lunar surface. One idea is to use a lunar shuttle with a nuclear fission rocket engine. The hydrogen fuel would come from water on the moon. Fission rockets have twice the specific impulse of the best chemical rockets leading to low fuel consumption. In addition, they would leave the oxygen from the electrolysis of water available for the lunar settlements.

Stanley Borowski of NASA/GRC is co-author of a paper giving the status of nuclear fission rockets:

NTR Technology Development and Key Activities Supporting a Human Mars Mission in the Early-2030 Timeframe

Fission rockets were developed in the 1970’s but the technology was never tested in flight. We used his paper to create a fission rocket. A 3D model based on a drawing the paper is shown below:

NTP

We built the launch vehicle using a single script in the Spacecraft Control Toolbox for MATLAB:

Spacecraft Control Toolbox 2015.1

The script uses a bilinear tangent steering law to estimate the required two way delta-v. The lander flies to 12 km where it meets a freighter. The crew is housed in an Orion spacecraft. The vehicle is shown below:

FissionNL

The landing legs are based on the Apollo Lunar Module. The liquid hydrogen is stored in the 4 spherical tanks. The nuclear thermal engine is hidden by the box to which the legs are attached. The lander lifts the Orion spacecraft and 6000 kg more of payload which would include helium-3 mined on the moon.

The Orion model was created by Amazing3DGraphics. Amazing3D is really good at creating low polygon count models that are useful for simulation and disturbance modeling.

The script and new supporting functions will be available as part of SCT Release v2015.2.

Comparing disturbance models

A customer recently asked us for help comparing the disturbance analysis available in the CubeSat Toolbox with the full model in SCT Professional, for a 3U CubeSat. That is, to compare CubeSatDisturbanceAnalysis to Disturbances. The CubeSat Toolbox uses a simplified model of the spacecraft geometry as a single set of fixed areas, nominally for a rectangular prism. The full model in SCT Pro allows for articulated and rotating bodies built of individual components. The CubeSat Toolbox has a subset of the environment and disturbance functions in Pro but includes drag and optical disturbances in Earth orbit. Given enough options for the two models it should be possible to get the exact same results. We will sketch out our process and discoveries in this post.

Creating a model of a 3U CubeSat in Pro was the easy part, as it is just a single box component, and verifying that the areas match those from CubeSatFaces is accomplished by a trivial command-line printout. Comparing the results of the optical and drag analysis is much more complex as there are so many variables:

  • Atmosphere model: exponential, scale height, J70
  • Solar flux model
  • Earth albedo model
  • Earth radiation model
  • Satellite optical properties (specular, diffuse, absorptive)
  • Attitude pointing (LVLH vs. ECI)

In order to compare the results, we had to call both disturbance models in a script and generate plots overlaying the resulting forces and torques. Here is the code defining the model for SCT Pro.

m = CreateComponent( 'make', 'box', 'x', 0.1, 'y', 0.1, 'z', 0.3,...
'name', 'Core', 'body', 1, 'mass', mass, ...
'faceColor', 'gold foil', 'emissivity', thermal.emissivity,...
'absorptivity', thermal.absorptivity, 'sigmaT', optical.sigmaT,...
'sigmaA', optical.sigmaA, 'sigmaD', optical.sigmaD, 'sigmaS', optical.sigmaS,...
'inside', 0);
BuildCADModel( 'add component', m );

You can see the thermal and optical properties that must be specified as well as the mass and dimensions. The spacecraft is inertially fixed and put into an equatorial orbit, so we would expect zero drag along the z axis and the x/y forces to oscillate at orbit rate. Then to call the disturbance model we generate a low Earth orbit, get the Earth environment and run the analysis, like so:

[r,v] = RVOrbGen(el,t);
e = EarthEnvironment( r, v, jD, d );
hD = Disturbances( 'init', g, e );
[fD,tD] = Disturbances( 'run', g, e, hD );

The EarthEnvironment function is where the guts of the space environment modeling occurs. This includes specifying albedo and radiation constants, calculating the atmospheric density over the orbit, computing the sun vector and solar flux magnitude, checking for eclipses, computing the Earth magnetic field, and correcting the inertial velocity for the rotation of the atmosphere for drag calculations. In the CubeSat toolbox, this data is computed inside the dynamical model, RHSCubeSat. The same steps of creating the model and calling the disturbance function are shown below.

c = RHSCubeSat;
c.mass = 3;
c.inertia = InertiaCubeSat( '3U', c.mass );
[a,n,rho] = CubeSatFaces('3U',1);
c.surfData.nFace = n;
c.surfData.area = a;
c.surfData.rFace = rho;
for k = 1:6
% Radiation coefficients [absorbed; specular; diffuse]
c.surfData.sigma(:,k) = [optical.sigmaA;optical.sigmaS;optical.sigmaD];
end
c.atm = [];
q = QZero*ones(1,size(r,2));
[tC, fC, h, hECI, fr, tq] = CubeSatDisturbanceAnalysis( c, q, r, v, jD );

Let’s look at drag first, as it proved to be the easiest to verify. The primary difference between the CubeSat model and full disturbance model for drag initially was the atmosphere model itself: CubeSat uses the Jacchia 1970 model by default, while EarthEnvironment specifies a scale height atmosphere. The Jacchia 1970 model accounts for changes in density with the time of day, resulting in an orbit rate periodicity; however it is computationally more intensive and not needed in preliminary analysis. The scale heights model depends only on altitude and is very quick. The CubeSat dynamic model already had an option to switch to the scale height atmosphere if desired, so we added that same option to the CubeSat disturbance analysis function. This promptly resulted in a close result between the models for the drag force.

Comparison of the drag force between the CubeSat and full disturbance models for a 3U CubeSat

Comparison of the drag force between the CubeSat Toolbox and SCT Pro disturbance models for a 3U CubeSat

A slight variation remains due to a difference in the transformation between the inertial frame and Earth-fixed frame between the two models. This transformation is used to account for the rotation of Earth’s atmosphere, as drag depends on the relative velocity between the surface and the air. CubeSat uses a fast almanac function, ECIToEF, to compute this matrix for a given Julian date. This model accounts for nutation but not as accurately as TruEarth does in Pro. The EarthEnvironment function in Pro, however, uses a simpler transformation using Earth’s rotational rate about the inertial z-axis, ignoring nutation. This accounts for the nonzero Z force in the CubeSat output, which can be seen to be four orders of magnitude less than the X/Y forces. Both approaches are valid for a preliminary analysis so we accept this small remaining difference.

Producing an equally close comparison for the optical forces unearthed a few bugs in the CubeSat version as well as differences in fidelity that are intentional. First, recall that there are three main contributions to optical forces in Earth orbit: solar flux; Earth albedo – that is, reflected flux; and Earth infrared radiation. The fluxes can be modeled simply as constants, or at higher fidelity by accounting for distance from the radiating body. The full disturbance model accounts for the change in solar flux over the year as the Earth moves in its orbit, which amounts to about 100 W/m2 or 7% of the average flux. The CubeSat environment model was not doing this, but since it was already calling the sun vector function which calculates the needed data, we decided to add it. The sun vector itself can be modeled a number of ways, with CubeSat providing a low fidelity almanac version and Pro a higher fidelity almanac option as well as an option for JPL ephemerides.

Making a temporary change in CubeSat to use the higher fidelity sun almanac provided closer results, but there were still differences in the optical forces. A check on the optical coefficients revealed that Disturbances assumed 100% diffuse reflection for planetary infrared radiation while the CubeSat version assumed 100% absorption. This was result of a misunderstanding of the model when the CubeSat version was created. The intention of the model is to assume 100% absorption, but the radiation has to be reemitted or the temperature of the body would increase to infinity. Hence the diffuse coefficient in the Pro model. So, the CubeSat version’s coefficients were corrected to match. This improved the agreement between the models but still not completely.

Further investigation uncovered a bug in the calculation of the planetary radiation flux and Earth albedo flux. We call this a copy/paste bug, as the distance scaling factor was in the code, but applied to the wrong variable when the CubeSat version was created from the Pro version – so the scaling was unitized out. Correctly applying the scale factor to the flux input (and not its unit vector) to the base solar force function, SolarF, resulted in exact agreement between the two models. The figure below shows the final differences with the CubeSat sun almanac function restored to its default value.

Optical force comparison

Comparison of the optical force between the simplified and full disturbance models for a 3U CubeSat

To summarize, we identified the following subtleties in making a direct comparison between these two disturbance models:

  • Selecting the same atmosphere model, AtmDens2, for both cases
  • Understanding the sun vector model in each case (SunV1 vs. SunV2)
  • Adding the solar flux scaling with the Earth’s distance from the sun to the CubeSat model (1 line)
  • Updating the optical coefficients for planetary radiation in the CubeSat model (to be diffuse)
  • Correcting a bug in the scaling of the Earth albedo and radiation fluxes with altitude in CubeSat
  • Accounting for Earth’s nutation in transforming velocity to the Earth-fixed frame

This provides a great illustration of how careful one needs to be in performing any disturbance analysis, as there are so many subleties in even simple models that must be recorded if results are to be replicated. Every constant, every scale factor, every coefficient, every source of ephemeris needs to be understood.

We were happy to provide the results of our analysis to our customer so that he could decide which models he wants to use in his analysis. This is an example of the technical support we can provide to all our SCT customers; if you have a question, just ask!