//A rossler iterator, Tim Jones
#include <stdlib.h>
#include <math.h>
#include <fstream.h>
#include <iostream.h>
#include <strstream.h>
#include "nr.h"
#include "nrutil.h"
#include "nrutil.c"
#include "rk4.c"
// \dot{x} = -y - z
// \dot{y} = x + ay
// \dot{z} = b + z(x-c)
//We use RK4 to iterate this system and recreat an attractor for analysis

int main(){
 int i1,i2,i3,j,n1,n,N;  double x,a,b,c; double h=0.01; double *y; 
 double *dydx; double *yout;
//INPUT NEEDED FROM USER HERE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 n=3;
//THIS WILL BE THE NUMBER OF DIMENSIONS INVOLVED  !!!!!!!!!!!!

 y = dvector(1,n);  dydx = dvector(1,n); yout = dvector(1,n);

//NOW YOU WILL HAVE TO SET YOUR INITIAL CONDITIONS  !!!!!!!!!!
 y[1]=0.25; y[2]=0.0; y[3]=0.1;
//ONE COULD MAKE THIS A RANDOM PROCESS SO LONG AS x \in A !!!!

//this first main loop is here to remove original transients *

 for(i2=0; i2 <=10000; i2++){
   x = i2 * h;
   derivs(x,y,dydx);
   rk4(y,dydx,n,x,h,yout,derivs);
   for(j=1; j<=n; j++){y[j]=yout[j];}}
 //this ends the first main loop

 //this second loop is where the real operation is ***********

 x = 0.0;
 //INPUT NEEDED FROM USER HERE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ofstream outf("doodle.dat");
 //THIS WILL BE THE FILE YOU OUTPUT TO  !!!!!!!!!!!!!!!!!!!!!!!

 for(i2=0; i2<=1024; i2++){
  for(i3=0; i3<=9; i3++){//We iterate ten times and output only 1/10th
    x = x + h;
    derivs(x,y,dydx);
    rk4(y,dydx,n,x,h,yout,derivs);
    for(j=1;j<=n; j++){y[j]=yout[j];}}  //ends i3 loop
  outf<<i2<<" "<<y[1]<<" "<<y[2]<<" "<<y[3]<<"\n";
}  //ends i2 loop

double derivs(x,*y,*dydx){ a=0.398; b = 2.0; c=4.0;
 dydx[1] = -y[2]-y[3];
 dydx[2] = y[1] + a*y[2];
 dydx[3] = b + y[3]*(y[1]-c); }