/* This file contains all the necessary functions to solve an n-th order */
/* ordinary differential equation using the Runge-Kutta method. */
/* The equation must be converted to a system of n first order differential */
/* equations, which are then aproximated using a fourth-order Runge-Kutta. */
/* The equations are supposed to be in ascending order of derivative, i.e. */
/* the first equation should correspond to the first derivative of the wanted */
/* solution function, the second equation to the second derivative, etc. */
/* Thus, "rungekutta" returns the solution for the first equation, which */
/* is also the solution of the n-th order differential equation. */

/* Example : */
/* W''(x) + aW'(x) + bW(x) = cF(x) */
/* Y1(x) = W(x) */
/* Y2(x) = W'(x) */
/* System of equations (in correct order) : */
/* Y1'(x) = Y2(x) */
/* Y2'(x) = -aY2(x) - bY1(x) + cF(x) */

# include <stdio.h>
# include <stdarg.h>

# include "rungekutta.h"

/* IMPLEMENTATION */

/******************************************************************************/
/* initialvalues : this function MUST always be called before any use of      */
/*      the Runge-Kutta algorithm. It initializes the order of the equation   */
/*      to be solved (and thus of the system of 1. order eqs.), the initial   */
/*      value of the independent variable, and the vector of solutions (or    */
/*      initial values) for the equations, called 'vectsols'.                 */
/*      If it is not called, execution will be halted with an error report.   */
/* Inputs :                                                                   */
/*      order : the order of the equation to be solved (al least 1).          */
/*      initx : the initial value of the independent variable of the equation.*/
/*              This value is the left limit of the interval over which       */
/*              solution is desired : [initx, ?].                             */
/*      ... : a list of exactly (order) initial values for the (order)        */
/*            first order equations of the system. These values are assumed   */
/*            to be of type 'double'.                                         */
/******************************************************************************/
void initialvalues(int order, double initx, ...)
{
  va_list pt;
  int i;
  struct lnode *t;

  va_start(pt, initx);
  xi = initx;
  n = order;
  if ((vectsols = t = (struct lnode *) malloc(sizeof(struct lnode))) == NULL)
    {
      fprintf(stderr, "initialvalues: error in call to malloc.\n");
      exit(1);
    }
  for (i=1; i<=n; i++)
  {
    t->y = va_arg(pt, double);
    if (i == n)
      t->next = NULL;
    else if((t->next = (struct lnode *)malloc(sizeof(struct lnode))) == NULL)
           {
             fprintf(stderr, "initialvalues: error in call to malloc.\n");
             exit(1);
           }
    t = t->next;
  }
  va_end(pt);
}

/* NOT YET IMPLEMENTED */
double checkinterval(double h, double y)
{
}

/******************************************************************************/
/* yi : should be called from the individual functions that code the n first  */
/*      order differential equations. Let xi be the independent variable, and */
/*      y1, y2, ..., yn be the solutions at xi for the Y1,...,Yn DE. When "yi"*/
/*      is called, it returns a specified value inside the solutions vector   */
/*      y1,...,yn.                                                            */
/* Input : i : index between 1 and n used to access the desired yi.           */
/* Output : the value of the asked yi.                                        */
/******************************************************************************/
double yi(int i)
{
  int j;
  struct lnode *p;

  if (i < 1 || i > n)
    {
      fprintf(stderr, "yi: index of asked value is out of bounds : %i.\n",i);
      exit(2);
    }
  j = 1;
  p = vectsols;
  while (j++ != i)
    p = p->next;
  return(p->y);
}

void chgyi(int i, double v)
{
  int j;
  struct lnode *p;

  if (i < 1 || i > n)
    {
      fprintf(stderr, "yi: index of asked value is out of bounds : %i.\n",i);
      exit(2);
    }
  j = 1;
  p = vectsols;
  while (j++ != i)
    p = p->next;
  p->y = v;
}

