## Worked solution for PLANET example

Planet.cc is a worked solution for the task:
 Write code that implements the Leapfrog method. Get it going for one initial condition, such as x0 = (1.5,2), v0 = (0.4,0).

Every few iterations it prints out
t, x1, x2, v1, v2, T, V, T+V
where T is the kinetic energy and V is the potential energy. When it's printed all those outputs to a file, you can plot things from the file using gnuplot.
Directories: The first C++ program, Planet.cc, will put its output in a directory called output (provided you first create that directory and then run the program following the instructions written in the comments at the top of Planet.cc). The "gnu" files (gnu | gnu10 | gnu11) go in the upper directory.
```// Planet.cc
//  Simulates planet (inverse square law)
//    using leapfrog method.
//  example usage:
//       mkdir output
//       ./Planet > output/1
// recommended gnuplot commands:
//  set size ratio -1
//  plot   'output/1' u 2:3 w linesp lt 3 pt 4
//  replot 'output/1' u 2:3:4:5 w vector lt 5

#include <iostream>
#include <cmath>
using namespace std;

#define D  2  // number of dimensions
struct particle {
double x[D] ; // (x,y) coordinates
double p[D] ; // momentum
double F[D] ; // force
double im   ; // inverse mass
double GMm  ; // gravitational parameter of this particle
double v[D] ; // velocity
double T    ; // kinetic energy
double V    ; // potential energy
double r    ; // absolute distance from origin
} ; // Note the definition of a structure ends with a semicolon

struct control {
int verbose ; // program verbosity
int printing ; // period with which to print
};

void Force( particle &a )
// sets the force vector and the potential energy
{
double R = 0.0 ;
for ( int i = 0 ; i < D ; i++ ) {
R += a.x[i] * a.x[i] ;
}
double r = sqrt(R) ;
a.V = - a.GMm / r ;
a.r = r ;
for ( int i = 0 ; i < D ; i++ ) {
a.F[i] = - a.GMm * a.x[i] / (r*R) ;  // inverse sq
}
}

void PositionStep ( particle &a , double dt )
{ // the "&" indicates 'this is a pass by reference'
for ( int i = 0 ; i < D ; i++ ) {
a.x[i] += dt * a.p[i] * a.im ;
}
}

void MomentumStep ( particle &a , double dt )
{ // the "&" indicates 'this is a pass by reference.
for ( int i = 0 ; i < D ; i++ ) {
a.p[i] += dt * a.F[i] ;
}
}

void v2p( particle &a )
// propagate a change in velocity into the momentum vector
{
for ( int i = 0 ; i < D ; i++ ) {
a.p[i] =  a.v[i] / a.im ;
}
}

void p2v ( particle &a ) {
for ( int i = 0 ; i < D ; i++ ) {
a.v[i]  = a.p[i] * a.im ;
}
}

void pv2T ( particle &a ) {
a.T=0.0;
for ( int i = 0 ; i < D ; i++ ) {
a.T    += 0.5*a.v[i] * a.p[i] ;
}
}

void showState ( particle &a )
{
for ( int i = 0 ; i < D ; i++ ) {
cout << "\t"<<a.x[i];
}
for ( int i = 0 ; i < D ; i++ ) {
cout << "\t"<<a.v[i];
}
cout << "\t" << a.T << "\t" << a.V << "\t" << a.T+a.V << endl;
cout << "#" ;
for ( int i = 0 ; i < D ; i++ ) {
cout << "\tx["<<i<<"]\t";
}
for ( int i = 0 ; i < D ; i++ ) {
cout << "\tv["<<i<<"]\t";
}
cout << "\tT\t\tV\t\tT+V" ;
cout << endl;
}

void LeapfrogDynamics( particle &a ,
double dt , // step size
double &t , // current time
int N ,  // number of steps
control &c  // parameters controlling output, etc
)
{
for ( int i=0 ; i < N ; i ++ ) {
if( c.printing && !(i%c.printing) ) {
// this iteration we will print out the state
p2v(a);   pv2T(a) ;     Force( a ) ;
cout << t ; showState( a ) ;
}
// each iteration we move the position a half-step
PositionStep( a , dt*0.5 ) ;
Force( a ) ; // (computes the force at this position)
// then we move the momentum full-step
MomentumStep( a , dt )     ;
// then we move the position a half-step
PositionStep( a , dt*0.5 ) ;
t += dt ;
}
}

void showVelocity ( particle &a )
{ // the "&" indicates 'this is a pass by reference.
for ( int i = 0 ; i < D ; i++ ) {
cout << "v["<<i<<"] = "<<a.v[i]<<endl;
}
}

int main()
{
particle   a ;
control    c ;
c.verbose=1;
c.printing=5 ; // how many iterations to go
//               between state-printings

a.im   = 1.0 ;
a.v[0] = 0.4;
a.v[1] = 0.0;
a.x[0] = 1.5;
a.x[1] = 2;

a.GMm = 1.0;
v2p(a) ; // make sure the momentum is set up

double dt = 0.011 ; // time per step
double t = 0.0 ; // start time
int N = 1000 ;  // number of steps to make

LeapfrogDynamics( a , dt , t , N , c ) ;

return 0;
}

```

