Comet Dynamics and Implications For Orbiting Spacecraft

So it's been a while since my last update. I've completed my masters program, gotten a job, moved to Tucson Az, and haven't had much time to MatLab. Also, my MatLab subscription ran out! OH NO! I use it at work but I can't use it for my own uses. So I have been learning c++ and will continue to update this blog with whatever I create using that.

So I thought I would add my thesis to this blog. Below is a link to download the PDF file. Let me know what you think! (Honestly, I thought I skimped out a little bit on the orbiting satellite analysis. If I had another month, it would have been epic.)

Comet Dynamics and Implications For Orbiting Spacecraft

Quick Update

Thesis is done and has been defended. Once it has been published I will include a link here. Now I just need to finish the semester and find a job...

A few things have been left out of the thesis analysis since I ran out of time, such at the stability of the comet's rotation and conduction between discretized faces. In time I will include them into the comet simulation and ultimately into SpaceSim, which I plan on restructuring to allow for more modularity and less glitchiness.

I haven't written anything about the satellites interaction with the outgassing jets. Below is a figure of the orbiting satellite, but the body of the text is in the thesis, I don't feel like writing it all out right now. There is also a video that goes along with it, but it isn't optimized since it wasn't needed in the thesis (damn paper not displaying moving images) and I just needed a snapshot (as seen below).

Orbit of a satellite with/without interaction with outgassing jets. 

Other then that, I helped some colleagues develop an app to diagnose anemia by looking at the color of eyes/nails. I did the coding and created a GUI for it, but I can't discuss more then that since there is some no-talking-about-it thing, but it's pretty cool and works well. 

Thesis and things

For the last while I have been hard at work on my thesis, which mainly encompasses what what I have written about. Right now I am sitting at a solid 29 pages (written, 40 total including all the random bull-stuff), with 3 of the 5 chapters are done, with chapter 1 being the intro, 2 is the comet model with sublimating jets and whatnot, and 3 is the dynamic expressions used. The 4th chapter will be the simulated results and the 5th is analyzing a satellite orbiting the comet.

So, I will only talk about the new stuff. 

Going off the last post, I added in a variable shape to go along with the change in inertia and mass. Essentially I just made the face vertices into a variable. It's a little buggy right though. The way it is written, the faces translate, but do not change area. Therefore after many orbits, it becomes very inaccurate. I have been playing around with ways to fix that, like changing the area as a ratio of volume to surface lost, however there are implications for the other faces which use that vertex. So that is going on the back burners for now.

 I also have written a way to describe the outgassing jets, which can be seen in the figures below. 

Left: Density/Pressure ratio throughout the outgassing jet 

Right: Geometry of the outgassing jet

To "get around" the fact that the faces are triangles, a circle of equal area is placed at the center, and is considered the pseudo-face. Beta functions (alpha=beta=2) are then used to describe the distribution through the jet. While the density at the face reaches 150%, it does not make a difference in the far field. It can be remedied with a beta function varying with distance from the face (starting at alpha=beta=1 which is a square wave), but I decided it was more trouble than it is worth. Using these definitions, the effect of outgassing jets on a orbiting satellite can be found with a different simulation I have also been working on, but that results from that are for a different time.

Comet Dynamics II

So I finished the comet model for now (for a second time). I have reduced how long the model takes to run again from 33 seconds to 14 seconds by taking out the cross function and replacing it with the direct equation.

Changes from the last post:

• Variable inertia matrix

The inertia matrix now changes with the ejected mass. This makes the simulation more realistic. However, the shape of the comet does not change (ie. the change in mass does not change the comet thickness). 

• Random off-gassing jet placement and variable composition

Instead of off-gassing jets being able to be produced on every face, now the chance that a jet will be on a specific face is randomized (I currently have it at a 5% chance). The composition of those faces are then randomized at a maximum of 5% of the face being volatilizes. This also effects the initial inertia matrix.

One implication of this is with constant face jets, the comet usually will spin up along a specific rotation axis. This can be seen in the following video. The great thing is that the axis it settles on does not have the be a major or minor axis due to the off-gassing thrust.

Black line: Rotation axis.

