About Stephanie Thomas

Ms. Thomas is vice president of Princeton Satellite Systems. She manages technical programs for the Air Force and participates in corporate management of the company. Ms. Thomas has been with PSS since her first internship as an MIT undergraduate in 1996!

Opening for a part-time bookkeeper/administrator

We are looking for an energetic, results-oriented person able to combine the responsibilities of financial manager and coordinator of back office operations, working part time.

Responsibilities include but are not limited to:

  • financial planning and budgeting;
  • accounting and bookkeeping;
  • customer invoicing and cash management;
  • procurement;
  • HR record keeping and benefits management;
  • coordinating payroll working with an external service provider;
  • managing office facilities.

A candidate must have financial education and a proven record of at least 3 years independently performing on a job with similar responsibilities, possession of skills in QuickBooks or a similar financial package along with MS Office. Excellent communication skills are anticipated.

Please send resumes to info@psatellite.com! Also see our post on LinkedIn.

NIAC Pluto mission talk now available online

On Tuesday, August 23rd I had the privilege of giving my talk on our Fusion-Enabled Pluto Orbiter and Lander at the 2016 NIAC Symposium. The video of the LiveStream is now archived and available for viewing. My talk starts at 17:30 minutes in, after Michael VanWoerkom’s NIMPH talk.

The talk was well-received and we had some good questions from the audience and the LiveStream. In retrospect I did wish I had added a slide on our overall program plan in terms of the PFRC machine and temperature and field strength, since I got quite a few questions on those specifics at the poster session. PFRC-1 demonstrated heating electrons to 0.3 keV in 3 ms pulses. The goal of the current machine – PFRC 2 – is heating ions to 1 keV with a 1.2 kG field. The next machine I refer to in the talk, PFRC 3, would initially heat ions to 5 keV with a 10 kG field, and towards the end of its life we would push the field to 80 kG, heat ions to 50 keV, and add some helium-3 to get actual fusion events. The final goal would be 100 second-duration plasmas with a fusion gain between 0.1 and 2. A completed reactor would operate in steady-state.

Thank you NIAC for this opportunity!!

NIAC Orientation

I had a great time at the NIAC orientation in Washington DC last week, where I got “mugged” with program manager Jason Derleth:

Stephanie Thomas and Jason Derleth posing with a NIAC mug

Stephanie receiving her NIAC mug from Jason

The meeting was at the Museum of the American Indian, which was a great venue with so much beautiful art to see, and a cafe featuring unusual native foods from across America (elderberry sauce on the salmon). I had the opportunity to meet the other NIAC Fellows, and put names and faces to the other creative projects selected, as well as meet the illustrious NIAC external council. These experienced folks provide advice and encouragement throughout the NIAC process from their experience as physicists, engineers, biologists, science hackers, and even science fiction authors.

I have to say, my poster on the fusion rocket engine was popular, and everyone wanted to know how it works, why it hasn’t been funded already, and how soon the engine can be ready. Of course, we have yet to actually demonstrate fusion using Dr. Cohen’s heating method, but that is why we need the NIAC study – to flesh out the science and engineering of the rocket application to help bring in funding for building the next generation machine. And yes, let’s get to Pluto in only 4 years the next time! I’m really looking forward to working on the project in the next few months and presenting it at the NIAC symposium in August!


NASA Innovative Advanced Concepts (NIAC) Selection

We are very pleased to announce that Ms. Stephanie Thomas of Princeton Satellite Systems has been selected to be a 2016 NIAC Fellow. This Phase I study, entitled “Fusion-Enabled Pluto Orbiter and Lander,” will explore the possibility of using Direct Fusion Drive (DFD) to deliver an orbiter to Pluto complete with a lander. DFD is a fusion propulsion concept built upon a small, clean field-reversed configuration fusion reactor with a naturally linear geometry. The reactor becomes a rocket engine when additional propellant flows through, providing power as well as propulsion in one integrated device. This engine could halve the transit time to Pluto to 5 years from the nearly 10 years needed for New Horizons, while delivering 1000 kg worth of payload into orbit and providing up to 2 MW of power. This will enable remarkable data collection such as high-definition video and drilling into the planet’s surface. The technology provides a path to terrestrial fusion as well as eventual human missions across the entire solar system. The Phase I study will focus on creating higher fidelity models of the engine performance to enable optmization of possible mission trajectories and better quantification of the predicted specific power.

Continue reading

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

Elements for a 2000 satellite constellation? No problem!

A customer recently asked if we had any constellation design tools in the toolbox, which at that point we did not. This customer had been tasked with simulating a constellation of nearly 2000 CubeSats, and had available only a GUI-based tool that required him to enter each satellite’s elements individually. He was able to clone them, but still had to open up each one to edit the elements. This means he first had to make a table of the elements he wanted, then painstakingly enter them. In the Spacecraft Control Toolbox, creating these elements and simulating or visualizing the satellites is as easy as writing a for loop.

