David Yaylali

Aerospace and Physics

Included here is a selection of neat computational tools and analyses, the main purpose of which is to show a broad selection of topics I have experience with. I'll continue to upload new material as I'm able. I've presented this material in a very pedagogical manner so that it can be understood by non-experts in each specific field.



Dynamics of a Launch to Orbit

To launch a spacecraft into orbit using chemical rockets, the primary goal is to raise the horizontal speed of the vehicle such that the vehicle, under the influence of gravity, falls around the Earth. This was nicely conceptualized by Newton using the cannonball thought experiment: Newton's Cannonball So, despite the fact that rockets are typically launched vertically, it is crucial that the rocket at some point "pitches over" so that it can provide thrust parallel to the Earth's surface, and pick up this necessary "horizontal" speed.

It so happens that this "pitching over" is naturally accomplished by the gravitational dynamics of spherical bodies: as long as there is some small initial pitch angle, the rocket will naturally increase its pitch as it accelerates away from the surface. This is known as a gravity turn maneuver. (Note that the dynamics here are completely intuitive: If a ball is thrown in the air in some non-vertical direction, the velocity vector of the ball will "pitch-over" as it follows the arc of its ballistic trajectory.) Typically, the initial pitch angle will be initiated by some onboard controller either early on, when the vehicle velocity is low, or later in the launch once the rocket is high enough to substantially reduce atmospheric drag; otherwise, the change in direction caused by the pitch-over will induce immense aerodynamics stresses on the vehicle, which could cause catastrophic failure.

The main point, however, is that to achieve a certain orbit (say, a circular orbit at 300km altitude), the launch vehicle must achieve a high-enough "horizontal" velocity (in this case, 7.726 km/s). This involves a certain amount of thrust, and also a certain pitch-over maneuver. As we now illustrate, the actual trajectory is extremely sensitive to these choices.

This simple MATLAB program serves to simulate the launch dynamics of a single-stage rocket which pitches over from vertical to some initial prescribed pitch-over angle, at some prescribed altitude. For simplicity, this code approximates the Earth as non-rotating so that the dynamics can be captured in two dimensions; that is, whatever direction we choose to pitch over determines the plane of the Keplerian conic section. I've given a more detailed description of the dynamics involved (including some details on the full 3D picture and how to incorporate the rotation into the simulation) in this set of notes. To illustrate the results of running this code, we consider a rocket with wet-to-dry mass ratio $n=15$, thrust of $9.33\times 10^5~{\rm N}$, specific impulse $I_{\rm SP}=390~{\rm s}$, and some further assumptions on aerodynamic parameters of the launch vehicle.

We first choose a specific pitch-over angle, $\gamma_{\rm po}=89.92^\circ$, and vary the altitude at which this pitch-over occurs. The altitude and speed of the resulting trajectories are shown below. Launch Trajectories It is clear that varying the pitch-over altitude by meters has a drastic influence on the trajectory. We see that a pitch-over altitude of $h_{\rm po}=\{130~{\rm m},135~{\rm m}\}$ gives trajectories which essentially escape the atmosphere ($h\gtrsim60~{\rm km}$) with velocities exceeding orbital velocity. Of course the flight-path angle associated with this velocity is important, since the velocity must be mainly horizontal in order to achieve orbit. It turns out (not shown) that velocity is indeed mainly horizontal for these two cases; we show the resulting orbit for $h_{\rm po} = 130~{\rm m}$ is illustrated below: Launch Trajectories (It turns out that this is not, in fact, a stable orbit, since at perigee the trajectory sinks substantially into the atmosphere and decays due to drag. So in this case an "orbit-raising" maneuver must be performed near apogee.)

We now illustrate what happens when we keep the pitch-over altitude constant ($h_{\rm po} = 125~{\rm m}$) and vary the initial pitch-over angle $\gamma_0$: Launch Trajectories Similarly, we see that even the slightest of differences in this initial angle (i.e., a quarter of a tenth of a degree) makes a drastic difference in the resulting trajectory. This sensitivity to the initial conditions suggests non-linearity in the dynamics. In fact, the dynamics are indeed non-linear, as outlined in the notes I have linked above. So, although these dynamics are not particularly complicated, optimizing the trajectory is a far more challenging problem. There are some analytical ways of doing this in special cases, but typically this is done numerically using various software tools. I hope to include a discussion of launch optimization soon, so stay tuned.

Back to Top

Generating Porkchop Plots!

