/*
find_truc.c
Copyright (C) 2009 Timothy Jones
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
//////////////////////////////////////
#include
#include
#include
#include
#include
#include
#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; i20)
{
if(flag1==true)
{ //If we have hit PS at least period N times
hit=-1; flag1=false;
///////////////////////////////
if(yout[2]ymin){
if(y[2]ymin){
///////////////////////////////
if(level==10){outf2<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<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]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 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