#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'
#include <TApplication.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TF1.h>
#endif

#include <iostream>
#include <string>
#include <cmath>
using namespace std;

void fatal(char *text);
double fval(double time);
double fdotval(double time);

#define PI 3.1415926
#define NBPT 500    // Number of plot points

double sth = 1.0/double(17.0);
double sq2 = 1.0/sqrt(2.0);
double sq  = sqrt(2.0);

int main(int nbargs, char* args[])
{
  double  m=2.0;  // mass 
  double  c=3.0;  // spring rate 
  double  k=4.0;  // damping (decay)
  double  T=20.0; // length of the interval
  double  x0=1.0; // initial condition
  double  v0=0.0; // Starting velocity

  // if (nbargs!=2) fatal("Provide the starting path");
  // double start=atof(args[1]);

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

  double nextxdot;
  double nextx;
  double time=0;
  double t[NBPT];
  double x[NBPT];
  double xdot[NBPT];
  double h = T/double(NBPT);

  x[0] = x0;               // Starting value for the excursion 
  xdot[0] = v0;
  // Exact solution:
  for (int i=0;i<NBPT;i++){
    t[i]=time;
    x[i] = fval(time);           // exact solution
    xdot[i] = fdotval(time);     // derivative of the exact solution
    time=time+h;
    // Output: standard terminal, for plotting by 'gnuplot' e.g. 
    cout << '\t' << t[i] << '\t' << x[i] << '\t' << xdot[i] << '\n';
  }

  TGraph *grph1=new TGraph(NBPT,t,x);
  grph1->SetLineColor(2);
  TGraph *grph2=new TGraph(NBPT,t,xdot);
  grph2->SetLineColor(3);
  TGraph *grph3=new TGraph(NBPT,x,xdot);
  grph3->SetLineColor(2);

  TCanvas *c1 = new TCanvas("c1","Damped Oscillator",800, 400);
  c1->cd();
  grph2->Draw("");
  grph1->Draw("AC");
  grph2->Draw();
  c1->Update();

  TCanvas *c2 = new TCanvas("c2","Phase portrait of the Oscillator",400, 400);
  c2->cd();
  grph3->Draw("AC");
  c2->Update();

  theApp.Run();
  return 1;
}

//============================================================

double fval(double time)
// exact solution of the damped, inhomogeneous oscillator
{
  return exp(-time)*(12.0*sth*sq*sin(sq2*time)+16.0*sth*cos(sq2*time))
               + sth*(4.0*sin(time) + cos(time));
}

//===========================================================

double fdotval(double time)
// derivative of the solution
{
  return -exp(-time)*(20.0*sth*sq*sin(sq2*time)+4.0*sth*cos(sq2*time))
               + sth*(4.0*cos(time) - sin(time));
}

//============================================================

void fatal(char *text)
{
  cerr << text << "\n";
  cerr << "Program terminates now\n";
  exit(-1);
}
//============================================================