//A rossler iterator
#include <stdlib.h>
#include <math.h>
#include <fstream.h>
#include <iostream.h>
#include <stream.h>
//#include <strstream.h>
#include "nr.h"
#include "nrutil.h"
#include "nrutil.c"
#include "rk4.c"

int main(int argc, char **argv){
  int i,j,n; double x,a,b,c; double h=0.01; double *y, *dydx, *yout;
  n=3;  //dimensions
  y = dvector(1,n);  dydx = dvector(1,n); yout = dvector(1,n);
  y[1]=0.25; y[2]=0.0; y[3]=0.1;  //initial starting seed
  ofstream outf("data.dat"); 
  ofstream outf2("report.dat");
  //  a=0.2;
  //for the boss
  //a=0.432;
   a=atof(argv[1]);
 double xx=0;  double xy=0; double xz=0;
 double yy=0;  double yz=0; double zz=0;
 double RR=0;  double dxx=0;  double dyy=0; 
 double dzz=0; double Lx=0; double Ly=0; double Lz=0;
 // printf("%f \n",a);
  for(i=0; i<=10000; i++){
      x = i * h; 
      derivs(a,x,y,dydx);
      rk4(a,y,dydx,n,x,h,yout,derivs); 
      for(j=1; j<=n; j++){y[j]=yout[j];}}

  //  while(a<=0.432){
 xx=0;  xy=0; xz=0;
 yy=0;  yz=0; zz=0;
 RR=0;  dxx=0; dyy=0; 
 dzz=0; Lx=0;  Ly=0; Lz=0;


   x=0;
   for(i=1; i<=50000; i++){ 
     x = x + h;
    derivs(a,x,y,dydx);
    rk4(a,y,dydx,n,x,h,yout,derivs);
    for(j=1;j<=n; j++){y[j]=yout[j];}
    outf<<" "<<y[1]<<" "<<y[2]<<" "<<y[3]<<"\n";
    xx+=y[1]*y[1];
    xy+=y[1]*y[2];
    xz+=y[1]*y[3];
    yy+=y[2]*y[2];
    yz+=y[2]*y[3];
    zz+=y[3]*y[3];
    RR+=xx+yy+zz;
    dxx+=dydx[1]*dydx[1];
    dyy+=dydx[2]*dydx[2];
    dzz+=dydx[3]*dydx[3];
    Lx+=y[2]*dydx[3]-y[3]*dydx[2];
    Ly+=y[3]*dydx[1]-y[1]*dydx[3];
    Lz+=y[1]*dydx[2]-y[2]*dydx[1];//}
    outf2<<" "<<x<<" "<<(dxx+dyy+dzz)/(2*h*50000)<<" "<<"\n";

    //outf2<<" "<<a<<" "<<xx/50000<<" "<<xy/50000<<" "<<xz/50000<<" "<<yy/50000<<" "<<yz/50000<<" "<<zz/50000<<" "<<RR/50000<<" "<<(dxx+dyy+dzz)/(2*h*h*50000)<<" "<<Lx/(h*50000)<<" "<<Ly/(h*50000)<<" "<<Lz/(h*50000)<<"\n";
   }
   // a=a+0.001;}  //ends a while loop
}  //ends main

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