MatLab Programs (Under Construction...)

Description: This is a program writen by Luigi Giaccari which computes triangulations (by two simplices) of point cloud data in three-dimensions. Very useful for plotting two dimensional surfaces given by point data.

Matlab Central Page

Author's Page(Find other cool matlab programs here as well.)

This first collection of programs are associated with the introduction note set for the N-Body problem.

nBodyWpar.m: This is the main files needed for the simulation of the N-Body problem. This file codes the vector field for the problem as a matLab function. The number of bodies, as well as their masses and the gravitational constant must be passed along with the current position, to the function.

nBodyPlayGround.m: A test program for the program 'nBodyWpar.m'. Before running this program set the number of bodies 'N', the gravitational constant 'G', and the length of the simulation 'simulationTime'. The program then generates random masses and initial conditions and integrates the resulting system. The output is several graphs illustrating the results. The program is fun to play with but also demonstrates that 'nBodyWpar.m' works properly with any number of bodies and parameter values that may occur in other work.

NbodyOne.m:: The program is used in the introduction set of notes to look at the five body problem with symmetric initial conditions. It calls the program 'rstTOijk.m' to carry out some coordinate transformations described in those notes, and the program 'nBodyWpar.m' for the integration.

rstTOijk.m: This function computes a coordinate transformation used in NbodyOne.m.

twobodyDesign.m: The program is used in the introduction note set to design elliptic orbits in the two body problem with desired properties (eccentricity ect).

sitnikovMap.m: The program simulates the Sitnikov problem (actually the continuation of the Sitnikov problem into the three body problem, but with very small perturbation parameter). It is used to produce the results in the introduction set of notes that pertain to the Sitnikov problem. The output is a couple of Poincare maps and plots of the height of the third body and the configuration of the entire system in physical space (R^3). The program lets the user define both the number of initial conditions to be included in the map, and the number of intersections of each trajectory with the Poincare section.

NOTE: This program replaces the program 'sitnikov.m' which use to be posted here. That program had several errors and more importantly did not allow the number of crossings to be set at run time. One simpley integrated for 'a long time' and hoped to find enough.

sitnikovNewton.m: This function is used in the program 'sitnikovMap.m'. The caller passes a pair of points, one on each side of the Poincare section where neither are with in the prescribed tolerence which would allow them to be considered intersections. The function implements a Newton algorithm which takes this as the 'initial guess' and then converges to a true intersection between them.

Under construction




This collection of programs are associated with the note set 4. They can be thought of as a general MatLab toolbox for the CRTBP.

CRTBP.m: This function codes the vector field for the CRTBP. Any program which integrates the CRTBP in forward time will call this.

BackwardCRTBP.m: This is the vectorfield for reversing time in the CRTBP. This function and the previous one are the essential functions required to integrate the CRTBP.

librationPoints.m: Given the value of mu for the system, the fucntion computes the coordinates of all five libration points.

jacobiConst.m: When passed a value of mu and a state of the CRTBP, the function computes the Jacobi constant of the state.

magnitudeVelocity.m: When passed a value of mu, a position and the value of the Jacobi integral, the function computes the magnitude of the corrosponding velocity vector. This is typically used when computing initial velocities where you know all of the above and what direction you want to go. This information determines uniquely the magnitude of the velocity.

stateTransCRTBP.m: Given and initial time, a final time, a state of the CRTBP (reference orbit) and the value of mu, the function computes the State transition matrix (differential of the flow). This function calls the next two functions as well.

sysSolveCRTBP.m: The vector field for the equations of first variation for the CRTBP. Necessary for numerically computing the differential of the flow.

G_CRTBP.m: This function computes the submatrix 'G' for the CRTBP. (See note set 2, 4, and 5 for details).

vectorField_CRTBP.m: When passed a state of the CRTBP this function computes the vector field at that point. Since this is also the time derivative of the flow at that state/time it is sometimes a usefull function to have handy. (Function is almost identical to CRTBP, except that no integration information is passed along, which makes it cleaner in application).

CRTBPpoincareNewton.m: This function is handed , the value of mu, a tolerence, two points on either side of the y=0 plane, and the times at which they occur. It uses a Newton Method to find the crossing to a number of decimal places specified by the variable 'tollerence'. The name comes from the fact that it is a Newton Method for computing Poincare section crossings. All the poincare section programs call this. (Actually the sections are computed in the Planar CRTBP as this is a problem with a two dimensional section after the value of the energy has been fixed).

interpCrossingData.m: When the third body gets too close to the primaries then the Newton method above can really slow down. When this happens some of the programs below use this function to simply interpolate the crossing. This is a quadradic interpolation which matches the derivatives at the points to a porabola. However when this happens there is no garuntee that the tolerence is met, so it's good to keep count of how many times this happens in programs.

The functions above are the components of the programs that follow. The programs that follow use the tool kit above to explore various aspects of, and solve problems concerning, the CRTBP.

CRTBPplayGround.m: The user specifies a time (stopping time for the integration) and the program generates a random mu, and random initial conditions. The CRTBP is integrated untill the specified time from the intial conditions.
The program is fun to play with but also tests the CRTBP tool box on randomly chosen data.

CRTBPpoincareGrid.m: This is a program which computes a Poincare section for the planar CRTBP. The section is computed on a grid of points instead of for points only along the x-axis. This 'fills in' a lot of the 'holes' that one sometimes sees in this kind of picture (such holes are cuased by trajectories on invariant tori that never cross the x axis with xdot=0). If one wants to build the xdot=0 map then this program could be run with a grid of one line. The program could also be used to build maps that refine certian regions in the Poincare maps domain.

