Computational Physics, 1B, Term 2

Worked solution for PLANET example USING CLASSES

| README | planet.cc | gnu | sun | all.zip |

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).

It is written in advanced C++, going beyond what's expected in the 1B course. It defines classes (see page 61 of the first tutorial) for particles. It uses a structure called xyvect to hold 2-dimensional vectors. This structure definition defines how to do arithmetic with these objects, thus permitting the simulation code to contain elegant statements such as
    x += v * dt;
and
    v += a * dt;
Such human-friendly statements are not possible in many older languages such as C.
In the class approach to the particle, all the functions associated with the particle are put within the definition of the particle.
The program writes the trajectory to stdout, so the intended usage is something like:
 planet > tmp
The printout contains
t, x1, x2, v1, v2, r, theta, J, 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.
For example
gnuplot
load 'gnu'
 


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

# planet > tmp
 set style line 5 lw 0.5 lt 5
 set style line 1 lw 0.5 lt 1
 set size ratio -1
 set nokey
 plot 'sun' u 1:2 w p ps 3 pt 3 , \
  'tmp' u 2:3 w linesp ps 4 pt 4

 pause -1 'ready'
   
 plot 'sun' u 1:2 w p ps 3 pt 3 , \
  'tmp' u 2:3 w linesp ps 3 pt 4,\
  'tmp' u 2:3:4:5 w vector lt 5

 pause -1 'ready'
 plot 'sun' u 1:2 w p ps 3 pt 3 , \
   'tmp' u 2:3 w linesp ps 3 pt 4,\
   'tmp' u 2:3:4:5 w vector ls 5,\
   'tmp' every 1::::1 u 2:3:4:5 w vector lt 1
# the 'every' above selects just the first line 
# of the file and displays a vector with line type 1


			     

The file sun contains the coordinates of the sun
0 0
// planet.cc     by George Saklatvala 
//
// Shows how to implement dynamics
// in an elegant advanced C++ way.
//
// The structure xyvect, declared at
// the top of the file, uses arithmetic operator overloading
//
// Usage:
//   planet > tmp
//   set size ratio -1 ; plot 'tmp' u 2:3 w linesp 3 4, 'tmp' u 2:3:4:5 w vector lt 5

#include<iostream>
#include<fstream>

#include<cfloat>
#include<cstdlib>
#include<cmath>

using namespace std;

struct xyvect{

  double x,y;
  xyvect(): x(), y() {}
  xyvect(double xx, double yy): x(xx), y(yy) {}
  double normsq() const {return x*x + y*y;}
  double norm() const {return sqrt(normsq());}
  double ang() const {return atan2(y,x);}

  xyvect& operator+=(const xyvect& r)
  {
    x += r.x;
    y += r.y;
    return *this;
  }

  xyvect& operator-=(const xyvect& r)
  {
    x -= r.x;
    y -= r.y;
    return *this;
  }

  xyvect& operator*=(double r)
  {
    x *= r;
    y *= r;
    return *this;
  }

  xyvect& operator/=(double r)
  {
    x /= r;
    y /= r;
    return *this;
  }

};

xyvect operator* (double d, const xyvect& v)
{
  xyvect a(v);
  a*=d;
  return a;
}

xyvect operator* (const xyvect& v, double d)
{
  xyvect a(v);
  a*=d;
  return a;
}

xyvect operator/ (const xyvect& v, double d)
{
  xyvect a(v);
  a/=d;
  return a;
}

struct impact{}; //exception class

// This is where the planet bit starts....

const double G = 1; // gravitational constant
const double M = 1; // mass of sun

class particle{

private:

  const double m;
  xyvect x;
  xyvect v;
  double t;
  xyvect a;

  double r() const {return x.norm();}
  double theta() const {return x.ang();}
  double J() const {return m*(x.x*v.y-x.y*v.x);}
  double T() const {return m*v.normsq();}
  double V() const {return -G*M*m/r();}
  double E() const {return T() + V();}

public:

  particle(double mass, xyvect x0, xyvect v0, double t0=0): 
    m(mass), x(x0), v(v0), t(t0), a(){}

  void x_inc(double dt) // step x in the direction of v.
  {
    x += v * dt;
    if(r()<2*dt*v.norm() || r() < DBL_EPSILON)throw impact();
  }

  void v_inc(double dt) // step v in the direction of the acceleration.
  {
    v += a * dt;
  }

  void t_inc(double dt)
  {
    t += dt;
  }

  void a_calc() // find the acceleration
  {
    a = -G*M*x/pow(r(),3);
  }

  ostream& write_state(ostream& of) const
  {
    of<<t<<" "
      <<x.x<<" "<<x.y<<" "
      <<v.x<<" "<<v.y<<" "
      <<r()<<" "<<theta()<<" "
      <<J()<<" "<<T()<<" "
      <<V()<<" "<<E()<<endl;
    return of;  
  }
};

void euler_step(particle& p, double dt)
{
  p.a_calc();
  p.x_inc(dt);
  p.v_inc(dt);
  p.t_inc(dt);
}

void leapfrog_step(particle& p, double dt)
{
  p.x_inc(0.5*dt);
  p.a_calc();
  p.v_inc(dt);
  p.x_inc(0.5*dt);
  p.t_inc(dt);
}

void simulate(particle& p, void(*step)(particle&, double),
	      double t, int n, int frames, ostream& of)
{ // number of steps made is n (total duration is t)
  // number of frames written out is frames
  const double dt = t / n;
  for(int i=0;i<frames-1;++i){
    p.write_state(of);
    for(int j=0;j<n/frames;++j)
      step(p,dt);
  }
}

int main()
{
  double time = 11.0; // time to simulate
  int n = 1000;      // total number of steps to make
  int frames = 100; // number of frames to record
  double m = 1;
  const xyvect x0(1.5,2.0);
  const xyvect v0(0.4,0);
  particle p(m,x0,v0);
  try{
    simulate(p,leapfrog_step,time,n,frames,cout);
    // edit the above command to change the step style
  }catch(impact&){
    cout << "Collision occurred\n";
  }
  return 0; 
}
example


David MacKay
Last modified: Wed Nov 14 10:11:32 2007