//////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fstream>
#include <iostream>
#include <strstream>
#include "nr.h"
#include "nrutil.h"
#include "nrutil.c"
#define NRANSI
using namespace std;

#define period 1


void derivsy(double, double[], double[],double, double, double);

///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////


void derivsy(double x, double y[], double dydx[],double a, double b, double c){ 

dydx[1] = -y[2]-y[3]; 
dydx[2] = y[1] + a*y[2]; 
dydx[3] = b + y[3]*(y[1]-c);}

///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////

void mrk4(double y[], double dydx[], int n, double x, double h, double yout[],
	 void (*derivs)(double, double [], double [], double, double, double),double a, double b, double c)
{

  REMOVED DUE TO NUMERICAL RECIPES COPYRIGHT ISSUES, please see 
  http://www.gnu.org/software/gsl/ for alternative version of 
  Runge Kutta.

}        

///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////

void poin(double ymin, double ymax, double init[], double n,ofstream& outf2, double **close,int NUMLOOP, bool flaginit, int level, double a, double b, double c)
{

double h=0.01; double *y; double *dydx; double *yout;
double x; double yold,yo,lclose,lcloseold; 
double *trail;
bool flag1 = false;  //flag1:  have we had a first hit on the Poincare section?
bool flagclose = true;
y = dvector(1,n);  dydx = dvector(1,n); yout = dvector(1,n);
trail=dvector(1,8);


for(int i1=1; i1<=7; i1++){ trail[i1]=0.0;}

int closecount=1;

//Attractor specific data
y[1]=init[1]; y[2]=init[2]; y[3]=init[3];

//b=2.0; c=4.0;
//a=0.555999999999983;

double fixedx = c/2-sqrt(c*c -4 * a *b)/2;
double fixedy = -c/(2*a) + sqrt(c*c -4 *a*b)/(2*a);

//Run through first loop to get into attractor

if(flaginit==true)
{

  for(int i2=0; i2<10000; i2++)
    {
      x += h;
      derivsy(x,y,dydx,a,b,c);
      mrk4(y,dydx,n,x,h,yout,derivsy,a,b,c);
      for(int j=1; j<=n; j++){y[j]=yout[j];}
    }

flaginit=false; 
}



x=0;
yold=0.0;  
lclose=0.0;
lcloseold=0.0;
//Real loop is here
int hit=0;

double numdiv;

if(level<10){numdiv =1;}
else{numdiv=4;}

for(int i2=0; i2<NUMLOOP; i2++)
{

 x += h;
 yo = y[1];
 derivsy(x,y,dydx,a,b,c);
 mrk4(y,dydx,n,x,h,yout,derivsy,a,b,c);

 // If yo*yout[1] < 0, then it has just changed signs

if(((yo-fixedx)*(yout[1]-fixedx) < 0) && (yout[2]< fixedy) && dydx[1]>0)
{
       if(flag1==true)
       {  //If we have hit PS at least period N times
	 hit=-1; flag1=false;
       ///////////////////////////////
	 if(yout[2]<ymax && yout[2]>ymin){
	   if(y[2]<ymax && y[2]>ymin){
       ///////////////////////////////
            if(level==10){outf2<<trail[8-period]<<" "<<yout[2]<<" "<<fabs(yout[2]-yold)<<"\n";}
	     // outf2<<trail[8-period]<<" "<<yout[2]<<" "<<fabs(trail[8-period]-yold)<<"\n";

	     if(fabs(yout[2]-trail[8-period])<0.2/(1.0*numdiv*period) )
	     {   //Now zoom in on candidates 
     
	     //          lclose=yold;    
	        lclose=trail[8-period];
		if(level==1)
		{
		    for(int i=1; i<closecount; i++){if(fabs(lclose-close[i][2])< 0.2/(1.0*period)){flagclose=false;}}
		    if(flagclose==true){
		      for(int j=1; j<=n; j++){close[closecount][j]=yout[j];}
		      closecount++;
		    } 
                flagclose=true;
	        }
          //For higher iterations, collect all around each ~intersection
		if(level==10)
		{ 
		  for(int j=1; j<=n; j++){close[closecount][j]=yout[j];} 
		  closecount++;
		}
             }
       ////////////////////////
	  }}
       ////////////////////////
       }
      
      for(int b1=1; b1<7; b1++){ trail[b1]=trail[b1+1];}
      trail[7]=yout[2];
      yold=yout[2];  hit++;
      //  printf("hit: %i, period: %i \n",hit,period);
      if(hit==period){flag1=true;}//After our first hit, we can start recording

     }
 
 //Reassign for next loop
for(int j=1; j<=n; j++){y[j]=yout[j];}
//outf<<y[1]<<" "<<y[2]<<" "<<y[3]<<"\n";
}

flaginit=true;
free_dvector(y,1,n);
free_dvector(dydx,1,n);
free_dvector(yout,1,n);
free_dvector(trail,1,n);
}