Yellow line: Sun vector .

Comet Color: Proportional to temperature (Red is hot).

Dark blue line: Total force on body.

Rings: Torque direction (green is +, red is -) and normalized magnitude. 

Light blue lines: Direction and normalized magnitudes of off-gassing jets.

Off-Gassing Dynamics

The off-gassing of an asteroid also has implications in satellite dynamics. However, different size and velocity particles are created and can effect a spacecraft in different ways, ie. it can totally destroy it, puncture a solar panel, provide a torque, etc.. All of which would make precise data acquisition difficult. The second part of my thesis project will be analyzing this.

• Off-Gassing

On top of the trust off-gassing provides to the asteroid, it also provides a pressure on spacecraft, as well as accelerates dust particulate off the surface of the comet. To be able to analyze this, I have made a simple model to find the velocity that specific size dust would travel at. The velocity of each particle is a function of the radius of the comet, the coefficient of drag, the solar flux density, a diffusion coefficient, radius of the comet, the gas ejecta density and velocity, and the particle density and radius. 

In most literature I have seen, out-gas jets are taken to have a half-angle and not diverge from that. However, since the pressure will not be constant through the coma, I created the diffusion coefficient. There is no empirical data behind this except that a high pressure jet would not stay at a constant half-angle due to pressure differentials between the jet and a low density area / vacuum, ie. it would defuse. Also, literature says the majority of mass will be at 20% of the velocity of the jets and usually be between 20 and 100 m/s. My off-gassing model adheres to that for dust larger than 1 mm. I believe this makes sense, but that's because I made it this way. 

The following plot is the results from this model:

Dust Radius: [0.0001;0.0005;0.001;0.0025;0.005;0.01;0.02;0.03;0.04;0.05] (m)

Gas Velocity: 400 (m/s)     Gas Density: 0.004 (kg/m^3)     Dust Density: 300 (kg/m^3)     

Diffusion Ratio: 1.5           Coefficient of drag: 0.5              Comet Radius: 500 (m)

Left: Velocity of 10 different particles

Middle: Distance from initial location of particles

Right: Median acceleration of the particles

Conclusions:

While I showed one run of the model here, it can be seen that the particles settle on a certain velocity relatively quick (withing 15 seconds for any "reasonable" parameter combination). I need to do further research or develop a formula to find the velocity from the above parameters without having to run a quick simulation. On top of that, at close quarters to the comet (within 1000 meters, the acceleration due to the gas dominates, whereas the acceleration due to gravity and solar radiation pressure gain prominence after a while. 

So whats next:

• Satellite Dynamics

Here is a quick model of the satellite I will be using (it's the same as past satellite pictures, but now it has an "engine" in the back and a "camera" in the front). Behold the JimSat in all it's glory:

JimSat

I will be using this to analyze a satellites ability to acquire data in highly dynamic offgassing/asteroid environment.

• Coma and Tail Model

I'm not too sure if I will be able to do this, but with the given trajectory and composition I want to predict the tail and brightness of the coma. Maybe even predict if the comet will give the Earth a meteor shower, I don't know.

Comet Dynamics

So I finished the comet model for now. It works well for me to move on and work on other parts of the model. I have verified it with literature and it gets similar results (which was a nice thing to see), but it still needs some work. It should be noted that the only "volatile" taken into account so far is ice/water vapor. 

Before I get into the bulk of this post, I found a few edits which make the model run magnitudes faster (~33 seconds opposed to 40 minutes). 

   - Due to the differential equation solver (ode113), if the axes lengths of the ellipse are equal or if the rotation axis is along a ellipse axis, the model will take forever (about 40 minutes, depending on simulation time) to run. This is because the differential equation solver wants to keep the angular velocity along the same axis, or having an angular acceleration of [0;0;0]. By making the rotation axis slightly offset from an ellipse axis, this is avoided since it neither wants to keep a constant rotation axis as well as constant angular velocity. I currently have the rotation axis at [pi/1000;pi/1000], or 0.18 degrees, and it works amazingly.

• Comet Model (Asteroid model with off-gassing)

The difference between the asteroid model and comet model is the thrust and energy use from off-gassing. The temperature of the face is used to predict the sublimation rate. Since the comet is in a vacuum, the velocity of the gas should be mach 1, which is found easily using values for water vapor for molar mass and heat capacity ratio. The force is also found and is shown to greatly effect the rotation of the comet.

So without further explanation, here are the same plots and videos from the last two posts, but with off-gassing. The first plot is a comet with axes a=5000, b=3000, and c=2000;, the second is scaled down by 100. 

Large Comet:

Small Comet: 

Conclusions: 

The time scale for both the first and second simulations was the same (just under 1 day's time). Notice the relatively constant temperature of the fast rotating comet, as shown in the video, this makes the force tend to be in the direction of travel, instead of along the anti-sun vector. This has some large implications for small comets because this means that they would (rather quickly) increase their semi-major axis until the off-gassing reduced to a negligible amount. Also, it can be concluded that small comets either will break apart due to their spin rate, become tumblers, or break apart from internal pressure build up. 

YORP Update

This will be quick post. I have found a few errors in the YORP model and corrected them. I was wondering why the force plot looked as it did, the reason is because it was only recording one face. As you can see below, this has been fixed.

Following are a few things to note and the new output plot as well as a video I made to validate some aspects of the model.

Things to note:

• The torques and forces are in the body frame. To make it useful it needs to be put into the orbital frame (check the video).
• The average temperature reaches steady state as the number of faces increases and the asteroid becomes circular about the rotation axis (if the rotational and orbital angular momentum axes also align).
• The warmest time on the asteroids is in the afternoon (around 2-4).

For the video below, the yellow line is the direction of solar rays, the red line is the direction of the force due to solar radiation pressure and thermal emittance. The color of the asteroid is relative to it's temperature with the colormap [1 0 0] (red) being the total hottest temperature and [0 0 0] being absolute zero. 

Asteroid Rotation Video viewed along the Z-axis

(Sorry about the resolution, uploading it requires compressing)

New output plot

Asteroid Dynamics

Finals are done and I am on winter break, so I finally have time to work on my thesis model.

Before getting into this, here are a few definitions:
    -Solar Radiation Pressure: Momentum transfer from photons which provide a force onto any body.
    -Yarkovsky Effect: This effect takes place due to the radiation pressure that a rotating body emits. Due to it's rotation, the force applied to the body can, over large periods of time, alter the orbital parameters of the asteroid.
    -YORP Effect: Second order Yarkovsky effect where, if the force is off axis, it provides a torque to the body and affects the rotational state of the asteroid.

Over the last few days, I have been putting a YORP model together, which I will then put into a Yarkovsky / bulk model, which I have yet to create (I had to make one first, I decided I would make the YORP). On top of that, I also need to add a Cometary function into the YORP model, which will be done after the bulk model. This model will include off-gassing. The total model will be broken up into bulk and body parameters. The bulk parameters being orbital based (radius and velocity) while the body parameters being rotation based (angle and angular velocity).

The way I currently see it going, the YORP model will be nested inside the Yarkovsky. I hope to find a way to simplify the model in an individual time dependent equation for force and torque (change in body and bulk parameters), but I don't think this will be possible for a more accurate asteroid model (not an ellipsoid). On that note, right now I use an ellipsoid for ease of computation. The lower the number of faces, the quicker it works. Once I have the model validated for the simplified model, I will run it on a more complex shape model.

Above is current output. Going counter-clockwise from the top left, the individual subplots are as follows:

-Picture of the asteroid and the rotation axis (red)
-Average temperature
-Temperature of every face
-Unit quaternion describing the angle
-Angular velocity in x-y-z directions
-Total force on surface (x,y,z)
-Torque on body

Asteroid Gravity Models

I did this a while back but decided only now to post it. Right now I have 3 gravity models (2 working). The problem with gravity models is that for concave structures, a tetrahedron from a face to the origin may have elements external to the body or overlap other tetrahedra. This "creates" mass/volume outside of the actual object volume which can cause huge errors for simulating an orbiting satellite.

A few ways around this are:
   1). to create two "origins" and set the origin for each face as the one where the inner product (dot product) of the CCW face normal and radius vector are positive (if both are positive take the one with the smaller radius).
   2). Define the mass to exclusively be on the surface.
   3). Do some math and define the gravity as a field affected by the body volume.