gnu is a simple gnuplot file for plotting the output of Planet.cc. It displays the positions and the velocities in a single graph.

```set size ratio -1
set nokey
r = 2
r2= 2.5
set xrange [-r:r2]
set yrange [-r:r2]
plot 'output/sun' w p ps 3 pt 5, \
'output/1' u 2:3 w l lt 2

plot 'output/sun' w p ps 3 pt 5, \
'output/1' u 2:3 w l lt 2, \
'output/1' every 10 u 2:3 w p ps 4 pt 3

plot 'output/sun' w p ps 3 pt 5, \
'output/1' u 2:3 w l lt 2,\
'output/1' every 20 u 2:3:4:5 w vector lt 5

plot 'output/sun' w p ps 3 pt 5, \
'output/1' u 2:3 w l lt 2,\
'output/1' every 2 u 2:3:4:5 w vector lt 6, \
'output/1' every 10 u 2:3 w p ps 4 pt 4

pause -1 'and finally here is the energy versus time'

set autos xy
set size noratio
plot \
'output/1' u 1:8 w l lt 2

```

gnu10 and gnu11 are a fancier pair of gnuplot files for making a movie from the output of Planet.cc. gnu10 sets things up, then calls gnu11, which does the plot and calls itself recursively.

```# This gnuplot movie
# shows six evolving snapshots
set size ratio -1
set nokey
r = 4
r2= 4.5
set xrange [-r:r2]
set yrange [-r:r2]
# time range for line
tmin1 = 0
tmax1 = 100
# time range for points
tmin2 = 50
tmax2 = 100
# gap between points
gap2 = 10
step = 3
LastLine = 200
call 'gnu11'
```

```plot 'output/sun' w p ps 3 pt 5, \
'output/1' every gap2::tmin2::tmax2 u 2:3 w p ps 4 pt 4,\
'output/1' every 1::tmin1::tmax1 u 2:3 w l lt 2
pause 0.1
tmin1=tmin1+step
tmin2=tmin2+step
tmax1=tmax1+step
tmax2=tmax2+step
if ( tmax1 < LastLine ) call 'gnu11'
```

