// Programm simuliert den Van der Pol - Oszillator
// unter Verwendung der Integrationsmethode Runge Kutta 
// 4. Ordnung (rk4 der Numerical Recipes).
// Thomas Maiwald, 10.11.2003

#include <iostream.h>
#include <stdlib.h>
#include <math.h>
#include "/home/maiwald/include/nrutil.h"
#include "/home/maiwald/include/nr.h"


// Globale Variable, damit auch die Funktion derivs
// auf den Wert von mu zugreifen kann.
double my;

// Funktion derivs wird von der Recipes-Methode rk4
// benoetigt. Hier nur die Signatur, Implementierung ist
// unterhalb der main-Routine.
void derivs(double, double *, double *);

//----------------------------------------------------------------------------

int main(int argc, char **argv)
{

  // Pruefe Anzahl der Uebergabeparameter
  if(argc != 6)
    {
      cerr << " Aufruf: ./vanderpol my y1 y2 h N " << endl;
      exit(1);
    }
  
  // Uebergabeparameter uebergeben
  double* y    = dvector(1,2); // Recipes-Vektor mit Komponenten y[1] und y[2]
  double* dydx = dvector(1,2);
  my           = atof(argv[1]);
  y[1]         = atof(argv[2]);
  y[2]         = atof(argv[3]);
  double h     = atof(argv[4]);
  int N        = atoi(argv[5]);

  // x ist hier die Zeit. Wir starten zum Zeitpunkt 0.
  double x=0;

  // Die Funktion mit den Ableitungen muss einmal aufgerufen werden:
  derivs(x,y,dydx);

  // Kern des Programms: N Integrationsschritte werden mit Runge-Kutta 
  // berechnet. Siehe Funktion rk4(...) der Numerical-Recipes.
  for(int i=1; i<=N; i++)
{
  // Ausgabe der Zeit und des Wertes
  cout<<x<<"\t"<<y[1]<<endl;

  // Zeit wird um h vergroessert
  x=x+h;

  // Runge-Kutta wird aufgerufen
  rk4(y,dydx,2,x,h,y,derivs);
}

  // Aufraeumen
  free_dvector(y,1,2);
  free_dvector(dydx,1,2);

  return 0;
} 

//-----------------------------------------------------------------------------

// Hier stehen die Ableitungen des Van-der-Pol-Oszillators:
void derivs(double x, double *y, double *dydx){

  dydx[1]=y[2];
  dydx[2]=my*(1-y[1]*y[1])*y[2]-y[1];

}