|
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:
plot(PoincareGrid(:,1),PoincareGrid(:,2),'b.')
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.
|