/******
 * A program to look at
 * a chaotic differential equation
 #    compile with "cc -lm chaos.c -o chaos"
 #    run with "chaos > chaos.out"
 ******

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

/* diffeq to solve, expressed as system of linear eqns 
    Here I'm doing a second order eqn, namely
          x''(t) = - A sin( B x) + C sin( t)  - B x'(t)
    where A,B,C are constants.
    We change this into a system of two first order eqns, 
    namely 
         x' = dxdt( x, v, t) = v
         v' = dvdt( x, v, t) = -A dsin( B x) + C dsin(t) - B v
    and we have functions that calculate dxdt and dvdt.
    I'm going to write this all out, as if we were doing
    the general case.  */

/* These are global variables - you can 
   see them or set them from anywhere in this file.
   But that means you can't use these variables anywhere else. */
double A;
double B;
double C;

double dxdt( double x, double v, double t ) {
  return(v);
}
double dvdt( double x, double v, double t) {
  return( -A*dsin(x) + B*dsin(t) - C*v );
}

/* Now, here's a simple engine to look for 
   a solution with a fixed timestep dt, for N points,
   starting at a given x0, v0, t0 */
void simpleSolver( int N, double x0, double v0, double t0, double dt, 
		   double x[], double v[] ){
  int i;
  x[0] = x0;
  v[0] = v0;
  for (i=0;i<(N-1);i++) {
    x[i+1] = x[i] + dt * dxdt(x[i], v[i], t0+dt*i);
    v[i+1] = v[i] + dt * dvdt(x[i], v[i], t0+dt*i);
  }
}


/* Stuff to try: 
     - implement the runge-kutta algorithm in 
       the text, pg 62, and see how the steps compare
       when you try to get similar results. */

#define NPOINTS 1000

main() {
  double x[NPOINTS];
  double v[NPOINTS];
  double dt;  
  int N = 500;
  double x0, v0, t0;
  int i;

  /* Set the constants to some values */
  A = 0.5;   /* =   omega^2/forcing_freq^2  */
  B = 0.2;   /* =   amplitude of forcing term  */
  C = 0.1;   /* =   friction force  */

		     x0 = 0.2;
  v0 = -1;
  dt = 0.01;;
  simpleSolve( N, x0, v0, t0, x, v);

  for (i=0;i<N;i++){
    printf( "%g %g \n", x[i], v[i]);
  }

}