Computational Physics, 1B, Term 2

Worked solution for PLANET example

| README | Planet.cc | gnu | gnu10 | gnu11 | output/sun | Planet2.cc | gnu20 | gnu21 | output/gnu0 |

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

pause -1 'ready'

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

pause -1 'ready'
   
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
 
pause -1 'ready'
   
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


example1

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 <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib>
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"<<a.x[i];
  }
  for ( int i = 0 ; i < D ; i++ ) { 
    fout << "\t"<<a.v[i];
  }
  fout << "\t" << a.T     // KE
       << "\t" << a.V     // PE
       << "\t" << a.T+a.V // Total E
       << "\t" << a.J     // Angular momentum
       << endl;
  fout << "#" ;
  for ( int i = 0 ; i < D ; i++ ) { 
    fout << "\tx["<<i<<"]\t";
  }
  for ( int i = 0 ; i < D ; i++ ) { 
    fout << "\tv["<<i<<"]\t";
  }
  fout << "\tT\t\tV\t\tT+V\t\tJ" ; 
  fout << endl;
}

void showStateByArrow ( particle &a , ostream &fout )
{
  fout << "set arrow 101 from " << a.x[0] << ","
       << a.x[1]  << " to "
       << a.x[0]+a.v[0] << ","
       << a.x[1]+a.v[1] << " lt 5 " << endl ; 
}

void PossibleOutput( particle &a , control &c , int i , double t ) {
  if( ( c.printing && !(i%c.printing) )
      || ( c.commanding  && !(i%c.commanding) ) ) {
    p2v(a) ;  pv2T(a) ;  Force( a ) ;
  }
  if( ( c.printing && !(i%c.printing) ) ) {
    c.fout << t ; showState( a ,  c.fout ) ;
  }
  if ( c.commanding  && !(i%c.commanding) ) {
    // SEND COMMANDS TO gnuplot live via 'cout' pipe
    c.fout.flush() ; // make sure file is out
    system("tail -39 /tmp/1 > /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
Last modified: Mon Sep 24 16:06:61 2012