/*************************    Program 5.A    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define g1(y1,y2,t) (y2)
#define g2(y1,y2,t) (-q*y2-sin(y1)+b*cos(w*t))
#define NMAX 65536
double q,b,w;

main()
/* Program for the power spectra of a driven pendulum under
   damping with the fourth order Runge-Kutta algorithm. Given
   parameters: Q, B, and W (omega_0).  Copyright (c) Tao Pang 1997. */
{
int i,n,m,l,np;
double pi,rn,f1,h,od,t,y1,y2;
double dk11,dk21,dk12,dk22,dk13,dk23,dk14,dk24;
double y[2][NMAX],ar[NMAX],ai[NMAX],wr[NMAX],wi[NMAX],o[NMAX];
extern double q,b,w;
void fft();

n = NMAX;
l = 128;
m = 16;
pi = 4*atan(1);
rn = n;
f1 = 1/sqrt(rn);
w  = 2.0/3;
q  = 0.5;
b  = 1.15;
h  = 2*pi/(l*w);
od = 2*pi/(n*h*w*w);

y[0][0] = 0;
y[1][0] = 2;

/* Using the Runge-Kutta algorithm to integrate the equation */

for (i = 0;  i < n-1; ++i)
  {
  t = h*(i+1);
  y1 = y[0][i];
  y2 = y[1][i];
  dk11 = h*g1(y1,y2,t);
  dk21 = h*g2(y1,y2,t);
  dk12 = h*g1((y1+dk11/2),(y2+dk21/2),(t+h/2));
  dk22 = h*g2((y1+dk11/2),(y2+dk21/2),(t+h/2));
  dk13 = h*g1((y1+dk12/2),(y2+dk22/2),(t+h/2));
  dk23 = h*g2((y1+dk12/2),(y2+dk22/2),(t+h/2));
  dk14 = h*g1((y1+dk13),(y2+dk23),(t+h));
  dk24 = h*g2((y1+dk13),(y2+dk23),(t+h));
  y[0][i+1] = y[0][i]+(dk11+2*(dk12+dk13)+dk14)/6;
  y[1][i+1] = y[1][i]+(dk21+2*(dk22+dk23)+dk24)/6;

/* Bring theta back to the region [-pi,pi] */
 
  np = y[0][i+1]/(2*pi)+0.5;
  y[0][i+1] = y[0][i+1]-2*pi*np;
  }
for (i = 0; i < n; ++i)
  {
  ar[i] = y[0][i];
  wr[i] = y[1][i];
  ai[i] = 0;
  wi[i] = 0;
  }
fft (ar,ai,n,m);
fft (wr,wi,n,m);

for (i = 0; i < n; ++i)
  {
  o[i] = i*od;
  ar[i] = pow(f1*ar[i],2)+pow(f1*ai[i],2);
  wr[i] = pow(f1*wr[i],2)+pow(f1*wi[i],2);
  ar[i] = log10(ar[i]);
  wr[i] = log10(wr[i]);
  }
for (i = 0; i < l*m; i = i+4)
  {
  printf("%16.8lf %16.8lf %16.8lf\n", o[i],ar[i],wr[i]);
  }
}

void fft (ar,ai,n,m)
/* An example of the fast Fourier transform subroutine with
   N = 2**M.  AR and AI are the real and imaginary part of
   data in the input and corresponding Fourier coefficients
   in the output.  Copyright (c) Tao Pang 1997. */
int n,m;
double ar[],ai[];
{
int i,j,k,l,n1,n2,l1,l2;
double pi,a1,a2,q,v,u;

pi = 4*atan(1);
n2 = n/2;

n1 = pow(2,m);
if (n1 != n)
  {
  printf("Indices do not match\n");
  exit(1);
  }

/* Rearrange the data to the bit reversed order */

l = 1;
for (k = 0; k < n-1; ++k)
  {
  if (k < l-1)
    {
    a1 = ar[l-1];
    a2 = ai[l-1];
    ar[l-1] = ar[k];
    ar[k] = a1;
    ai[l-1] = ai[k];
    ai[k] = a2;
    }
  j = n2;
  while (j < l)
    {
    l = l-j;
    j = j/2;
    }
  l = l+j;
  }

/* Perform additions at all levels with reordered data */

l2 = 1;
for (l = 0; l < m; ++l)
  {
  q = 0;
  l1 = l2;
  l2 = 2*l1;
  for (k = 0; k < l1; ++k)
    {
    u =  cos(q);
    v = -sin(q);
    q =  q+pi/l1;
    for (j = k; j < n; j = j+l2)
      {
      i = j+l1;
      a1 = ar[i]*u-ai[i]*v;
      a2 = ar[i]*v+ai[i]*u;
      ar[i] = ar[j]-a1;
      ar[j] = ar[j]+a1;
      ai[i] = ai[j]-a2;
      ai[j] = ai[j]+a2;
      }
    }
  }
}