/* 
************************************************************************
*  カオス: Duffing方程式
************************************************************************
  t：時間、x：位置、v：速度、dt：時間の刻み幅
  t0、x0、v0：初期値
  f0：速度に比例する抵抗の係数
  tlim：時間の最大値

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

#define N 2                                   /* number of equations */

double a, b, f0, omega;
 
main()
{
   void rungen(double t, double x[], double dt);
   double f(double t, double x[], int i);
   
   double t, x[N], v0, dt, h;
   double pi,tlim;
   int i, j;
   
/* 出力ファイルの指定 */
   FILE *output;
   FILE *output2;
   output = fopen("duffing.dat","w");
   output2 = fopen("poincare.dat","w");

   pi=3.1415926535897932385;
   a=0.05;
   b=1.0;
   omega=1.0;
   dt=2*pi/1000;

/*  強制力の振幅 */
   printf("Input: f0\n");
   scanf("%lf", &f0);
   printf("f0=%g\n",f0);
   printf("Input: tlim/(2*pi)\n");
   scanf("%lf", &tlim);
   printf("tlim=%g\n",tlim);
   tlim=2*pi*tlim;

/* 初期値 */
   i=0;
   t = 0;
   v0 = 1.0;
   x[0] = 0.0;                                /* initial position  */
   x[1] = v0;                                /* initial velocity  */
   
   fprintf(output, "%f\t%f\t%f\n", t, x[0], x[1]); 
   fprintf(output2, "%f\t%f\n", x[0], x[1]);


   while(t <= tlim)
   {
     i++;
/* ルンゲ・クッタ法ルーチンを呼び出す */
      rungen(t, x, dt);
      t += dt;    
      fprintf(output, "%f\t%f\t%f\n", t, x[0], x[1]); 
      if(i == (i/1000)*1000)  fprintf(output2, "%f\t%f\n", x[0], x[1]);
   }
   printf("data stored in duffing.dat\n");
   fclose(output);
}
/*-----------------------end of main program--------------------------*/

/* Runge-Kutta subroutine */
void rungen(double t, double x[], double step)
{
   double f(double t, double x[], int i);
   double h=step/2.0,                         /* the midpoint */
          t1[N], t2[N], t3[N],                /* temporary storage */
          k1[N], k2[N], k3[N],k4[N];          /* for Runge-Kutta  */
   int i;
 
   for (i=0; i<N; i++) t1[i] = x[i]+0.5*(k1[i]=step*f(t, x, i));
   for (i=0; i<N; i++) t2[i] = x[i]+0.5*(k2[i]=step*f(t+h, t1, i));
   for (i=0; i<N; i++) t3[i] = x[i]+    (k3[i]=step*f(t+h, t2, i));
   for (i=0; i<N; i++) k4[i] =                 step*f(t + step, t3, i);

   for (i=0; i<N; i++) x[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
/*--------------------------------------------------------------------*/

/* definition of equations - this is the harmonic oscillator */
double  f(double t, double x[], int i)
{
   if (i == 0) return(x[1]);               /* RHS of first equation */
   if (i == 1) return(-a*x[1]-b*pow(x[0],3)+f0*cos(omega*t));   /* RHS of second equation */
}