///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////



void regular(double init[], int n, double output[], ofstream& outf, bool printflag, double a, double b, double c){


  double h=0.01; double *y; double *dydx; double *yout; double *trail;
double x; double yold,yo; 

 int safety=1;

bool flag1 = false;  //flag1:  have we had a first hit on the Poincare section?
bool flagclose = true;
bool flagreg;

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

int closecount=1;

//Attractor specific data
y[1]=init[1]; y[2]=init[2]; y[3]=init[3];



//b=2.0; c=4.0;
//a=0.555999999999983;

double fixedx = c/2-sqrt(c*c -4 * a *b)/2;
double fixedy = -c/(2*a) + sqrt(c*c -4 *a*b)/(2*a);

flagreg=false;
int hit=0;
while(flagreg==false){
 
 x += h;
 yo = y[1];
 derivsy(x,y,dydx,a,b,c);
 mrk4(y,dydx,n,x,h,yout,derivsy,a,b,c);
 for(int i=1; i<=n; i++){y[i]=yout[i];}
 //Which side of the poincare section is it starting on?  Shouldn't matter.
 //Safety says must iterate certain amount before counting crossing.

 if(safety>100){if(((yo-fixedx)*(yout[1]-fixedx) < 0) && (yout[2]< fixedy)  && dydx[1]>0 ){hit++; if(hit==period){for(int i=1; i<=n; i++){output[i]=y[i];}flagreg=true;}}}

 if(y[1]*y[1] + y[2]*y[2] + y[3]*y[3] > 1000){flagreg=true; printf("ERROR:  Orbit strayed too far.  Exiting now. \n"); exit(1);}
 safety++;

 if(printflag==true){outf<<y[1]<<" "<<y[2]<<" "<<y[3]<<"\n";}

 if(safety>100000){flagreg=true;  printf("ERROR:  Orbit iteration exceeded maximum.  Exiting now. \n"); exit(1);}


 
 // if(flagtrue==true){for(int i=1; i<=n; i++){output[i]=yout[i];}}

 }


free_dvector(y,1,n);
free_dvector(dydx,1,n);
free_dvector(yout,1,n);

}


///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////