Originally, I was only using a point mass tetrahedra model. This model took the faces, connected them to the origin (creating a tetrahedron), found the volume, and put a point mass at the center of the tetrahedron proportional to the volume and density. This model does not work for very irregularly shaped objects, however it works best for the majority of asteroids.

The other working model is a shell model. It takes all of the mass at the surface of the asteroid where each unit of surface area contains the same amount of mass. While this model is better for more irregularly shaped objects, by putting all the mass at the surface it creates a gravity singularity at the surface and can mess up data if trying to analyze an orbit close to the asteroid. This model is best for analyzing the gravity field for very irregularly shaped asteroids at a reasonable distance from the surface.

The following plots are for validating these two models. I will discuss the differences in the plots and what denotes a good "fit" at the end.

It should be noted that the scale at the bottom is the number of points in the analysis direction and does not correlate to the size of the objects.

Eros
    

        

Left: Point Gravity Tetrahedra Model

Right: Surface Area Model

Churyumov-Gerasimenko (67P)

Left: Point Gravity Tetrahedra Model

Right: Surface Area Model

What to take from this:

Point Mass Tetrahedra

• Underestimates radius since it takes the center of the tetrahedra.
• This conserves mass while avoiding the gravity singularity at the surface
• For very irregularly shaped objects (67P), the volume is disconfigured and does not give accurate information about the shape/gravity of the object.

Surface Area

• Preserves shape of irregularly shaped objects (67P)

Conclusions

• The Point Mass Tetrahedra model should be used for regular shaped objects.
• The Surface Area model should be used for irregularly shaped objects. But may create large errors is simulating a close approach.
• Another model needs to be created to better represent this.
• A mixture of the two may work. I am also looking into Spherical Harmonics as well as a triangulation method I read about a while ago. But for now these will do.

Not SpaceSim

I will discuss 2 different programs in this post: a 2D finite element analysis program and a Robot arm. In my next post I will discuss my progress for my thesis.

I will not go too in-depth, so if there are any questions feel free to comment. 

Here is a quick picture of Dawn. I found the .stl file online and wrote a quick parser to make my program able to plot it.

2D FEA:

Without getting to deep into it, I made an FEA program for my Intro to Finite Element's class (my partner is doing the 3D analysis in Abaqus). The object being analyzed is the Philae Lander's leg. We are doing static analyses for the stress, strain, internal forces, and deformation for a resting Philae as well as when it was compacted before the first bounce. Assuming no energy dissipation (pretty bad assumption since there was lots of kinetic dissipation), the maximum internal energy is proportional to the square of the velocity after the bounce, which is estimated to be around 38 cm/s.

Here is what it looks like:

2D FEA Program

Robotics:

While we have not yet completed the project, my main contribution so far has been this robot which is able to get out of a randomly generated enclosure. It does this by always turning right, after the first left. It saves where it has been as an artificial "memory" (very basic AI), and if it "sees" the same spot over and over it randomizes it's movements until a new wall has been found.

I created the simulated robot since we were having trouble interfacting with the real robot (PhantomOmni).

Oh yeah, on that note I figured out how to take complex videos on matlab, the buffer needs to be set as "zbuffer" otherwise it does not do something that it has to do.

Anyways, here is a video. For best results try to play both videos at the same time:

Robot Arm

Robot Artificial Memory

Thesis

I'm back with a quick update!

I have been sidetracked from the GUI for the time being with my thesis project, which will eventually turn into the Asteroid/Comet Toolbox along with the other asteroid functions I have made.

As of now, nothing has been coded yet. However, the method is set and the math is almost all mathed out.

Comet 67P Churyumov-Gerasimenko

There will be functions for:

