// Loesung Blatt 2 AR-Prozesse und Van der Pol Oszillator
// 5.11.2002 Kilian Bartholomé

#include <time.h>
#include <iostream>
#include <fstream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "/home/maiwald/include/nr.h"
#include "/home/maiwald/include/nrutil.h"

#define PI 3.141592653589

#ifdef WIN32
    #include <sys/timeb.h>
#else
    #include <sys/time.h>
#endif

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

  //Zugriff auf Systemzeit

  long millisekunden = 0;
#ifdef WIN32
  struct __timeb64 timebuffer;
  _ftime64( &timebuffer );
  millisekunden = timebuffer.millitm;
#else 
  struct timeval tp;
  gettimeofday(&tp, NULL);
  millisekunden = tp.tv_usec;
#endif
  long idum = - millisekunden;


  //Berechnung AR[1] Prozeß
  ofstream AR1out("AR1.txt",ios::out);
  int N=8192;
  double tauAR1=5.0;
  double aAR1=exp(-1.0/tauAR1);
  double xAR1[N];
  double startwertAR1=0.1;
  xAR1[0]=startwertAR1;
  
  for(int i=1;i<N;i++)
    {
      xAR1[i]=aAR1*xAR1[i-1]+gasdev(&idum);
      AR1out<<xAR1[i]<<endl;
    }
  
  //Berechnung AR[2] Prozeß    
  ofstream AR2out("AR2.txt",ios::out);
  double tauAR2=20.0;
  double TAR2=20.0;
  double a1AR2=2.0*cos(2.0*PI/TAR2)*exp(-1.0/tauAR2);
  double a2AR2=-exp(-2.0/tauAR2);
  double xAR2[N];
  double startwert1AR2=0.1;
  double startwert2AR2=0.11;
  xAR2[0]=startwert1AR2;
  xAR2[1]=startwert2AR2;
  
  for(int i=2;i<N;i++)
    {
      xAR2[i]=a1AR2*xAR2[i-1]+a2AR2*xAR2[i-2]+gasdev(&idum);
      AR2out<<xAR2[i]<<endl;
    }
  
  
  //Berechnung van der Pol Oszillator
  N=50000;
  ofstream VDPout("VDP.txt",ios::out);
  double x1VDP[N];
  double x2VDP[N];
  double mu=1.0; //1.0 3.0 5.0
  double deltat=0.001;
  double sampling=0.1;
  double counter=0;
  double startwertVDP1=1;
  double startwertVDP2=1;
  x1VDP[0]=startwertVDP1;
  x2VDP[0]=startwertVDP2;
  
  for(int i=1;i<N;i++)
    {
      x1VDP[i]=x1VDP[i-1]+deltat*x2VDP[i-1];
      x2VDP[i]=x2VDP[i-1]+deltat*(mu*(1.0-x1VDP[i-1]*x1VDP[i-1])*x2VDP[i-1]-x1VDP[i-1])+sqrt(deltat)*gasdev(&idum);
      if(counter>=sampling)
	{
	  VDPout<<x1VDP[i]<<" "<<x2VDP[i]<<endl;
	  counter=0;
	}
      counter=counter+deltat;
    }
  
  
  return(0);
}