When planning an interplanetary mission, say, from Earth to Mars, it's important to plan your route wisely. It's also important to plan your departure time. At any given time, Earth and Mars will be in different places along their elliptical orbits, and move along their orbits at different speeds. As a result, the needed fuel and time-of-flight for the trip vary drastically depending on the this departure time. This information is efficiently encoded in the porkchop plot: I guess they're called this because they can vaguely resemble a porkchop... To me it looks more like an ahi tuna plot:

The idea behind porkchop plots is simple: along the x-axis is the departure time (usually in units of days or months for interplanetary missions), and along the y-axis is arrival time at the destination planet. Upon these axes, one can plot a variety of things, including time-of-flight for the journey and the amount of fuel needed. The latter quantity is typically plotted indirectly; it is often convenient to separately plot the characteristic energy \({\rm C}_3\), which describes the needed energy to initiate the trip, and \(v_{\infty}^{-}\), which describes the arrival speed with respect to the target planet. This is convenient since the trip is typically initiated by large launch systems which in many cases expend their fuel completely (and are jettisoned) by the time the spacecraft leaves the influence of the Earth. The quantity \(v_{\infty}^{-}\) then characterizes the amount of on-board fuel needed (i.e., the internal fuel within the spacecraft) in order to "stop" at the destination planet.

Ideally, it would be nice to reduce all three quantities (TOF, \({\rm C}_3\), and \(v_{\infty}^{-}\)) simultaneously. Since fuel spent and time-of-flight are typically inversely proportional, however, there is some trade-off that is invariably necessary. Porkchop plots are key to understanding these necessary trade-offs in order to plan the mission.

Attached here is a simple MATLAB code used to generate porkchop plots for trajectories between any two planets within the Solar System. The date for the trip is a free parameter that can be set. For the default choices within the code (Earth-to-Mars, 2005), this code outputs... This compares well to the plot shown on the Porkchop Plot Wikipedia page. Note that to run this code, the MATLAB Aerospace Toolbox along with the JPL Ephemerides support package is needed. The code works by solving Lambert's problem at each grid-point in the porkchop plot. The specific method of solution to Lambert's problem used in this code is outlined in Section 5.3 of Orbital Mechanics for Engineering Students, by H. Curtis (2005).

Back to Top

Consensus Control for Spacecraft Formation Flying

In cooperative consensus control, vehicles share information with one another, and this shared information can be used to build control laws which bring the system as a whole to some consensus in position and velocity.

LISA Formation Animation

For instance, for the above system, SC1 receives information from both SC2 and SC3, while SC2 and SC3 only receive information from SC1. For this communication topology, a simple consensus control law takes the form \begin{align} u_1 &= a_{12}(x_2-x_1) + a_{13}(x_3-x_1) \\ u_2 &= a_{21}(x_1-x_2) \\ u_3 &= a_{31}(x_1-x_3) \end{align} where \(x_i\) is the state of spacecraft \(i\). The \(a_{ij}\) coefficients capture the communication flow in the system, and in this case only the non-zero \(a_{ij}\) factors are shown. Since SC2 and SC3 do not communicate, SC2 does not know the position/velocity of SC3 and vise versa... that's why there's no term in the control law proportional to \((x_2-x_3)\).

This control law therefore feeds back the position and velocity differences to the dynamics of the system (when this information is available), so that as the system evolves these position/velocity differences are driven to zero. That is, this control protocol drives the system to consensus.

Under control laws of the above type, the closed-loop dynamics can be written in a very concise way: $$\dot {\pmb x} = (I_N \otimes A - L \otimes BK)\pmb{x}.$$ Here, \(A\) captures the natural dynamics of each vehicle, \(B\) captures Newton's second law (that forces cause accelerations), \(K\) captures the relative weights of the position and velocity feedback terms (taken to be 1-to-1 in the \(u_i\) control laws written above), and \(\pmb{x}\) is a vector containing the collection of all vehicle states. The graph Laplacian matrix \(L\) captures the communication topology.

Consensus control laws can be used to build vehicle formations. As discussed in my masters thesis (see § 3.3), one way to do this is to define virtual leaders, and build a leader-follower consensus control law for each vehicle. If the virtual leaders follow some prescribed trajectory, the physical vehicles will then be drawn to these desired trajectories, regardless of initial conditions.

As a neat example, consider the Laser Interferometer Space Antenna (LISA), a proposed space-based gravitational wave observatory. LISA Formation Diagram The idea for this mission is to place three spacecraft in a triangular formation, the center of which orbits the Sun. The spacecraft are therefore in a circular relative orbit about each other... that is, if you remain at the center of the formation, you will see the three spacecraft move along a circle, separated by 120°.

