Collisional Radiative Model paper Published in Review of Scientific Instruments 

Our new collisional radiative model paper is titled “Evaluation of a collisional radiative model for electron temperature determination in hydrogen plasma” DOI: It is part of the “Proceedings of the 24th Topical Conference on High-Temperature Plasma Diagnostics.”

This paper talks about a collisional-radiative (CR) model that extracts the electron temperature, Te, of hydrogen plasmas from Balmer-line-ratio measurements and is examined for the plasma electron density, ne, and Te ranges of 1010–1015 cm−3 and 5–500 eV, respectively. The first tests of the CR model on the Princeton Field Reversed Configuration-2 (PFRC-2) have been made, including comparisons with other diagnostics. These comparisons are informative as different diagnostics sample different parts of the electron energy distribution function.

Extracted Te(t) for the data for two values of Pc. The shaded region represents the statistical error bar.

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 )
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.

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!

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];
[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;

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

  % Start the second burn
  if( t(k) > tStart(2) && t(k) < tEnd(2) )
    inMnvr = true;
  if( inMnvr )
    dRHS.force = thrustE*QTForm(x(7:10),dC.body_vector)*nToKN; % kN
    dRHS.force = [0;0;0];
  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);

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.

A Third Planet Discovered Orbiting Proxima Centauri


A third planet, as large as 26% of the mass of Earth, has been discovered orbiting our nearest stellar neighbor, Proxima Centauri .Astronomer João Faria and his collaborators detected Proxima Centauri d using the Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations.

It would be exciting to send a spacecraft to enter the Alpha-Centauri system and orbit this planet. At Princeton Satellite System we’ve looked at interstellar flight using the Direct Fusion Drive nuclear fusion propulsion system.

Interstellar Fusion Propulsion

At the 2021 Breakthrough Energy Conference we presented findings for both flyby and orbital missions. Flyby missions are easier, but orbit entry would allow detailed study of the planet. A flyby gets your spacecraft close, but it is moving really fast!

The following charts give an outline of our talk. The first shows the optimal exhaust velocity based on sigma, the ratio of power to mass. Our designs have a sigma from 0.75 to 2 kW/kg. With 2 kW/kg, the optimal exhaust velocity is 4000 km/s. The mission would take about 800 years. Our current designs can’t get exhaust velocities higher than 200 km/s. We’d need another method to produce thrust.

Mission Analysis

The next plot shows a point mission that reaches Alpha Centauri in 500 years. This requires a sigma of about 20. The spacecraft accelerates and decelerates continuously. The mission could be improved by staging, much like on a rocket that launches from the Earth into orbit.

Selected Mission

The next figure shows how the starship would enter the Alpha Centauri system.

Alpha Centauri System Insertion

The final plot shows the orbital maneuvers that lower the orbit and rendezvous with the planet.

Lowering the orbit to rendezvous with the planet.

Even 500 years is a long time! This is over ten times the lifetime of Voyager, but much less than some engineering marvels built on the Earth.

We hope to someday be able to build fusion powered spacecraft that will head into interstellar space!

Moonfall – The Movie

Moonfall is a movie coming out in 2022. It creates a scenario where the Moon’s orbit is changed and set on a collision course with the Earth. It is fun to work out the orbital mechanics.

Let us assume that the Moon is in a circular orbit around the Earth. It is actually more influenced by the Sun than the Earth, but the circular orbit approximation is sufficient for our purposes. A mysterious force changes the orbit from circular to elliptical so that at closest approach it hits the Earth. The transfer orbit has an eccentricity of 0.9673 and a semi-major axis of 195000 km. The new orbital period is 9.9 days so it will hit the Earth in 5 days!

What kind of force is needed? The required velocity change is 0.83 km/s so a force of 6 x 1016 N applied over 10 seconds is required. To get an idea of how large that force really is, the Space Launch System (SLS) Block 2 vehicle produces about 10 million pounds of thrust [1], which is approximately 50 x 106 N (50 MN). Hence it would take 1.2 billion SLS rockets firing for 10 seconds to perform such a re-direction of the Moon! An image of the SLS is shown below (image from [1]).

As the Moon approaches the Earth it is going to raise the tides. A simple formula (really only valid when the Moon is far from the Earth) is

where is the gravitational constant for the moon, is the gravitational constant for the Earth, r is the distance between the Earth and Moon and a is the radius of the Earth. The distance during the approach and the wave height are shown in the following plot.

By around 3 days the tides started getting really big! We’d expect the Moon’s gravitational force also to pull on the solid part of the Earth’s surface, causing all sorts of trouble.



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!

The Optical Navigation Module for the Spacecraft Control Toolbox is Now Available

Space optical navigation employs a camera for attitude determination and a second high dynamic range camera on a pan/track mount for terrain and celestial body tracking. Navigation and attitude determination are performed in a Bayesian framework using anUnscented Kalman Filter with an IMU as the navigation and attitude base. The Optical Navigation Module provides MATLAB code for implementing optical navigation. Additional measurements can be added including a sun sensor for sun distance measurements in interplanetary space, Global Positioning System (GPS) measurements near the Earth, and range and range rate from ground stations or other spacecraft in deep space. The system is suitable for both lunar and Mars landing missions and icy moon and asteroid orbital missions such as Artemis, Lunar Orbital Platform Gateway, Orion Multi-Purpose Crew Vehicle, Europa Clipper, Lucy, Psyche. It is also applicable to any situation where GPS is not available.

