// $Id: sem.C,v 1.3 2006/12/05 08:59:12 jingliu Exp $

// include std C++ IO functions
#include <iostream>
#include <fstream>
using namespace std;

// include necessary header files provided by ROOT (the declarations
// of Int_t, Double_t, etc. are automatically included when you
// include one of these header files, eg. TMath.h)
#include <TMath.h>
using namespace TMath;

#include <TVector3.h>   // http://root.cern.ch/root/html/TVector3.html
#include <TCanvas.h>
#include <TGraph.h>
#include <TLegend.h>
 
const Int_t n = 3;              // Sun, Earth and Moon   

// change rate of the position
TVector3 fx(Int_t i, TVector3 *x, TVector3 *v)
{
   return v[i];
}

// change rate of the velocity
TVector3 fv(Int_t i, TVector3 *x, TVector3 *v)
{
   Double_t GM[n];              // G*M
   GM[0] = 4.0*Power(Pi(),2);   // mass of the Sun
   GM[1] = GM[0]*3.00e-6;       // mass of the Earth
   GM[2] = GM[1]*1.23e-2;       // mass of the Moon
   
   // by default all vector components are initialized to zero
   TVector3 r;                  // relative position vector
   TVector3 a;                  // acceleration
   
   for (Int_t j=0; j<n; j++) {
      if (j != i) {
         r = x[i]-x[j];
         a+=-GM[j]/r.Mag2()*r.Unit();
      }
   }
   return a;
}

// std C++ main function, please refer to
// http://wwwx.cs.unc.edu/~sparkst/howto/cpp_main.php
//int main(int argc, char *argv[])
int sem(Double_t dt=0.001, Double_t Tmax=1)
{
   // read arguments and initialize dt, Tmax accordingly
   Double_t dt, Tmax;   
/*
   if (argc == 2) { // the first argument, argv[0], is the name of the program
      dt   = atof(argv[1]); // http://www.cppreference.com/stdstring/index.html
      Tmax = 1;                 // 1 year
   }
   else if (argc >= 3) {
      dt   = atof(argv[1]);
      Tmax = atof(argv[2]);
   }
   else {
      dt   = 0.001;
      Tmax = 1;
   }
*/
   ofstream dat[n];
   dat[0].open("sun.dat");
   dat[1].open("earth.dat");
   dat[2].open("moon.dat");
   
   // initialize system
   Double_t t = 0;
   Double_t Ve=6*3.14159265;    // velocity of the Earth
   Double_t R=0.985;            // R = V_moon/V_earth 

   TVector3 x[n], v[n];         // position & velocity vectors
   TVector3 xp[n], vp[n];       // temporary vectors for calculation

   x[0].SetXYZ(0,0,0);          // Sun in Origin
   x[1].SetXYZ(0.2,0,0);        // Earth is 0.2 A.U. away from Sun (periapsis)
   x[2].SetXYZ(0.2026,0,0);     // Moon is 0.0026 A.U. away from Earth

   v[0].SetXYZ(0,0,0);          // Sun is still
   v[1].SetXYZ(0,Ve,0);         // Earth moves towards Y direction
   v[2].SetXYZ(0,Ve*R,0);       // Moon moves also towards Y direction
   
   while (t<Tmax) {             // main loop
      for (Int_t i=0; i<n; i++) // output results
         dat[i]<<x[i].X()<<" "<<x[i].Y()<<endl;
      
      // 2nd order Runge-Kutta Method
      for (Int_t i=0; i<n; i++) { // loop over sun, earth, moon
         xp[i] = x[i] + dt/2.0*fx(i,x,v);
         vp[i] = v[i] + dt/2.0*fv(i,x,v);
      }
                                
      for (Int_t i=0; i<n; i++) {
         x[i] += dt*fx(i,xp,vp); // dt*fx works well, fx*dt doesn't work
         v[i] += dt*fv(i,xp,vp);
      }
      t += dt;
   }

   for (Int_t i=0; i<n; i++) dat[i].close();

   // draw the results
   TCanvas *c1 = new TCanvas("c1","",0,0,800,800);
   TGraph *g0 = new TGraph("sun.dat");
   TGraph *g1 = new TGraph("earth.dat");
   TGraph *g2 = new TGraph("moon.dat");
   
   g2->SetTitle("the Sun, Earth and Moon");
   g2->GetXaxis()->SetTitle("X (A.U.)");
   g2->GetXaxis()->CenterTitle();
   g2->GetYaxis()->SetTitle("Y (A.U.)");
   g2->GetYaxis()->CenterTitle();

   g2->GetXaxis()->SetLimits(-2,0.5); // set real range of axis
   g2->GetXaxis()->SetRangeUser(-2,0.5); // set visible range of axis
   g2->GetYaxis()->SetLimits(-1.25,1.25);
   g2->GetYaxis()->SetRangeUser(-1.25,1.25);

   g2->SetMarkerStyle(20); // http://root.cern.ch/root/html/TAttMarker.html
   g2->SetMarkerSize(0.2);
   g2->SetMarkerColor(46); // http://root.cern.ch/root/html/TAttFill.html
   g2->Draw("ap");
   
   g1->SetMarkerStyle(20);
   g1->SetMarkerSize(0.2);
   g1->SetMarkerColor(38);
   g1->Draw("p");
   
   g0->SetMarkerStyle(29);
   g0->SetMarkerSize(1);
   g0->SetMarkerColor(2);
   g0->Draw("p");

   TLegend* leg = new TLegend(0.6, 0.7, 0.85, 0.85);
   leg->AddEntry(g0,"Sun","p");
   leg->AddEntry(g1,"Earth","p");
   leg->AddEntry(g2,"Moon","p");
   leg->Draw();

   // http://root.cern.ch/root/html/TPostScript.html#TPostScript:description
   c1->Print("sem.pdf");
 
   return 0;                    // C++ std
}