• A new gravity model 
• The point mass model does not work for convex structures such as 67P seen above)
• Radiation pressure and emissivity
• Model of comet off-gassing
• Composition matrix
• Volatilizes sublimate and have set constants/parameters
• Change in orbital elements and rotation due to SRP and off-gassing
• In attempt to analyze the yarkovsky effect
• Trajectory of off-gassed elements
• If they will impact an orbiting object
• Ability to predict a tail
• Would need to include ionization model

Lucky for me, ESA currently has the Rosetta spacecraft in orbit about the comet 67P Churyumov-Gerasimenko (seen above) and the oort-cloud comet, Sliding Spring, is making a close approach to Mars in the upcoming week. I plan on using data from these missions to validate my model.

Rockets

I decided to finish the rocket editing before going on to create the tool boxes. It is similar to the previous post on rockets (found here).

I have created a video (a friend pointed me in the right direction for free video capture software) to show how the rockets are created and edited. 

How it works: 

     The user names the rocket and indicates the number of stages and parallel thrusters (0 assumes no parallel stage). Then, a preset rocket program is used to give an editable initial rocket. All physical characteristics are able to be edited, including but not limited to: masses, dimensions, specific impulse, and mass flow rate. 

By clicking update, the structure containing the object is updated and then plotted again with any changes. The masses and total deltaV are also refreshed. 

Next on my list of things to tackle is the impulsive trajectory solver.

And the world keeps spinning

So I was trying to make a stand alone movie capture and snapshot program before I posted this. I was successful in making the snapshot, however the movie poses problems when two while loops had to run simultaneously and MatLab is single thread (I do not have the parallel toolbox). It has a fix but it would take time that I am not willing to currently put into it. Maybe one day once the tools are done or I get a strike of inspiration I will do it. Regardless, here is a quick snapshot-ed image of the Earth. If you want you are able to change the rectangle of screen capture.

Round

The reason I wanted the movie capture program was to make it so I would not have to capture the following section on my phone (I don't know of any good free screen recording software). But anyway, here it goes.

As you can see from the new play, pause, rewind buttons in the above image, it has the the ability to (of course) play, pause, and rewind any plot pictured in the axis with a time step set by the user. The planets spin with their actual angular velocity and are also able to be viewed with satellites orbiting about them as well as orbiting the sun.

Earth spinning with multiple satellites orbiting. There is a 10 minute (600 second) time step between steps.

Mercury, Venus, Earth, Mars, and the asteroid Itokowa orbiting the sun with a ~2.8 day step.

Unfortunately, the buttons work in a while loop so I cannot make a video of them without including the capture code in the same loop.


Time to start working on the tools. I think I may have to write them as "plugins" to avoid the need to edit the structure of some of the other functions.

SpaceSim (GUI) 1.0

So here is the framework of the GUI. If there is any feedback, I would love to hear it. Remember, I am an engineer, not a computer scientist. For now, I have not compiled it to an .exe file and MatLab must be installed to use it.

So starting from the get-go: the main file is run from anywhere, it get's its current location and makes that, and all internal folders, part of the path. From there, the GUI is opened and populated. I will explain what each part of the figure does, then show examples of the capabilities.

Initial start-up figure of SpaceSim.

Starting from the top: 

File menu: Allows the user to save, load, or exit the GUI. The save option saves the current structure of objects, as well as any populated tools, which can then be loaded for future use/analysis. 

Edit Menu: Right now this is only a placeholder, I do not have any editable program properties yet (except for the objects, but I will get to that).

Add Object: Menu to add an object to the object list. The current addable objects are planet, rocket, asteroid, and satellite. I will speak about each of these objects in further detail in a later section.

Tools: While I have not implemented any of these into the GUI (YET!), they will consist of (but are not limited to):  an orbit simulator / impulsive and finite transfer optimizer for multiple or single central bodies (including asteroids), design and launch of a rocket, a Walker Constellation tool along with analysis abilities such as access times with constraints, and a mission design tool which splices together the other tools to create a rocket to launch a satellite to a defined orbit where further analysis can be done.

Help: A user manual which I am writing as I go along.

Toolbar Menu: I have included the generic MatLab figure toolbar since it contains all the functions I need. Until a time when I can create my own toolbar or map zooming/rotating to keys I will leave it as this, since I have more important things to do.

Time: The time is taken at the start of opening the program. This time can be edited, as long as it stays in a form which MatLab can read. To view the bodies at a different time, the plots must be updated manually.

View Dropdown Menu: Lists all the current bodies as well as an option for "Orbit", which allows the user to choose from the Earth and Sun (will eventually be every planetary body object) and takes every object whose "parent" is one of those two and allows the user to check boxes and view their orbits about the main body. 

Current bodies moving about the central body, the Sun. 

A picture of the orbits of these planets at the date indicated on the top right assuming the initial orbits are from January 1, 2000 and extrapolated to their current position via. Kepler's equations (which is when the planet data is from, but this will change in the future). The planets are currently not listed, but there will be text right next to the dot saying which planet it is. Planet Awesome = Yellow, Earth = Green, Jupiter = Blue.

Axes: Displays whatever is prompted to by the view dropdown menu. These are 3D and can be controlled via the figure toolbar.

Object List: A list of objects which the user has input. All of these objects are able to be edited and deleted at will. Currently all input angles are in radians because I am lazy. It isn't that hard to convert from degrees. Deal with it.

Current object capabilities. The editor has the same options as the create panel.

     Planets: To create a planet objects, there are two options. The first is to choose a defined planet from a list which contains all 9 (including Pluto) planets, the Sun, and Earth's Moon. All of these have images included which can looked at by changing the "View" dropdown menu to be the name of that object. The object will appear in the list as soon as it is added to the object structure. Alternatively, one can create their own planet and define all of the variables themselves, including the name, parent body, orbital elements, and any other parameter. Also (this is my favorite part), the user is able to add a .jpg image to be mapped over the planet to visualize it.  

Left: Add planet pop up.

Right: as well as the custom planet options. There will be more options in the future such as atmospheric density.

Left: Jupiter, as viewed from the program.

Right: Planet Awesome, including a picture mapped over it for the user's pleasure (, a picture which my friend wouldn't allow on facebook so I decided to defy him elsewhere).

     Asteroids: Same thing as planets only there is no creating an asteroid, since it must contain defined point and connectivity list. The parameters are able to be edited though. There are currently 4 asteroids which are programmed in, however I have a database of a fuckton (including their connectivity lists and all parameters) of asteroids which I have not converted to be able to be read by the program. Here is the image of Itokowa which I included in the last post. Right now the only thing the asteroids can do is be pretty. However (and once again as showed in the last post), I need to put the code to orbit asteroids into the GUI.

