#ifndef _CINT_
#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

Vec_DP *xp_p;
Mat_DP *yp_p;
Mat_DP *fp_p;

const int n_ord = 3;
const double a = 8.0;   // global parameter
const double ya = 2.0;  // initial condition

double yex(double t)
{
  // analytical solution of the ODE 
  double a2;
  a2 = 1.0/(1.0+a*a);
  return( (ya - a2)*exp(-t) + a2*( cos(a*t) + a*sin(a*t) ) );
}

void deriv(const DP t, Vec_I_DP &y, Vec_O_DP &f)
  // RHS of ODE
{
  f[0] =  cos(a*t) - y[0];
}

void adamsbash(int k, Vec_IO_DP &y, const DP h, Vec_IO_DP xp,
           Mat_IO_DP yp, Mat_IO_DP fp, const int n_ord, Vec_O_DP &yout )
{
  const int NMAXORD = 5;
  static DP beta[NMAXORD][NMAXORD] = {{1.0, 0.0, 0.0, 0.0, 0.0},
              {(3.0/2.0), (-1.0/2.0), 0.0, 0.0, 0.0},
              {(23.0/12.0),(-16.0/12.0),(5.0/12.0), 0.0, 0.0},
              {(55.0/24.0),(-59.0/24.0),(37.0/24.0),(-9.0/24.0), 0.0},
              {(1901.0/720.0),(-2774.0/720.0),(2616.0/720.0),(-1274.0/720.0),
                                                            (251.0/720.0)}};
  int i,j;

  int nvar=y.size();
  Vec_DP sum(nvar);

  if (n_ord > NMAXORD || n_ord < 1) 
                  NR::nrerror("Formula not present");

  for (i=0; i<nvar; i++) sum[i] = 0.0;

  for (j=0; j < n_ord; j++) {
      for (i=0; i<nvar; i++) 
           sum[i] += beta[n_ord-1][j]*fp[i][k-j-1];
  }

  for (i=0; i<nvar; i++) {
          yout[i] = y[i] + h*sum[i];
          yp[i][j+1] = yout[i];
  }

}

int main()
{
   const int N=1;       // dimension of system
   int i, j, NN =100;   // ## plot points
   DP hh, tt;
   Vec_DP y(N), yout(N), yy(N), f(N);
   double y1[NN], y2[NN], y3[NN], t1[NN];

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

   hh = 2.0*PI/(double)NN;
   tt = 0.0;
   y[0] = ya; 

   xp_p = new Vec_DP(NN);
   yp_p = new Mat_DP(N,NN);
   fp_p = new Mat_DP(N,NN);
   Vec_DP &xp = *xp_p;
   Mat_DP &yp = *yp_p;
   Mat_DP &fp = *fp_p;

   // starting procedure
   xp[0] = tt; 
   deriv(tt,y,f);
   for (i=0;i<N;i++) {
      yp[i][0] = y[i];
      fp[i][0] = f[i];
   }
   cout << "# Start: " << "n_ord = " << n_ord << endl;
   //  cout << "j  = 0 " << "\t" << xp[0] << "\t" << yp[0][0];
   //                                 cout << "\t" << fp[0][0] << endl; 
 
   for (j=1; j < n_ord; j++) {
      deriv(tt,y,f);
      NR::rk4(y, f, tt, hh, y, deriv);
      tt += hh;
      xp[j] = tt;
      deriv(tt,y,f); // function value of solution point
      for (i=0;i<N;i++) {
          yp[i][j]=y[i];
          fp[i][j]=f[i];
      }
  
   }
   // cout << " -> end of starting procedure ..." << endl;

   // multistep method ( Adams-Bashforth, order <n_ord> )
   for (j=n_ord; j<NN; j++) {
      tt = xp[j-1];
      for (i=0;i<N;i++) {
         y[i] = yp[i][j-1];
         f[i] = fp[i][j-1];
      }
      // multistep method: Adams-Bashforth, order 1-5
      adamsbash(j, y, hh, xp, yp, fp, n_ord, yout);
                        
      tt += hh;
      xp[j] = tt;
      deriv(tt,yout,f);
      for (i=0;i<N;i++) {
          yp[i][j]=yout[i];
          fp[i][j]=f[i];
      }
      // cout << "j = " << j << "\t" << xp[j] << "\t"
      //                           << yp[0][j] << "\t" << fp[0][j] << endl;
   }

// end of the main procedure
// Simple output for GNUPLOT:

   cout << "#-> analytical solution: " << endl;
   cout << "#  j |   t     |   yex(t)      |  f(t,yex(t))  " << endl;
   cout << "#======================================================= " << endl;
   tt = 0.0;    
   for (int j=0;j<NN;j++) {
     yy[0] = yex(tt);
     deriv(tt, yy, f);
     cout << j << setw(12) <<  tt << setw(14) << yy[0] << setw(14);
                                                  cout << f[0] << endl;
     tt += hh;
   }

   cout << endl << "#-> numerical solution;   n_ord = " << n_ord << endl;
   cout << "#  j |   t     |   yp(t)      |  f(t,yp(t))      " << endl;
   cout << "# ======================================================= " << endl;
   for (int j=0;j<NN;j++) {
     cout << j << setw(12) << xp[j] << setw(14) << yp[0][j] << setw(14);
                                                  cout << fp[0][j] << endl;
   }

   cout << endl << "#-> error:   n_ord = " << n_ord << endl;
   cout << "#  j |   t    |   yp(t)     |  yex(t)    | yp(t) -  yex(t)";
   cout << endl;
   cout << "# =======================================================" << endl;
   for (int j=0;j<NN;j++) {
      yy[0] = yex(xp[j]);
      cout << j << setw(12) << xp[j] << setw(14) << yp[0][j] << setw(14);
                 cout << yy[0] << setw(15) << abs(yp[0][j]-yy[0]) << endl;
   }

// ROOT Graphics
   for (int j=0; j<NN; j++) {
     t1[j] = xp[j];
     y1[j] = yex(t1[j]);
     y2[j] = yp[0][j]; 
     y3[j] = fabs(y1[j]-y2[j]);
   }

   TGraph *grph1=new TGraph(NN,t1,y1);
   grph1->SetLineColor(2);
   TGraph *grph2=new TGraph(NN,t1,y2);
   grph2->SetLineColor(3);
   TGraph *grph3=new TGraph(NN,t1,y3);
   grph3->SetLineColor(2);

   TCanvas *c1 = new TCanvas("c1","Exercise 7.3",800, 400);
   c1->SetGrid();
   c1->cd();
   grph2->Draw("");
   grph1->Draw("AC");
   grph2->SetMarkerStyle(22);
   grph2->Draw("CP");
   grph1->SetTitle("Ex. 7.3: y'(t) = cos(a*t) - y, y(0)=2");
   c1->Update();

   TCanvas *c2 = new TCanvas("c2","Exercise 7.3",800, 400);
   c2->SetGrid();
   c2->cd();
   grph3->Draw("AC");
   grph3->SetTitle("Ex. 7.3: error [abs(y(t) - y_k)]");
   c2->Update(); 

   theApp.Run();

   delete xp_p, yp_p, fp_p;
   return 1;
}