/******************************************************************************/
/* initialize : called by "rungekutta" once every step. It creates a copy 'ys'*/
/*      of the solutions vector 'vectsols', which contains the values         */
/*      y1,...,yn. 'ys' is used in the Runge-Kutta algorithm. When creating   */
/*      it, the function also checks for errors, such as having more initial  */
/*      values than equations or viceversa. "initialize" also creates two     */
/*      other vectors of size n, 'kj' and 'phij', used in the algorithm.      */
/******************************************************************************/
void initialize(void)
{
  struct lnode *t1, *t2, *t3, *t4;
  int i;

  t1 = vectsols;
  if ((ys = t2 = (struct lnode *) malloc(sizeof(struct lnode))) == NULL)
    {
      fprintf(stderr, "initialize: error in call to malloc.\n");
      exit(1);
    }
  if ((kj = t3 = (struct lnode *) malloc(sizeof(struct lnode))) == NULL)
    {
      fprintf(stderr, "initialize: error in call to malloc.\n");
      exit(1);
    }
  if ((phij = t4 = (struct lnode *) malloc(sizeof(struct lnode))) == NULL)
    {
      fprintf(stderr, "initialize: error in call to malloc.\n");
      exit(1);
    }
  for (i=1; i<=n; i++)
    {
	      /* creates the vector ys, used to hold a copy of vectsols */
     if (((t1 = t1->next) == NULL) && (i < n))
	{
          fprintf(stderr, "initialize: %i equations and %i initial values.\n",n,i);
          exit(2);
        }
      if (t1 != NULL)
        {
          if ((t2->next = (struct lnode *)malloc(sizeof(struct lnode))) == NULL)
            {
              fprintf(stderr, "initialize: error in call to malloc.\n");
              exit(1);
            }
          t2 = t2->next;
        }
      else
        {
          t2->next = NULL;
        }
            /* creates the vectors to hold values of Kj and Phij */
      if (i == n)
        break;
      if ((t3->next = (struct lnode *) malloc(sizeof(struct lnode))) == NULL)
        {
          fprintf(stderr, "initialize: error in call to malloc.\n");
          exit(1);
        }
      t3 = t3->next;
      if ((t4->next = (struct lnode *) malloc(sizeof(struct lnode))) == NULL)
        {
          fprintf(stderr, "initialize: error in call to malloc.\n");
          exit(1);
        }
      t4 = t4->next;
    }
  if (t1 != NULL)
    {
      fprintf(stderr, "initialize: %i equations and %i initial values.\n",n,i);
      exit(2);
    }
  t3->next = t4->next = NULL;
}

/******************************************************************************/
/* freevalues : MUST be called after having finished using the Runge-Kutta    */
/*      algorithm. It frees the memory occupied by the vector of solutions    */
/*      'vectsols', which is first created by "initialvalues", and kept in    */
/*      memory during the execution of the algorithm.                         */
/******************************************************************************/
void freevalues(void)
{
  struct lnode *p;

  for (; vectsols != NULL; vectsols = p)
    {
      p = vectsols->next;
      free(vectsols);
    }
}