Itokowa.

     Satellites: The editable parameters of a satellite are similar to that of the previous two. The SA stands for solar array and the rest is self explanatory. Satellites are also able to be viewed, however there is only one model (rectangular everything).

Satellite edit options

Possible Satellite Configurations.

Rockets: Read this post for a full rundown of the rocket capabilities. While I have updated the code (uses inertia matrices and thrust vectoring instead of assuming a constant angular acceleration emulating a thrust vector), and implemented some of it into the GUI, this will give the gist of everything I have done with it and will put in. In short, you can create, view, and launch rockets. So far I have no put the ability to edit in rockets as I have with the other objects since it requires more effort to write the GUI and I would rather work on other parts because I already know it works. Following is a picture of the Falcon 9 Heavy, with projected staging data, as seen in the program.

Falcon 9 Heavy

I think that might be it for now ... I am pretty sure I missed something but I can always cover it some other time. 

Update

Hello!

I have not been posting about it for nearly a year due to school and nobody really reading this. While I have not worked on the Solar System Sandbox program for some time, I have still been working on a continuation of this project. My current project will be explained in the next post.


But before that, a brief update on what I have been working on:

Last semester I created a program which I dubbed "Asteroids", which does just that. Visualization and simulation of orbits about asteroids.

Visualization of Itokowa, which was visited by the Hayabusa spacecraft. 

Itokowa is one of the many asteroids already programmed in.