Imagine how such a formation is initially built. One would assume that all three SC are launched by a single launch vehicle to the proper Earth-trailing heliocentric orbit: LISA Formation Animation Taking the position of the initially bundled spacecraft as the center of relative-orbit coordinate system, one can now utilize three virtual leader-follower consensus control laws to properly position each spacecraft along a circular relative orbit. The formation of the LISA constellation under this control protocol is illustrated below: LISA Formation Animation We can also plot the LVLH position vector magnitudes and interspacecraft distances as a function of time. The former should approach a constant (indicating circular relative orbits); the latter should approach a constant as well, indicating that the formation becomes a constant equilateral triangle. LISA Maneuver Geometry

Again, for a more detailed description of this maneuver, see § 3.3 of my masters thesis. The MATLAB code used to simulate this system (and generate the above animation) can be downloaded here.

Back to Top

Simulating Rigid Body Motion

The spinning "T-handle" in space!

The natural motion of a rotating body can be completely surprising to most people (including myself!) when first viewed in the proper setting. For instance, when the "falling" motion of a body is removed from the picture (as it is in orbit, on the ISS), a rotating body such as this "T-handle" behaves very strangely:

This motion of course arises from the coupling between the orthoginal angular velocity components in Euler's equations (with zero external torque), \begin{align} I_1 \dot{\omega}_1 + (I_3-I_2) \omega_2 \omega_3 &= 0 \\ I_2 \dot{\omega}_2 + (I_1-I_3) \omega_1 \omega_3 &= 0 \\ I_3 \dot{\omega}_3 + (I_2-I_1) \omega_1 \omega_2 &= 0, \end{align} where \(I_i\) are the principle moments of inertia of the rotating body and \(\omega_i\) are the components of angular velocity (in the rotating frame fixed to the body) with respect to the principle axes. For a very nice exposition of this material, written by my undergraduate research adviser Nick Wheeler, see here. Long story short, a rigid body such as a book or this T-handle can stably tumble about its "shortest" and "longest" axes, but not about its intermediate axis.

Simulating this rotating "T-handle" numerically is a very illuminating exercise, as it involves conceptualization of and transformation between rotating reference frames, dynamics of rigid bodies, and numerical integration. This cute little code performs this within MATLAB in the formalism of the Lie group \({\rm SO}(3)\), the rotation group. The state of the rigid body at any time is represented by the group members \(R\), which are simply the rotation matrices representing its current orientation. In terms of \({\rm SO}(3)\), Euler's equations (without torque) can be expressed as \[ J \dot{\boldsymbol{\omega}} = - \boldsymbol{\omega}^{\times}J \boldsymbol{\omega}. \] Using this code, one finds a very nice replication of the strange behavior seen in the above video (the black, red, and blue rods represent the (orthogonal) principle axis corresponding to the T-handle): animated T-handle

Back to Top

Decay-Chain Kinematics in Python

...or, building a Monte Carlo particle accelerator simulator

The Standard Model is our best description of the quantum mechanical nature of our universe. Along with the formalism of quantum field theory, the Standard Model (SM) allows us to predict the results of any particle physics experiment, and allows us to do so with incredible precision and (so far) accuracy. That said, we know the SM is incomplete. One particularly dramatic failing is that it does not account for dark matter, which has been indirectly determined (by considering several independent phenomena) to account for 80% of the matter content in our universe. I will not attempt to describe the beauty of the Standard Model or explain dark matter here (follow the above links for good introductions). And though I will use some standard particle physics terminology, I intend the following discussion to be understandable by anyone with a mathematical or engineering background.

Imagine that, in addition to the matter particles of the SM, there exists a new particle $\phi$, along with a whole family (or, tower) of new particles $\chi_i$, $i = 0,1,\ldots,n$. (To be a more technically precise, we assume the $\phi$ to be a scalar particle charged under SU(3) color, and the $\chi_i$ are assumed to be spin-1/2 fermions). We will assume that the indices of the $\chi_i$ particles are ordered by mass (so that $m_{\chi_0}< m_{\chi_1}< \cdots $), and also that $m_{\chi_n} < m_{\phi}$. We introduce a new fundamental interaction between these states and the Standard Model (SM) quarks, which is given by where at left we show the Feynman diagram and at right we show the interaction Lagrangian. The details of this Lagrangian are unimportant for our current purposes; it allows us to compute, in the formalism of quantum field theory (QFT), the probability of interaction between the different particles $q$, $\phi$, and $\chi_i$. This interaction therefore connects the "normal" matter we are familiar with to the exotic new particles $\phi$ and $\chi_i$.

