#include <stdio.h> 
#include <math.h>

#define PI2 6.2831853
#define N_NAME1_AUX 30

struct Par {
	double gamma;
	double k22; double k21; double k11;
	double w1; double w2; double delta1; double delta2;
	double A;} a;


void rk4(void deri(int , double [], double [], double ), \
double h[], int n, double t, double dt)
{
#define naux 260 

int i;
double k1[naux],k2[naux],k3[naux],k4[naux],h0[naux];
double dt2, dt6;

dt2=dt/2.;
dt6=dt/6.;

for (i = 0 ; i<n; i++)
	h0[i] = h[i];

deri(n,h0,k1,t);
for (i =0 ; i<n ; i++)
	h0[i]=h[i]+dt2*k1[i];

deri(n,h0,k2,t+dt2);
for (i =0 ; i<n ; i++)
	h0[i]=h[i]+dt2*k2[i];

deri(n,h0,k3,t+dt2);
for (i =0 ; i<n ; i++)
	h0[i]=h[i]+dt*k3[i];

deri(n,h0,k4,t+dt);
for (i = 0 ; i<n ; i++)
	{h0[i]=h[i]+dt*k4[i];};

for (i =0; i<n ; i++)
	h[i]=h[i]+dt6*(2.*(k2[i]+k3[i])+k1[i]+k4[i]);

return;
}




void takens(int n, double v[], double dv[], double t)
{
  double x,y;
  
                                                                                                 
  x=v[0]; y=v[1]; 
                                                                                                 
  dv[0]  = (a.w2*y-a.delta2*x)+(a.gamma/2.0)*(1-x*x+y*y)+(a.k22/2.0)*x*(1-x*x-y*y)+(a.k21/2.)*sqrt(1-(2*a.delta1/a.k11))*(cos(a.w1*t)-cos(a.w1*t)*(x*x-y*y)-2*x*y*sin(a.w1*t)); 
  dv[1]  = -(a.w2*x+a.delta2*y)+(a.gamma/2.0)*(-2*x*y)+  (a.k22/2.0)*y*(1-x*x-y*y)+(a.k21/2.)*sqrt(1-(2*a.delta1/a.k11))*(sin(a.w1*t)+sin(a.w1*t)*(x*x-y*y)-2*x*y*cos(a.w1*t)); 

                                                                               
return;
}




void nombrar1_par1(char *name_space1,int par1){
  sprintf(name_space1,"mapa.%d.dat",par1);
return;
}

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

  int ik,i,counter;
  double v[2],t,dt;
  double auxilia0[10000], auxilia1[10000];
  FILE *ptr1;
        char nombre_archivo1[N_NAME1_AUX];



for(ik=1;ik<=120;ik++){

nombrar1_par1(nombre_archivo1,ik);
ptr1=fopen(nombre_archivo1,"w");

 a.w1 = 4.5148;
 //a.k21= 4.44;
 a.k21= 4.75;
     a.w2 = 5.8;
    a.gamma=1.83+6*4/10.;
    a.k22=3.+7*4/10.;
    a.k11=5.;
	a.delta1=0.01;
	a.delta2=1.;
        dt = 0.01;  
      t=0; counter=0;

for(i=0;i<10000;i++){auxilia0[i]=-0.5+0.8*cos((PI2/10000.)*i);auxilia1[i]=0.8*sin((PI2/10000.)*i);}



 		 for(i = 0; i < 10000; i++){
		 v[0]=auxilia0[i]; v[1]=auxilia1[i];
		 t=0;
    while(t < 0.1*ik*PI2/a.w1){
                 rk4(takens, v, 2, t, dt);
	         t = t + dt;
		 }
                 fprintf(ptr1,"%lg\t%lg\t%lg\t%lg\t%lg\n",auxilia0[i],auxilia1[i],v[0],v[1],t);
 		   }
//	fprintf(ptr1,"\n");
//	fprintf(ptr1,"%d\n",counter);
//	fprintf(ptr1,"\n");

        counter++;
   
 
  fclose(ptr1);
}
}