I included some functions from the Solar System Sandbox for orbit transfers and whatnot and used this program in a paper I co-authored titled "The Near Earth Object (NEO) Scout Spacecraft: A low-cost approach to in-situ characterization of the NEO population" (which can be found here: Neo-Scout Report). I worked specifically on the trajectory, proximity maneuvers, and attitude control sections, while editing it with a few others. Examples of these functions can be seen below (taken straight from the paper so my lazy ass doesn't have to re-run the simulations).

In the paper, the program is used to analyze the trajectory of a satellite orbiting about an asteroid with given face/vertices and density. It includes functions for station keeping and landing/taking off from the asteroid.

Other than that, I have been working realistic rocket program utilizing inertia matrices and thrust vectoring instead of assuming a set angle of turn.

And my current big project, I am making a GUI to put all of these functions together with improved graphics and usability. The frame for editing and visualization is currently finished and I will begin to add tools in very soon. But I will show all of this at a later time.

Orbit Transfer (Part 1)

Julian Date: 2456513

Last time: Plotted current and target orbit specified by the user and indicated initial/target locations of orbit

This time:

I finished the function for an orbit transfer. While I have only included the functions for an orbit transfer, the great part about the functions, namely the Lambert's problem function, is that it can be utilized for all of the transfer functions (rendezvous also require Wiltshire-Clohessey and deep space/moon requires compensation for the gravitational pull of the celestial bodies). This being said, the accuracy of this function decreases with the distance from Earth due to the increased gravitational pull of the other celestial bodies.

So on to the show:

So from the last post, the program has either launched a rocket or the orbital elements of the current orbit have been input. For this example, the elements for both orbits are as follows.

Left: Orbital elements of the two orbits. The top is printed from past inputs and the lower was just input in this frame.

Right: Current and target orbits and departure/arrival locations.

After this, there is an option to either choose the transfer time, or optimize the time for the smallest change in velocity possible. First we will go over the first option.

Choose the time:

After asking if the user wants to launch at a different anomaly, the program prompts for the time of transfer (for this example I intend to create a hyperbola, so I input 100000 seconds, or 27 hours 46 minutes and 48 seconds).

The Lambert's problem function is run and the trajectory is then plotted over the previous plot.

 Plot of the orbits and target trajectory. 

The program then runs through a bunch of equations to get the outputs shown below.

 As you can see, the transfer orbit is a hyperbola. The negative efficiency of the second turn means that the change in velocity was greater than the final velocity itself. I will fix this problem in the future. The rest is pretty self explanatory.

Optimization of the total Delta-V:

After choosing to optimize for a delta-V, the program iterates through Lambert's problem (starting at 10000 seconds, or 2 hours 46 minutes and 28 seconds, and taking 500 second steps) until the delta-V stops decreasing. I will confirm this is the most efficient path later with a function I wrote specifically for that reason. 

The trajectory is chosen from the absolute min and then plotted.

Left: Orthographic view of the trajectory and orbits.

Right: View in plane with original orbit and slightly out of plane with the transfer orbit.

It then plots the same relevant information as before.

Notice the elliptic eccentricity. All efficient transfer orbits should be an ellipse. 

So that is it for this function. I need to find an already developed program to run my results against to see how they stand up.

Proof of the delta-V efficiency:

I edited the code slightly to iterate for a certain amount of time and then plotted the delta-V and eccentricity against time, and the eccentricity vs. the delta-V. Originally I decided to optimize based on eccentricity. I made this decision based on the intuition that a hyperbolic orbit will vertically asymptote at t_transfer=0 and a an orbit would horizontally asymptote as t_transfer approaches infinity. However, and as you will see, the optimum delta-V occurs slightly after the minimum transfer orbit's eccentricity**.

Here is the data, I will zoom in and talk about each in detail after this figure.

Top Left: Delta-V vs. Time. 

Top Right: Eccentricity vs. Time

Bottom: Eccentricity vs. Delta-V

Top Left*: The absolute minimum on the delta-V  vs. time plot is at about 3.08e5 seconds.

Top Right*: The absolute minimum on the eccentricity vs. time plot is at about 2.72e5 seconds.

Bottom**: Shows the minimum delta-V, and therefore the eccentricity at which the efficiency switches direction. This will be discussed more thoroughly in the next few figures.

*As you can see, the absolute minimum for the delta-V occurs slightly after the minimum for the eccentricity. The eccentricity at 3.08e5 seconds is 0.8148, a change of in eccentricity of 0.0042.

As seen above, the minimum delta-V is around 2.051 at an eccentricity of around 0.081495. Notice again how the minimum delta-V is not at the minimum eccentricity (but it is very close).

**The plot on the right should remind you of a hyperbola (v-asymptote at x=0 and h-asymptote at y=0) and a logarithmic function (v-asymptote at x=0 and h-asymptote at y=1). This is because, as I alluded to it earlier, at lower transfer times, the trajectory is hyperbolic. Therefore the eccentricity approaches infinity as the time approaches zero (+). Also, as to not escape the gravitational pull of the earth, the eccentricity cannot exceed 1, therefore as time increases to infinity, the eccentricity would approach 1. The point where the trajectory is more efficient  should be at the minimum points of the plot because a higher eccentricity indicates a higher energy state. The lower the energy state, the less energy is needed to be added to the system. The reason the delta-V does not have a minimum has something to do with the change in inclination, and I will have to look more into that at a later time.

Well that was a lot. I will be working on the rendezvous function next, which shouldn't take long, followed by transfers to celestial bodies. I also plan on figuring out some kind of transcendental equation to optimize the orbit transfer, with variable arrival/departure anomalies, and time. But it is late and I am tired, so I will end this post riiigggghhhttt . . . here.

Rocket Launch (Continued: Part 3)

Julian Date: 2456509

Last time: Plots the trajectory of a rocket after the engines stop firing using an N-body approach.

This time:

Now, from the starting menu, the user has an option to do an orbital transfer. If this option is chosen, it asks if the user wants to either launch a rocket or input orbital elements to calculate the original location and velocity of the rocket. I will go through each in turn.


Orbital Elements:

As you can see, the program prompts for the 6 orbital elements (angular momentum found from eccentricity and semi-major axis). 

 It then plots the orbit and draws a line to the current location.

Rocket:

The program gives the same options to the user for launching a rocket as the rocket function. Once it is "launched" it takes the radius and velocity vector of the rocket and finds the 6 orbital elements from that. It then plots the projected orbit of the rocket, as seen below.

     The launch: Falcon 9 (light), 5000 kg payload, thrust vectoring at 20 seconds for 20 seconds at 30                degrees to east for 80 degrees.

Left: Projected orbit of the rocket in accordance with the orbital elements.

Right: Close up on the perigee end. Notice the line pointing to where the rocket currently is. Also, the rocket's trajectory does land it inside of the Earth (meaning it crashes) at some points on the perigee end. 

This figure shows the N-body version of an identical launch. Notice the similarities in  launch location and trajectory. If you look closely at the N-body projection, you can see the lines do not line up perfectly as they do with the orbital element plot. This is because of including the masses of all other celestial bodies instead of assuming a two body system and only using Earth's mass. While this model is more accurate in the long run, the orbital element model can be used to predict short term trajectories, and can be back checked with the N-body model and the trajectory can be adjusted if it gets off course.

Here the option is given for which type of orbital transfers the user wants. The options currently are a general orbit transfer, a rendezvous, a rendezvous with the moon, or a deep space mission which will include rendezvous with other planets, solar/celestial body orbiting, and gravity assists. Depending on the choice, the program will either prompt for the target orbit's orbital elements, target satellite/spacecraft's orbital elements, or plot where the moon is.

Left: General orbit transfer: Plot of Earth, current trajectory and location, and target orbit.

Center: Rendezvous: Plot of Earth, current trajectory and location, and target orbit and satellite location.

Right: Plot of Earth, current trajectory and location, and current trajectory and location of the Moon.

This is as far as I have written. Next the user will decide at what anomaly they wish to launch, anomaly to arrive at target orbit, and the time of the transfer. The orbital transfer, rendezvous, and moon rendezvous are going to have very similar coding so I plan on finishing those before beginning the deep space function.

There are currently 27 functions (soon to be more), so needless to say my folder is beginning to become clustered. Once I finish what I am currently working on, or when I feel like it, I am going to go through and delete redundant code and re-write some functions to increase versatility.