#include "odes.h"
/***********************************************************************
 *   FILENAME:      odes.c                                             *
 *                                                                     *
 *   FUNCTION:      euler                                              *
 *                                                                     *
 *   PROGRAMMER:    Izzi Urieli                                        *
 *                                                                     *
 *   DATE:          12/27/94                                           *
 *                                                                     *
 *   PROTOTYPE:                                                        *
 *     void                                                            *
 *     euler(void deriv(double x, double y[], double dy[]),            *
 *           int n, double *x, double dx, double y[], double dy[]);    *
 *                                                                     *
 *   INCLUDE:       "odes.h"                                           *
 *                                                                     *
 *   GLOBAL VARIABLES: none                                            *
 *                                                                     *
 *   DESCRIPTION: integrates the vector of first order ordinary        *
 *      differential equations over the increment dx by the            *
 *      simple first order Euler method.                               *
 ***********************************************************************/
void
euler(void deriv(double x, double y[], double dy[]),
      int n, double *x, double dx, double y[], double dy[])
{
   int i;
   double x0, x1;

   x0 = *x;
   deriv(x0, y, dy);
   for(i = 0; i < n; i++)
      y[i] = y[i] + dx * dy[i];
   x1 = x0 + dx;
   *x = x1;
}
/***********************************************************************
 *   FILENAME:      odes.c                                             *
 *                                                                     *
 *   FUNCTION:      rk4                                                *
 *                                                                     *
 *   PROGRAMMER:    Izzi Urieli                                        *
 *                                                                     *
 *   DATE:          12/27/94                                           *
 *                                                                     *
 *   PROTOTYPE:                                                        *
 *     void                                                            *
 *     rk4(void deriv(double x, double y[], double dy[]),              *
 *         int n, double *x, double dx, double y[], double dy[]);      *
 *                                                                     *
 *   INCLUDE:       "odes.h"                                           *
 *                                                                     *
 *   SYMBOLIC CONSTANTS:  NEQU 25                                      *
 *   (maximum number of first order differential equations in the set) *
 *                                                                     *
 *   GLOBAL VARIABLES: none                                            *
 *                                                                     *
 *   DESCRIPTION: integrates the vector of first order ordinary        *
 *      differential equations over the increment dx by the            *
 *      Classical fourth order Runge-Kutta method                      *
 ***********************************************************************/
void
rk4(void deriv(double x, double y[], double dy[]),
    int n, double *x, double dx, double y[], double dy[])
{
   int i;
   double x0, xm, x1;
   double y0[NEQU], dy1[NEQU], dy2[NEQU], dy3[NEQU], dy4[NEQU];

   x0 = *x;
   for(i = 0; i < n; i++)
      y0[i] = y[i];

   deriv(x0, y, dy1);
   for(i = 0; i < n; i++)
      y[i] = y0[i] + 0.5 * dx * dy1[i];
   xm = x0 + 0.5 * dx;

   deriv(xm, y, dy2);
   for(i = 0; i < n; i++)
      y[i] = y0[i] + 0.5 * dx * dy2[i];

   deriv(xm, y, dy3);
   for(i = 0; i < n; i++)
      y[i] = y0[i] + dx * dy3[i];

   x1 = x0 + dx;
   deriv(x1, y, dy4);
   for(i = 0; i < n; i++) {
      dy[i] = (dy1[i] + 2.0 * (dy2[i] + dy3[i]) + dy4[i]) / 6.0;
      y[i] = y0[i] + dx * dy[i];
   }
   *x = x1;
}