int main(){
   

  double a,b,c;
  
  b=2.0; c=4.0;
  a=0.555999999999983;
  
 
  system("clear");
  printf("\n");
  printf("////////////////////  Rossler Periodic Orbit Finder  ////////////////////\n");
  printf("\n");
  printf("Running global search for candiate matches (return map intersection)...\n");

  
  ofstream outf2("out2.dat");
  int n;
  double diff, mindiff,smallx,smally,smallz;
  double *init;  
  double **close;
  double **closed;
  int counter=0;
  double ymax=-c/(2*a) + sqrt(c*c -4 *a*b)/(2*a); //fixedy
  double ymin=-12.0;   //covers entire attractor
  int bignum=2000;
  int NUMLOOP=10000000 + 5000000*period;
  bool flaginit=true;
  bool printflag=true;
  double level=1;
  double dx,dy,dz;
  double *output;
  int mesh=20;
 

  n=3;
  init   = dvector(1,n);  
  output = dvector(1,n);
  
  
  bool flagfin=false;
  double xgmin,xgmax,ygmin,ygmax,zgmin,zgmax,yorg;
  
  close  = dmatrix(1,bignum,1,n);
  for(int i=1; i<=bignum; i++){for(int j=1; j<=n; j++){close[i][j]=0.0; output[j]=0.0;}}

  init[1]=0.25;
  init[2]=0.0;
  init[3]=0.1;
   
  printf("Starting (ymin,ymax):  %f %f \n",ymin,ymax);
 
  poin(ymin,ymax,init,n,outf2,close,NUMLOOP,flaginit,level,a,b,c);

  for(int i=1; i<=bignum; i++){if(close[i][2] != 0.0){counter++; 
                  printf("Candidate (x,y,z):  %f %f %f \n",close[i][1],close[i][2],close[i][3]);}}
  printf("===========>  %i candidates found.\n", counter);

  //Now scan

  //Clone the close list because we will need to reuse close
  closed  = dmatrix(1,counter,1,n);
  for(int i=1; i<=counter; i++){for(int j=1; j<=n; j++){closed[i][j]=close[i][j];}}
  
  ////////////////////////////////////////////////////////////////////////////////////////////////////

  int counted=counter;
  const char cc[] = "data00";
  char filename [ 10 ];
  for(int k=1; k<=counted; k++){
 
  printf("______________Candidate # %i ______________ \n", k);  
  printf("\n");
  
  for(int m=1; m<=n; m++){init[m]=closed[k][m];}
  
  //  (init[2]-0.02, init[2]+0.2)
  
  //Consider saving data, we end 
  init[1]=0.25;
  init[2]=0.0;
  init[3]=0.1;
 

  for(int i=1; i<=bignum; i++){for(int j=1; j<=n; j++){close[i][j]=0.0;}}

     level=10;
     NUMLOOP+=level*900000;
     counter=0;
     ymin = closed[k][2]-2/(1.0*level*period);
     ymax = closed[k][2]+2/(1.0*level*period);

     printf("Searching candidate with y value of %f.\n",init[2]);
     printf("Restarting with (ymin,ymax):  %f %f \n", ymin,ymax);
     poin(ymin,ymax,init,n,outf2,close,NUMLOOP,flaginit,level,a,b,c);      

     xgmin=close[1][1]; ygmin=close[1][2]; zgmin=close[1][3];
     xgmax=close[1][1]; ygmax=close[1][2]; zgmax=close[1][3];
     for(int ii=1; ii<=bignum; ii++){
                    if(close[ii][2] != 0.0)
		      {
                      counter++; 
		      if(close[ii][1]<xgmin){xgmin=close[ii][1];}
		      if(close[ii][2]<ygmin){ygmin=close[ii][2];}
		      if(close[ii][3]<zgmin){zgmin=close[ii][3];}
		      if(close[ii][1]>xgmax){xgmax=close[ii][1];}
		      if(close[ii][2]>ygmax){ygmax=close[ii][2];}
		      if(close[ii][3]>zgmax){zgmax=close[ii][3];}
		      printf("Close points (x,y,z):  %f %f %f \n",close[ii][1],close[ii][2],close[ii][3]);
                    }
     }
     
 printf("===========>  %i close points found. \n\n", counter);
 printf("Search grid obtained: \n");
 printf(" xmin: %f, xmax: %f \n ymin: %f, ymax: %f \n zmin: %f, zmax: %f \n\n",xgmin,xgmax,ygmin,ygmax,zgmin,zgmax);
        
 printf(".......Now search grid for best fit periodic orbit.......\n");
 printf("_________________________________________________________\n");

 if(xgmin*xgmin > 0.000){

 sprintf(filename, "%s%d", cc,k);
 ofstream outf(filename);
 dx=(xgmax-xgmin)/(1.0*mesh); dy=(ygmax-ygmin)/(1.0*mesh); dz=(zgmax-zgmin)/(1.0*mesh);
 diff=1000;
 mindiff=1000;
 printflag=false;

 for(int k1=0; k1<=mesh; k1++){
   for(int k2=0; k2<=mesh; k2++){
     for(int k3=0; k3<=mesh; k3++){
       init[1]=xgmin+dx*k1;
       init[2]=ygmin+dy*k2;
       init[3]=zgmin+dz*k3;
       regular(init,n,output,outf,printflag,a,b,c);
       //       printflag=false;
       diff=((xgmin+dx*k1)-output[1])*((xgmin+dx*k1)-output[1]) + ((ygmin+dy*k2)-output[2])*((ygmin+dy*k2)-output[2]) +((zgmin+dz*k3)-output[3])*((zgmin+dz*k3)-output[3]);
       //       diff=fabs((ygmin+dy*k2)-output[2]);
       if(diff<mindiff){mindiff=diff; smallx=output[1]; smally=output[2]; smallz=output[3]; printf("---------: output:  %f %f %f %f ------------: \n",output[1],output[2],output[3],diff);
}
       //       printf("Initial (x,y,z): %f %f %f ----------------> Final: %f %f %f %f \n",xgmin+dx*k1,ygmin+dy*k2,zgmin+dz*k3,output[1],output[2],output[3],diff);
 }}}

 printf("\n");
 printf("==============> Best fit orbit has initial value at (x,y,z): %f %f %f with geometric distance %f\n", smallx,smally,smallz,mindiff);
 
 if(mindiff<0.5){
 printf("Fit is within standard.\n");
 init[1]=smallx; init[2]=smally; init[3]=smallz;

 printf("...Now printing data for period [%i] orbit...\n\n",period);
 
 printflag=true;
 regular(init,n,output,outf,printflag,a,b,c);
 outf.close();
 }


 }
 
 }

 printf("Program exiting normally.\n");
 printf("\n");
          
          



}

#undef NRANSI