Our workflow using the SCT is mostly programmatic, using functions and scripts. Functions such as El2RV make it easy to go from Kepler elements to a Cartesian state for initializing a simulation. In just a few hours, I was able to write a function to generate elements for a Walker-Delta constellation of any size, and plot the results; to examples are shown below.

The WalkerConstellation.m function will be in our 2015.1 release, due out in a week. It can generate elements for a classic rosette or a polar star. If there is sufficient customer interest we may expand the constellation design tools available in the toolbox!

%   Form:
%   elements = WalkerConstellation( t, p, f, inc, sma, doStar )
WalkerConstellation( 720, 24, 2, 65*pi/180, 6800 )



% Create a polar star similar to Iridium
WalkerConstellation( 66, 6, f, 86.4*pi/180, 7150, true );



Full function header:
%   Compute orbital elements for a Walker constellation.
%   Generates a delta constellation be default, also called a rosette. For
%   a star geometry, pass in true for the optional doStar parameter. The
%   notation is i:n/p/f; f, the relative spacing, is an integer which can
%   take any value from 0 to (p-1).
%   Real-life examples include the Galileo, a delta geometry, and Iridium,
%   a star geometry.
%   Since version 2015.1
%   Form:
%   elements = WalkerConstellation( t, p, f, inc, sma, doStar )
%   -----
%   Input
%   -----
%   t         (1,1)    Total number of satellites
%   p         (1,1)    Number of orbital planes
%   f         (1,1)    Relative spacing between satellites in adjacent planes
%   inc       (1,1)    Inclination (rad)
%   sma       (1,1)    Semi-major axis (km)
%   raan      (1,1)    RAAN spread, optional. The default is 2*pi.
%   ------
%   Output
%   ------
%   elements   (6,t)    Kepler elements
%  Reference: Larson and Wertz, Space Mission Analysis and Design, second
%  edition (1996), "Constellation Patterns", p. 191

Simulating Magnetic Hysteresis Damping

CubeSats have caused a renewed interest in magnetic control of satellites, and passive hysteresis damping in particular. Modeling actual hysteresis rods on a satellite is not trivial, and generally requires empirical data on the properties of the rods selected. Our newest CubeSat simulation demonstrates damping using rods in LEO. A permanent magnet is modeled using a constant dipole moment, and we expect the satellite to align with the magnetic field and damp. We evaluate the results by plotting the angle between the dipole and the Earth’s magnetic field and the body rates.

First, let’s verify the magnetic hysteresis model in the toolbox using the bulk material properties in orbit. We use a dipole model of the Earth’s magnetic field. The nice hysteresis curves below confirms that we are computing the derivatives of the magnetic field correctly in the body frame, which requires careful accounting of rotating coordinates. Also we stay within the saturation limits which means our magnetic flux derivatives are correct too.

Hysteresis curves from simulating magnetic hysteresis in orbit

Hysteresis curves from simulating magnetic hysteresis in orbit

We will assume the rods are 1 mm radius and 95 mm length, with rods placed perpendicular to each other and the permanent magnet. Three rods are used per axis. The apparent rod parameters are taken from the literature. The actual rods will not reach saturation while in orbit, so we will see a minor loop.

Minor loops from damping rods

Minor loops from damping rods using apparent properties

The rods produce only a small amount of damping per orbit, so we have to run for many orbits or days to see significant damping – in some passive satellites, the total time allotted for stabilization is two months! In this case we test the rods’ ability to damp the torque induced by turning on a torque rod with a dipole of 1 AM2 and allowing the CubeSat to align itself with the magnetic field, starting from LVLH pointing.

Damping in LEO using hysteresis rods

Damping in LEO using hysteresis rods

Simulating the rods is time-intensive, with a timestep of about 4 seconds required – which makes a simulation of several days on orbit take several minutes of computation. Once performance of the rods has been verified, a simple damping factor can be substituted.

This new simulation along with the functions for hysteresis rod dynamics will be in the new version of our CubeSat Toolbox, due for release in June!


  1. F. Santoni and M. Zelli, “Passive magnetic attitude stabilization of the UNISAT-4 micro satellite”, Acta Astronautica,65 (2009) pp. 792-803
  2. J. Tellinen, “A Simple Scalar Model for Magnetic Hysteresis”, IEEE Transactions on Magnetics, Vol. 34, No. 4, July 1998
  3. T. Flatley and D. Henretty, “A Magnetic Hysteresis Model”, N95-27801 (NASA Technical Repoets Server), 1995

PSS featured on Time.com

Only two days after a visit by journalist Michael Lemonick, our DFD fusion drive was featured in his post on Time.com’s science section!

Going to Mars via Fusion Power

The article does misstate that Sam Cohen is a PSS engineer, when in fact he is the lead researcher on the PFRC at Princeton Plasma Physics Lab (PPPL). The proposed NASA mission to the Jupiter Icy Moons was JIMO – Jupiter Icy Moon Orbiter Mission. JUICE JUpiter ICy moon Explorer is a European Space Agency mission.

For more information on DFD, go to our fusion propulsion page.

Direct Fusion Drive