#ifndef _CINT_
// mit CERN-Grafik: siehe Mathematik-Pool, X4240-A.math.tu-cottbus.de
//     unter /opt/cern oder
//     http://root.cern.ch/drupal/ (fuer Download und Dokumentation)
//     1.) Setze die Umgebungsvariablen: 'source /opt/cern/bin/thisroot.sh'
//         (nur am Beginn einer Sitzung erforderlich!)
//     2.) Uebersetzen des C++-Programms, z.B.
//         'g++ `root-config --cflags` grxx011.cpp \
//                            `root-config --glibs` -o grxx011 -lm'
//     3.) bei Verwendung von NR: -I${NR_INC}, -L${NR_LIB_FLAGS}, -lnmr
//
#include <TApplication.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TF1.h>
#endif

#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
#include <cmath>
#include "nr.h"
using namespace std;

#define PI 3.141592653589793
#define N 4         // Dimension of the system (1-st order)
#define NP 40       // Number of discretization points
#define NBPT 200    // Number of plot points (analytic solution)

double z1(double t, double m);
double z2(double t, double m);

const double m = 2.5;
const double T = 6.0*PI;

void deriv(const DP t, Vec_I_DP &y, Vec_O_DP &f);

int main()
{
  int i, j;
  double hh, t0 = 0.0;
  DP t, h;
  Vec_DP y(N), yout(N), yerr(N), dydx(N);

  int fargc=-1;
  char **fargv;
  TApplication theApp("theApp",&fargc,fargv);

  double time;
  double tt[NBPT], ttnum[NP];
  double yy1[NP], yy2[NP], zz1[NBPT], zz2[NBPT];

  h = T/(double)NP;
  y[0] = y[1] = y[3] = 0.0;
  y[2] = 1.0/sqrt(m);
  t = t0;
  time = t0;
  cout << endl;
  cout << "# time" << "\t y1-exact" << "\t y1-num.   ";
            cout << "\t y2-exact" << "\t y2-num" << endl;
  cout << "# ";
  for (int i=0;i<70;i++) cout << "-";
  cout << endl;
  cout << time << setw(16) << z1(time,m) << setw(12) << y[0] << setw(18);
                      cout << z2(time,m) << setw(12) << y[1] << endl;
 
  for (int j=0;j<NP;j++) {
    ttnum[j] = t;
    yy1[j] = y[0];
    yy2[j] = y[1];  
    deriv(t, y, dydx);
    NR::rkck (y, dydx, t, h, yout, yerr, deriv);
    for (int i=0; i<N; i++)  y[i] = yout[i];
    t += h;
    time += h;
    
    cout << time << setw(16) << z1(time,m) << setw(12) << yout[0] << setw(16);
                        cout << z2(time,m) << setw(12) << yout[1] << endl;
  }

  // root graphics:
  tt[0] = t0;
  zz1[0] = z1(tt[0],m);
  zz2[0] = z2(tt[0],m);
  hh = T/(double)NBPT; 

  for (int j=1;j<NBPT;j++) {
    tt[j] = tt[j-1] + hh;
    zz1[j] = z1(tt[j],m);
    zz2[j] = z2(tt[j],m);
  }
  
  TGraph *grph1=new TGraph(NBPT,tt,zz1);
  grph1->SetLineColor(3);
  TGraph *grph2=new TGraph(NBPT,tt,zz2);
  grph2->SetLineColor(4);
  TGraph *grph3=new TGraph(NP,ttnum,yy1);
  grph3->SetMarkerColor(8);
  TGraph *grph4=new TGraph(NP,ttnum,yy2);
  grph4->SetLineColor(6);
  TCanvas *c1 = new TCanvas("c1","spring-mass system",800, 400);
  c1->SetGrid();
  c1->cd(); 
  grph1->Draw("AC");
  grph2->Draw();
  grph3->Draw("AC*");
  grph4->SetMarkerStyle(21);
  grph4->Draw("CP");
  grph3->SetTitle("Spring-mass system (with 2 masses)");
  c1 ->  Update();

  theApp.Run();
  return 1;
}



// =============================================
// RHS of the diff. equation

void deriv(const DP t, Vec_I_DP &y, Vec_O_DP &f)
{
  f[0] = y[2];
  f[1] = y[3];
  f[2] = -2.0*y[0]/m + y[1]/m;
  f[3] = y[0]/m - 2.0*y[1]/m;
}

// realization of the exact solution (test example)

double z1(double t, double m)
{
  return( sqrt(3.0)/6.0*sin(t*sqrt(3.0)/sqrt(m)) + sin(t/sqrt(m))/2.0 );
}

double z2(double t, double m)
{
  return( sin(t/sqrt(m))/2.0 - sqrt(3.0)/6.0*sin(t*sqrt(3.0)/sqrt(m)) );
}

// =====================================================================
// (vereinfachte) Gnuplot-Grafik:
// grxx054 > gr054.dat
// gnuplot:
// plot "gr054.dat" u 1:2 w l, "gr054.dat" u 1:3, \
//                     "gr054.dat" u 1:4 w l, "gr054.dat" u 1:5
//   2. und 4. Komponente: analyt. Loesung, 3. und 5. Komp: num. Loesung