The program takes a very long time to run. A day or two on my old desktop. If you run it on a UNIX system I suggest doing it in background mode. I also have a link below to the data from a run that I did. If one is really serious about this kind of thing, then this is about the place where you switch from MatLab to say FORTRAN or C++. I hope to have some such code up here soon along with runtime comparisons...

gridData1.mat(5.8meg): This is a large set of data I got out from the program above. The program found 1000 crossings of 100 different trajectories, four different times. The first time was the usual run with initial conditions on the x axis with xdot=0. The rest were done with different starting values of xdot>0. This filled in the large holes nicely. (The data comes from the default settings of the program 'CRTBPpoincareGrid.m'. Then the value of the Jacobi constant and other information that goes with the picture can be found there).

To see the plot, download that file (which is 5.8 megabites) and put it in your MatLab command path (as though it were an m-file). Once this is done type:

load gridData1;

at the command line prompt. After the prompt returns, type 'who' (withought the quotes). You should see that you have a variable called 'PoincareGrid' in the workspace (along with a dozen or so other variables). Now just type:


at the prompt. The resulting plot will be a rather dull looking solid shape. However if you experiment zooming in on it you will see lots of structure. Things are pretty busy for 0.1< x< 0.2. There is also a big resonance torus near x = 0.5. 'Static' in the picture is due to chaotic behavior. (See note set 4).

trajectoryDoctor.m: This is a companion program to the Poincare section programs. Pick a point on the Poincare section (an x and xdot pair. Just zoom into the picture and eye ball these as well as you can). Enter this data to 'trajectoryDoctor.m' along with the value of the Jacobi energy that was used to generate the Poincare Section, and how many crossings you want to see. 'trajectoryDoctor' will plot the orbit in configuration space for you, and show the Poincare section of just the orbit you asked about. This lets you see what certin structures in the Poincare section correspond to in configuration space.

This collection of programs are associated with the note set 5. They compute halo orbits, families of halo orbits, and stable and unstable manifolds of unstable periodic orbits.

unstableEarthMoonLyap.m: This program computes and plots the unstable manifold of an unstable periodic orbit. The branch that is plotted is the branch on the side of the second primary. If the default settings are not changed then the computation is caried out for a Lyapunov orbit near L1, in the Earth/Moon system.

If the user wishes to compute with a different orbit, then the initial conditions, and the period of the orbit must be input by the user. The program will compute the unstable manifold.

The program plots orbit segments of the manifold as it goes but does not keep the data. Then, the program is primarialy a visuilation tool as writen (This could easily be changed).

NOTE: If you run this program be carefull to read about the dependencies in the notes. The program calls several programs from the forth set of notes

stableEarthMoonLyap.m: This program computes and plots the stable manifold of an unstable periodic orbit. The branch that is plotted is the branch on the side of the second primary. If the default settings are not changed then the computation is caried out for a Lyapunov orbit near L1, in the Earth/Moon system. This is a slight modification of the program 'unstableEarthMoonLyap.m'. One could make a plot showing both the stabel and unstable manifolds by either combining the two prograrams or saving the data from each.

EarthSunHaloOrbit_NewtonMewhod.m: This program computes a halo orbit near the Sun/Earth L1 libration point, from an initial guess given in metric units. (The initial guess was provided by the instructor). This is the code that is used to obtain the results from section 3.1 in note set 5 (there the program was misnamed haloNewton3).

haloFamily.m: This program continues the Halo Orbit found in the previous program, computing the cylinder of Halos parametrized by energy. Used in 3.2 of note set 5.

haloFamilyL2.m: Begining from a planar Lyapunov orbit near L2 which is near bifurcation, the program continues the orbit out of plane, finding the cylinder of halo orbits. This is one way to compute halo orbits with no prior knowlege of their location. The program is used in section 3.3 of note set 5. (This program and the previous one take around an hour to run, depending on the number of orbits in the continuations, i.e the number of steps).


Description: A number of programs for computing stable and unstable manifolds in the standard map. The purpose of the programs is to illustrate the parametrization method and compare it with first order methods.

standardMap.m: M-File computes one iteration of the standard map.

inverse_standardMap.m: M-File computes one inverseiteration of the standard map.

standardMapSimulationII.m: Program produces a phase space plot of the standard map, and computes the stable and unstable manifolds by linear analysis.

parmEval.m: When passed a time, the matrix coefficents of the parametrizations, the desired order of the approximation, computes the polynomial at the given time.

derivEval.m: Similar to above but computes the tangent vector to the manifold.

homoclinicProof.m: Computes an Intersection of the stable and unstable manifolds by a Newton method.

parMethWS_and_WU.m: Computes the stable and unstable manifolds by parametrization method. User can set the order and parameter interval.

errorAnalysis.m: Error analysis of parametrization.

parMethVarOrder.m: An early version of the parametrization code. Just computes the unstable manifold.


Description: MatLab Tutorial for Financial Math REU summer school.
Financial math and matlab.pdf

Here is a link a good matlab reference page.
MatLab Tutorial

Description: Some more easy matlab examples.

Here is a link to some more mathematical biology matlab codes.
Mathematical Models in Biology

back to my homepage

back to celestial mechanics page