/******************************************************************************/
/* rungekutta : contains the Runge-Kutta algorithm for solution of a system   */
/*      of n first order ordinary differential equations. It is intended to be*/
/*      used to solve an n-th order DE, which has been decomposed into a      */
/*      system of n 1. order equations. It assumes that 'initvalues' has been */
/*      called once before. Every time it is called, it executes the algorithm*/
/*      for one step.                                                         */
/* Inputs :                                                                   */
/*      h : the value of the distance between steps. If xi is the value of    */
/*          the independent variable at the beginning of the step, then       */
/*          x(i+1) = xi + h is its value at the end of the step.              */
/*      ... : a list of exactly n pointers to functions, where n is the order */
/*          of the system. Assuming that Y1,...,Yn are the n first order DE   */
/*          that compose the system, there are supposed to be FY1,...,FYn     */
/*          functions that code the respective equations. In the list, the    */
/*          functions must be in ascending order of "derivative", i.e. FY1    */
/*          must correspond to the first derivative of the solution function  */
/*          Y, FY2 to the second derivative, etc.                             */
/*               The declaration of the FYi functions must be as follows :    */
/*                                                                            */
/*                     double FYi(double xi)                                  */
/*                                                                            */
/*          where xi is the value of the independent variable at the beginning*/
/*          of the step. FYi must return a 'double' value. When FYi needs to  */
/*          refer to a value in the solutions vector y1,...,yn, it should     */
/*          access it using "yi(i)", where i is the index of the desired      */
/*          value.                                                            */
/*               "rungekutta" should be called in the following way :         */
/*                                                                            */
/* x = rungekutta(h, (double (*)(double))FY1,..., (double (*)(double))FYn);   */
/*                                                                            */
/* Output : the computed solution of the n-th order DE at the next point,     */
/*          x(i+1) = x(i) + h.                                                */
/******************************************************************************/
double rungekutta(double h, ...)
{
  va_list pt;
  int j;
  double (*equation) (double);
  struct lnode *k, *phi, *tsols, *ty;

  initialize();

  va_start(pt, h);
  k = kj;
  tsols = vectsols;
  phi = phij;
  ty = ys;
  for (j=1; j<=n; j++)
    {
      equation = (double (*) (double)) va_arg(pt, void *);
      ty->y = tsols->y;
      phi->y = k->y = (*equation) (xi);
      k = k->next;
      tsols = tsols->next;
      phi = phi->next;
      ty = ty->next;
    }
  tsols = vectsols;
  k = kj;
  ty = ys;
  for (j=1; j<=n; j++)
    {
      tsols->y = ty->y + h*k->y/2;
      tsols = tsols->next;
      k = k->next;
      ty = ty->next;
    }
  xi += h/4.0;
  va_end(pt);

  va_start(pt, h);
  k = kj;
  phi = phij;
  for (j=1; j<=n; j++)
    {
      equation = (double (*) (double)) va_arg(pt, void *);
      k->y = (*equation) (xi);
      phi->y += 2*k->y;
      k = k->next;
      phi = phi->next;
    }
  k = kj;
  tsols = vectsols;
  ty = ys;
  for (j=1; j<=n; j++)
    {
      tsols->y = ty->y + h*k->y/2;
      k = k->next;
      tsols = tsols->next;
      ty = ty->next;
    }
    xi+=h/4.0;
  va_end(pt);

  va_start(pt, h);
  k = kj;
  phi = phij;
  for (j=1; j<=n; j++)
    {
      equation = (double (*) (double)) va_arg(pt, void *);
      k->y = (*equation) (xi);
      phi->y += 2*k->y;
      k = k->next;
      phi = phi->next;
    }
  k = kj;
  tsols = vectsols;
  ty = ys;
  for (j=1; j<=n; j++)
    {
      tsols->y = ty->y + h*k->y;
      k = k->next;
      tsols = tsols->next;
      ty = ty->next;
    }
  xi += h/4.0;
  va_end(pt);

  va_start(pt, h);
  k = kj;
  for (j=1; j<=n; j++)
    {
      equation = (double (*) (double)) va_arg(pt, void *);
      k->y = (*equation) (xi);
      k = k->next;
    }
  k = kj;
  phi = phij;
  tsols = vectsols;
  ty = ys;
  for (j=1; j<=n; j++)
    {
      tsols->y = ty->y + h*(phi->y + k->y)/6;
      k = k->next;
      phi = phi->next;
      tsols = tsols->next;
      ty = ty->next;
    }
    xi+=h/4.0;
  va_end(pt);

  for (; kj != NULL; kj = k)
    {
      k = kj->next;
      free(kj);
    }
  for (; phij != NULL; phij = phi)
    {
      phi = phij->next;
      free(phij);
    }
  for (; ys != NULL; ys = ty)
    {
      ty = ys->next;
      free(ys);
    }
  return (vectsols->y);
}