The Optical Navigation Module allows you to implement an optical navigation system for any of these applications. It includes dynamical models for cis-lunar and deep space missions along with measurement models for all of these sensors. Several scripts provide examples to get you going quickly.

This picture shows the camera aimed at the horizon and the stars that it can see during Earth reentry. The step counter gives the integration step. The star numbers are sequential from the file of stars but the stars come from the Hipparcos catalog.

This pictures shows the laboratory hardware for an optical navigation camera on a pan/tilt mount. Flexible cables eliminate the need for slip rings simplifying the design. The platform is driven by orthogonal stepping motors with harmonic drives.

Note the size. As with all of our toolboxes, full source code is provided.

What Makes a Reaction Wheel a Reaction Wheel?

A DC motor is the core of all momentum and reaction wheels. If you apply a voltage a, current will be produced which will cause the wheel to change speed. At the same time, the back electromotive force (EMF) will build up, eventually driving the motor torque to zero.

This is evident from the dynamical equation for a DC motor.

J\dot{\omega} = \frac{K_T}{R}\left(V - K_T\omega\right) + T_F

J is the inertia, K_T is the torque constant, V the voltage, T_F the friction torque, R the motor impedance and \omega is the angular rate of the shaft.

You can turn this into a reaction wheel by adding current feedback as shown in the following block diagram.

G is the forward gain. The input is the desired torque. This is divided by the torque constant to get the desired current. The difference between the motor current and the desired current is integrated. How do you pick the gain? If you work through the equations you will get this equation for the voltage, V

 \dot{V} + \frac{G}{R}V = G\frac{T_C}{K_T} + \frac{G}{R}K_T\omega

R/G is the time constant. The response is shown in the following plot. Even as the speed increases, the difference between the desired torque and motor torque is nearly zero.

Visiting Planet 9

In 2015, astronomers from Caltech determined that a giant ninth planet may be orbiting the Sun. It was called Planet X and then Planet 9. The discovery was based on perturbations in the orbits of TNOs, trans Neptunian Objects. The planet has about the mass of Neptune and is in a 10,000 to 20,000 year solar orbit. Jakub Scholtz of Durham University and James Unwin of University of Illinois at Chicago hypothesize that Planet 9 might be a black hole. The orbit of Planet 9 looks something like this.

We used a semi-major axis of 700 AU, an inclination of 30 degrees and an eccentricity of 0.6. The plot shows the full orbit of Planet 9, but the simulation only shows 150 years of the other planets.

It would be very interesting to visit Planet 9. One way is to use a solar sail. The sail would start on a trajectory aiming at perigee very close to the sun and then accelerate at high speed. Another approach is to use a spacecraft propelled by Direct Fusion Drive, a fusion propulsion system we’ve been working on for several years. A 26000 kg spacecraft with a 12 MW engine and 2000 kg of payload could rendezvous with Planet 9 (based on the above orbit) in just 11 years. This is the spacecraft trajectory

Direct Fusion Drive is based on the Princeton Field Reversed Configuration reactor invented by Dr. Samuel Cohen of the Princeton Plasma Physics Laboratory (PPPL). We have a experiment running at PPPL, funded by an ARPA-E OPEN grant, to perform critical ion-heating tests. Earlier work was funded by the NASA NIAC program. Hopefully we will be in a position to send a mission to Planet 9 in the not too distant future!

This analysis was done using the Spacecraft Control Toolbox v2020.1. Contact us for more information!

Lunar Helium-3 Return to the Earth

Helium-3 is available in the regolith of the moon and is a possible fuel for advanced nuclear fusion reactors on Earth. It would be extracted from the lunar regolith, packaged and returned to Earth. One question is how to return the helium-3 to the Earth. One approach is to use aerodynamic braking to return the helium-3 to a low Earth orbit where it would be picked up by the Space Rapid Transit (SRT) reusable launch vehicle and delivered to an airport where it would be shipped to power plants. SRT It is a two stage to orbit vehicle with a hypersonic air-breathing engine in the first stage.

The overall architecture is shown below.

One of the major advantages of SRT is that it can land and takeoff at any major airport. The first stage can be used as a transport vehicle. Since it is fully reusable and operates like an aircraft it is potentially much less expensive than vertical launch.

The return from the Earth involves launching the helium-3 tanker into orbit and then doing a departure burn that puts the spacecraft in an elliptical Earth orbit with a low perigee. As the return vehicle passes through perigee, aerodynamic drag lowers apogee until apogee and perigee are the same. This is shown in the following plots.

The first plot show the altitude from the Earth, the velocity magnitude and the drag force magnitude. The second plot shows the orbit. The last plot shows how apogee is reduced with each pass through perigee. It takes 10 weeks to enter the final orbit if the orbit perigee is 100 km. Note that perigee doesn’t change. The simulation uses a free-molecular aerodynamic flow model. For simplicity, it does not include lunar gravity perturbations.

Ideally, the lunar return vehicle would be brought back to Earth and reused.

The maneuver uses only drag. A lifting vehicle would have an additional degree of freedom since the force vector could be controlled.

This analysis was done with the Spacecraft Control Toolbox. The function will be available in Version 2020.2 available in early fall. Contact us for more information!