#include <stdio.h> 
#include <math.h>
#include"nrutil.h"
#include"nrutil.c"

#define PI2 6.2831853


double w;
double A; 


void takens(int n, double v[], double dv[], double t){
	double x,y;
	x=v[0]; y=v[1];;

 	dv[0] = x-x*x*x+y;
 	dv[1] = -0.1*(x-0.1) - A*cos(w*t);
}


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

  int n_it,n_pre,i,j,period,count,size;
  int imin, imax, jmin, jmax, total;
  double v[2],t,dt,t_pre,counter,told;
  double v2[2],v3[2],v4[2],v5[2],v0[2],v1[2], *series;
  FILE *ptr, *ptr2, *ptr3, *ptr4, *ptr5;

  ptr  = fopen("vanderpol.dat","wt");
  ptr2 = fopen("vanderpol1.dat","wt");
  ptr3 = fopen("vanderpol2.dat","wt");
  ptr4 = fopen("vanderpol3.dat","wt");
  ptr5 = fopen("vanderpol4.dat","wt");


  imin = 10;
  imax = 50;
  jmin = 10;
  jmax = 40;

  total = (imax-imin)*(jmax-jmin);


  for(i = imin; i < imax; i++){
	for(j = jmin; j < jmax; j++){

		A = i*3/400.0;
		w = j*2.28/400.0;

		period = 0;
		  told = 0;
		    dt = 0.01;  
		 t_pre = 100.;
		     t = 200.; 


		v[0] = 0.1; v[1] = 0.1;

		while(t > 0.){
		   rk4(takens,v,2,t,dt);
		     told = t;
		       t -= dt;
 		   if(t < t_pre){	 
			if((told > 1*PI2/w) && (t < 1*PI2/w)) {v1[0] = v[0]; v1[1] = v[1];};
			if((told > 2*PI2/w) && (t < 2*PI2/w)) {v2[0] = v[0]; v2[1] = v[1];};
			if((told > 3*PI2/w) && (t < 3*PI2/w)) {v3[0] = v[0]; v3[1] = v[1];};
			if((told > 4*PI2/w) && (t < 4*PI2/w)) {v4[0] = v[0]; v4[1] = v[1];};
			if((told > 5*PI2/w) && (t < 5*PI2/w)) {v5[0] = v[0]; v5[1] = v[1];};
		  }
		}

		if( sqrt((v5[0]-v1[0])*(v5[0]-v1[0]) + (v5[1]-v1[1])*(v5[1]-v1[1]))<0.01) period=4;
		if( sqrt((v4[0]-v1[0])*(v4[0]-v1[0]) + (v4[1]-v1[1])*(v4[1]-v1[1]))<0.01) period=3;
		if( sqrt((v3[0]-v1[0])*(v3[0]-v1[0]) + (v3[1]-v1[1])*(v3[1]-v1[1]))<0.01) period=2;
		if( sqrt((v2[0]-v1[0])*(v2[0]-v1[0]) + (v2[1]-v1[1])*(v2[1]-v1[1]))<0.01) period=1;

                
		if(!(total%1000)) {printf("%ld - ",total); fflush(stdout);} total--;

		if(period==1)  fprintf(ptr2,"%lg\t%lg\t%ld\n", w, A, period);
		if(period==2)  fprintf(ptr3,"%lg\t%lg\t%ld\n", w, A, period);
		if(period==3)  fprintf(ptr4,"%lg\t%lg\t%ld\n", w, A, period);
		if(period==4)  fprintf(ptr5,"%lg\t%lg\t%ld\n", w, A, period);

	}
}
fclose(ptr);
fclose(ptr2);
fclose(ptr3);
fclose(ptr4);
fclose(ptr5);

exit(0);
}