The file sun contains the coordinates of the sun

 Planet2.cc is a much fancier solution for the task. It can handle either of two force laws (inverse square and Hooke law); and while it is running, it can print out gnuplot commands to stdout, so if you pipe its output into gnuplot, it makes a live movie. It can be modified so that it performs an Euler simulation immediately after a Leapfrog one. These extra features are not very elegant. Directories: For the second program (Planet2.cc), "output/gnu0" should go in the output directory. The output of the program is written to files in the "/tmp" folder, and commands for gnuplot are written to standard output (stdout). See the head of the cc file for instructions. ```// Planet2.cc // see also gnu20, gnu21 // usage: // Planet2 1300 0.03 25 | gnuplot -persist // This program requires the directory 'output' to exist, and // expects it to contain a file gnu0. #include #include #include #include using namespace std; #define D 2 // number of dimensions struct particle { double x[D] ; // (x,y) coordinates double p[D] ; // momentum double F[D] ; // force on particle double im ; // inverse mass of particle double GMm ; // gravitational parameter of particle double v[D] ; // velocity double T ; // kinetic energy double J ; // angular momentum double V ; // potential energy double r ; // absolute distance from origin int law ; // which force law } ; // Note the definition of a structure ends with a semicolon struct control { int verbose ; // program verbosity int printing ; // period with which to print int commanding ; // period with which to issue commands // ofstream *fout ; // where we write data to int time_in_ms ; // time between frames (commands) clock_t latest_time; ofstream fout ; int euler ; // whether to do euler method after leapfrog }; // the "&" in "&a" indicates 'this is a pass by reference'. void Force( particle &a ) // sets the force vector and the potential energy { double R = 0.0 ; for ( int i = 0 ; i < D ; i++ ) { R += a.x[i] * a.x[i] ; } double r = sqrt(R) ; a.r = r ; if ( a.law == -2 ) { // inverse square law a.V = - a.GMm / r ; for ( int i = 0 ; i < D ; i++ ) { a.F[i] = (- a.GMm * a.x[i] / (r*R) ) ; // inverse sq } } else { // Hooke Law a.V = 0.5 * a.GMm * R ; for ( int i = 0 ; i < D ; i++ ) { a.F[i] = (- a.GMm * a.x[i]) ; // linear spring law } } } // Implement x += v dt void PositionStep ( particle &a , double dt ) { for ( int i = 0 ; i < D ; i++ ) a.x[i] += dt * a.p[i] * a.im ; } // Implement p += f dt void MomentumStep ( particle &a , double dt ) { for ( int i = 0 ; i < D ; i++ ) a.p[i] += dt * a.F[i] ; } void v2p( particle &a ) // propagate the changed velocity into the momentum vector { for ( int i = 0 ; i < D ; i++ ) a.p[i] = a.v[i] / a.im ; } void p2v ( particle &a ) { for ( int i = 0 ; i < D ; i++ ) a.v[i] = a.p[i] * a.im ; } void pv2T ( particle &a ) { // compute the kinetic energy and angular momentum a.T=0.0; for ( int i = 0 ; i < D ; i++ ) { a.T += 0.5*a.v[i] * a.p[i] ; } a.J = a.x[0] * a.p[1] - a.x[1] * a.p[0] ; } // Print out the state (x,v,Energies,J) in columns void showState ( particle &a , ostream &fout ) { for ( int i = 0 ; i < D ; i++ ) { fout << "\t"< /tmp/2"); showStateByArrow( a , cout ) ; // sends 'set arrow' command to gnuplot if(i==0) { cout << "load 'output/gnu0'" << endl ; // to gnuplot cout.flush() ; c.latest_time = clock() ; } else { c.latest_time += CLOCKS_PER_SEC * c.time_in_ms / 1000 ; clock_t endwait = c.latest_time ; clock_t time_to_go = endwait-clock() ; while( clock() < endwait ) {} ; // wait till the appointed time cout << "load 'output/gnu0'" << endl ; // to gnuplot // cout << "replot" << endl ; while( clock() < endwait+time_to_go/3 ) {} ; // add a little further // pause to give gnuplot a chance to do its thing before we edit // files again } } } // The Leapfrog method simulates the equations of motion // dx/dt = p/m // dp/dt = f // by alternately updating the position (using the latest // momentum), then the momentum (using the force at the // updated position). void LeapfrogDynamics( particle &a, double dt , // step size double &t , // current time int N , // number of steps control &c // parameters controlling output, etc ) { for ( int i=0 ; i < N ; i ++ ) { PossibleOutput( a , c , i , t ) ; // each iteration we move the position a half-step PositionStep( a , dt*0.5 ) ; Force( a ) ; // (computes the force at this position) // then we move the momentum full-step MomentumStep( a , dt ) ; // then we move the position a half-step PositionStep( a , dt*0.5 ) ; t += dt ; } } // The Euler method updates the position and momentum // simultaneously. For the force laws considered here, // we can get the effect of simultaneous update by // finding the force, then // updating the position, THEN the momentum. void EulerDynamics( particle &a, double dt , // step size double &t , // current time int N , // number of steps control &c // parameters controlling output, etc ) { for ( int i=0 ; i < N ; i ++ ) { PossibleOutput( a , c , i , t ) ; Force( a ) ; // (computes the force at this position) PositionStep( a , dt ) ; MomentumStep( a , dt ) ; t += dt ; } } int main(int argc, char* argv[]) { particle a ; control c ; // set defaults int N = 250 ; double dt = 0.1 ; double t = 0.0 ; c.verbose = 1 ; c.printing = 4 ; c.commanding = 1 ; // number of iterations between plots c.time_in_ms = 100; // real time between plots (in ms) c.euler = 0 ; // whether to do euler after leapfrog // read in any command-line arguments if(argc>1) { sscanf( argv[1], "%d", &N ) ; // put the first command-line argument in N } if(argc>2) { sscanf( argv[2], "%lf", &dt ) ; // put the 2nd argument in dt } if(argc>3) { sscanf( argv[3], "%d", &(c.time_in_ms) ) ; } // try to write output file char fileName[]="/tmp/1"; c.fout.open(fileName); if(c.fout.good()==false){ cerr << "can't write to file " << fileName << endl; exit(0); } // set initial conditions for a particle a.law = -2 ; // -2: inverse square a.im = 1.0 ; a.v[0] = 0.4; a.v[1] = 0.0; a.x[0] = 1.5; a.x[1] = 2; a.GMm = 1.0; v2p(a) ; // get the initial momentum from the velocity LeapfrogDynamics( a , dt , t , N , c ) ; if ( c.euler ) { cerr << "# SWITCHING TO EULER\n" ; c.fout << "\n\n\n" ; EulerDynamics( a , dt , t , N , c ) ; } return 0; } ```

David MacKay