/* sim_difeq.c  solve simultaneous differential equations             */
/*                by Runge-Kutta by Gills method                      */
/*                                                                    */
/* solve system of equations Y_i = Y'_i(X, Y_1, ..., Y_n) for i=1,n   */
/*                                                                    */
/* user sets N number of equations (in N dependent variables)         */
/*           Y array of N dependent variables with initial values     */
/*           F array of N first derivatives of Y's (computed in loop) */
/*           Q array of N for use by Runge                            */
/*           X independent variable with initial value                */
/*           XLIM limit value of X with value                         */
/*           H step size, delta X with value                          */
/*           M setp counter, initialize to zero                       */
/*           K returned 1=compute derivatives, 2=have Y values at X+H */
#include <stdio.h>

int Runge(int N, double Y[], double F[], double Q[],
          double *X,  double H); /* Gills method */

int main(int argc, char *argv[])
{
  const int N = 4;   /* number of equations      */
  double Y[N];       /* dependent variables      */
  double F[N];       /* computed derivatives, Y' */
  double Q[N];       /* temporaries for Runge    */
  double X, XLIM, H;
  int i, K;

  printf("sim_difeq.c solving %d equations \n", N);
  X = 0.0;      /* initial independent variable */
  XLIM = 1.0;   /* final value of independent variable */
  H = 0.1;      /* step size */
                /* initial condition for dependent variables */
  for(i=0; i<N; i++) Y[i] = 0.0; 
	
  printf("simple numeric check, answer 1, 2, 3, 4 \n");
  while(X < XLIM) /* simple numeric check */
  {
    K = 1;
    while(K == 1)
    {
      K = Runge(N, Y, F, Q, &X, H); /* X, M local, K returned */
      F[0] = 1.0;
      F[1] = 2.0;
      F[2] = 3.0;
      F[3] = 4.0;
    } /* X incremented by H */
    printf("X=%18.15f Y[0]=%18.15f Y[1]=%18.15f \n", X, Y[0], Y[1]);
    printf("                     Y[2]=%18.15f Y[3]=%18.15f \n", Y[2], Y[3]);
  } /* X at XLIM (done) */

  printf("\n");
  printf("fourth order check 1, 1/2, 1/6, 1/24 \n");
  X = 0.0;      /* initial independent variable */
  XLIM = 1.0;   /* final value of independent variable */
  H = 0.01;     /* step size */
                /* initial condition for dependent variables */
  for(i=0; i<N; i++) Y[i] = 0.0; 
	
  while(X < XLIM) /* fourth order check */
  {
    K = 1;
    while(K == 1)
    {
      K = Runge(N, Y, F, Q, &X, H);
      F[0] = 1.0;
      F[1] = Y[0];
      F[2] = Y[1];
      F[3] = Y[2];
    } /* X incremented by H */
  } /* X at XLIM (done) */
  printf("X=%18.15f Y[0]=%18.15f  Y[1]=%18.15f \n", X, Y[0], Y[1]);
  printf("                     Y[2]=%18.15f  Y[3]=%18.15f \n", Y[2], Y[3]);
  printf("Y[3]*24.0 = %18.15f \n", Y[3]*24.0);
  printf("\n");
  printf("sim_difeq.c finished \n");
  return 0;
} /* end main of sim_difeq */

int Runge(int N, double Y[], double F[], double Q[],
          double *X,  double H) /* Gills method */
{
  double A;
  int i, K;
  static int M = 0;

  M++;     /* integration step */
  K = 1;   /* exit condition */
  if(M == 1)
  {
    for(i=0; i<N; i++) Q[i] = 0.0;
  }
  else if(M == 2)
  {
    A = 0.5;
    *X = *X+0.5*H;
    for(i=0; i<N; i++)
    {
      Y[i] = Y[i] + A*F[i]*H;
      Q[i] = 2.0*A*H*F[i];
    }
  }
  else if(M == 3)
  {
    A = 0.2928932188134524756;
    for(i=0; i<N; i++)
    {
      Y[i] = Y[i] + A*(F[i]*H-Q[i]);
      Q[i] = 2.0*A*H*F[i] + (1.0-3.0*A)*Q[i];
    }
  }
  else if(M == 4)
  {
    A = 1.7071067811865475244;
    *X = *X+0.5*H;
    for(i=0; i<N; i++)
    {
      Y[i] = Y[i] + A*(F[i]*H-Q[i]);
      Q[i] = 2.0*A*H*F[i] + (1.0-3.0*A)*Q[i];
    }
  }
  else if(M == 5)
  {
    A = 0.2928932188134524756;
    for(i=0; i<N; i++)
    {
      Y[i] = Y[i] + H*F[i]/6.0 - Q[i]/3.0;
    }
    M = 0;
    K = 2;
  }
  else
  {
    printf("ERROR in Runge, bad M=%d \n", M);
    M = 0;
    K = 2;
  }
  return K;
} /* end Runge */
/* end sim_difeq.c */