What is the motivation for this "Beyond the Standard Model" (BSM) physics? Well, due to the fact that $m_{\chi_0} < m_{\phi}$, the $\chi_0$ particle is completely stable; conservation of energy, along with the interaction structure illustrated above, prevents the $\chi_0$ from decaying any further. Moreover, if the couplings $c_{iq}$ are small enough and/or the mass of the $\phi$ is large enough, then interactions between $\chi_{i}$ and SM particles (protons, neutrons, etc) will be extremely weak. The $\chi_0$ particle can therefore be identified as dark matter!

How would we look for experimental signatures of this model in the lab? Since $\phi$ carries color charge, a pair of $\phi$s can be created by colliding two protons together at a high enough energy, such that $E> (2m_{\phi}) c^2$. We can therefore attempt to look for signatures of this model at particle accelerators like the Large Hadron Collider. A particular signal which would arise is illustrated below: In words, colliding protons can produce a particle-antiparticle pair of $\phi$ particles. Each of these $\phi$ particles will then decay, in steps, down to the lighter states ($\bar{\chi}_a$ and $\chi_b$ in the figure), and in doing so will emit a number of quarks (which are visible in colliders, forming what are referred to as "jets"). Notice that the entire decay chain is constructed by "stitching together" copies of the 3-point fundamental interaction given in the first figure.

Our goal is to simulate events given by the above topology; in essence, this will allow the detector teams at the LHC to search for this "new physics" signal. Of course, since the individual interactions making up the decay chain obey quantum mechanics, the exact decay path taken by each $\phi/\phi^{\dagger}$ particle is not deterministic. For instance, if there are 3 particles ($\chi_0$, $\chi_1$, and $\chi_2$), there is some probability of the decay chain $\phi \rightarrow \chi_2 \rightarrow \chi_1 \rightarrow \chi_0$, some probability for $\phi \rightarrow \chi_2 \rightarrow \chi_0$, and so on. We therefore need to look at many events, and use the jet number or jet momenta distribution to search for evidence of our model. Due to this quantum mechanical (probabilistic) behavior, we require a Monte Carlo simulation.

In other words, we require a Monte Carlo "particle accelerator simulator". Many such simulators exist (such as MadGraph) and are commonly used in HEP research. These tools, however, are ill-suited for the process illustrated above, since the combinatorics of such long decay chains quickly become far too computationally expensive. To deal with these types of events, where large "swaths" of jets are produced from a large number of decays, I created PySwath.

PySwath can be downloaded here. This Python program takes as input the model parameters and initial momenta of the $\phi \phi^{\dagger}$ pairs, then generates the decay chain based on the analytical probability distributions which can be computed from the interaction Lagrangian $\mathcal{L}_{\rm int}$. The program outputs a file which contains the identity and four-momenta of each particle in the decay chain, for all events generated during the simulation run. For illustrative purposes, this program also outputs histrograms (as .pdf files) showing the jet number distribution and the transverse momenta ($p_T$) distribution of the highest-$p_T$ jet: PySwath Graphs Note that only an even number of jets are created. The reason for this can be directly seen from the decay diagram shown above. First, both the $\phi$ and $\phi^{\dagger}$ will decay to some $\bar{\chi}_i$ and $\chi_j$ state, with each decay contributing a single quark-jet (a total of two jets). Any additional decay process within either leg will then give an additional two jets, e.g., $\chi_j \rightarrow q \bar{q} \chi_k$, with $k < j$. The total number of jets from the complete process must therefore be even.

Some additional notes:

  • This code is written in Python 2.7. It is packaged such that, when extracted, one can simply do:
        $ python PySwath.py
    in terminal. With default choices, there are 90 $\chi_i$ states ($\chi_0, \chi_1, \ldots, \chi_{89}$) and 1000 events are generated. Event generation takes about 15 seconds on my laptop.
  • A far more detailed description of the model illustrated above, along with a parameter-space analysis based on cosmological constraints and collider searches, will be published shortly (and I will link to the paper here when this paper is completed).
  • The initial production of the $\phi \phi^{\dagger}$ pair must be computed in MadGraph. I have included production data for a specific case (where $m_{\phi}=1$ TeV) so that the code can be run out-of-the-box. I've decided not to include production in this code since MadGraph has no trouble computing a simple 2-to-2 scattering process such as $p+p \rightarrow \phi + \phi^{\dagger}$.
  • This code does not take into account certain quantum mechanical effects (such as quantum interference between different processes). For our purposes these effects are essentially negligible.
  • The output file storing the particle identities and four-momenta is formatted according to the Les Houches Accords, so that it can be passed to other industry-standard programs for hadronization and detector simulation.
  • This code is actually a skeleton of the actual code we are using for our analysis, and just for illustrative purposes. The full code is structured slightly differently and contains a more extensive event analysis.

Back to Top