/* Runge-Kutta simulation of cardiovascular hemodynamics.
 * Timothy Davis, May 1989.  Previous versions by R.Sah and G.Moody.
 * Copyright 1988, 1989 Timothy L. Davis
 * $Source: /mit/hst/src/sim/RCS/simulate.c,v $
 * 
 * Run the CV simulation for a specific number of cycles.  Writes simulation 
 * data to extern procedure SimDataEnq: (6 pressures + time + Cl + Cr + 
 * capacitance array index)
 * simulate() takes original pressure/time/capacitance vector (x),& # of beats.
 * This simulator expects parameters to be set up (eqns.c:initeqnvars()),
 * and for initialization of the pressure array via estimate.c:estimate() or
 * previous call to simulate().
 */
#include <signal.h>
#include <math.h>
#include "CVDefs.h"
#include "../event/Event.h"
#define ADAPT_RK /* Use adaptive step size (quality control) Runge-Kutta */

typedef void (*sighandler_t)(int);

extern double Tsystole(), Cl(), Cr(), Cinc(); /* all from file eqns.c */
extern double getval();
extern void eqns();

/*PUBLIC*/
int simulate(x, ticks)
double x[]; /* starting values of 6 pressures, time, Cl, Cr, and
		        current index into bthr capacitance table */
int ticks; /* number of consecutive times steps desired */
{
  int j;
  double xsav[10]; /* private copy of pressure data for speed */
  double dxdt[6];		/* dx/dt[], evaluated by "eqns" */
#ifndef ADAPT_RK
  double dt;			/* time increment */
#endif
#ifdef SVR4
  void (*old_handler)(int) = 
    sigset(SIGFPE, fpe_in_simulate); /* from Event.c */
#else /* assuming BSD */
  int (*old_handler)()=(void *)signal(SIGFPE,(sighandler_t)fpe_in_simulate); /* from Event.c */
#endif

  for (j=0; j<10; j++) xsav[j] = x[j];

  for (j = 0; j < ticks; j++) { /* do ticks Runge-Kutta segments */
    initvalves(x);          /* set valves depending on pressures */
    fixvolume(x);		/* correct at sys venous reservoir */

#ifndef ADAPT_RK
    dt =  6*Tsystole(getval(x, SYSTEM, HEARTRATE))/SYSTOLESEGS;
    sanode(x, dt); /* determine onset of contraction and update constants */
    eqns(dxdt, (int) x[9], x);
    rk4(x, dxdt, x[9], dt, x, eqns);
    x[6] += dt;
    x[9] += Cinc(dt);
    x[7] = Cl((int) x[9]);
    x[8] = Cr((int) x[9]);
    if (fabs(x[0])<1000 && fabs(x[1])<1000 && fabs(x[2])<1000
	&& fabs(x[3])<1000 && fabs(x[4]) < 1000 && fabs(x[5]) < 1000)
      SimDataEnq(x);
    else {
#ifdef SVR4
    sigset(SIGFPE,old_handler); /* from Event.c */
#else /* assuming BSD */
    signal(SIGFPE, old_handler);
#endif
      for (j=0; j<10; j++) x[j] = xsav[j];
      return(0); /* FALSE to restart simulation attempt */
    }
#else
    {
      static double hdid = 0.000;
      static double htry = 0.003;
      static double yscale[] = {1, 1, 1, 1, 1, 1}; 
      double hnext;

      sanode(x, hdid); /* determine onset of contraction */
      eqns(dxdt, (int) x[9], x);
      if (!rkqc(x, dxdt, &x[9], htry, 0.01, 
		yscale, &hdid, &hnext, eqns)) {
#ifdef SVR4
	sigset(SIGFPE,old_handler);
#else
	signal(SIGFPE,(sighandler_t)old_handler);
#endif
	for (j=0; j<10; j++) x[j] = xsav[j];
	return(0); /* Convergence failure! */
      }
      htry = hnext > 0.01  ? 0.01 : hnext; /* Keep arterial diastole smooth */
      x[6] += hdid;
      x[7] = Cl((int)x[9]);
      x[8] = Cr((int)x[9]);
      SimDataEnq(x); 
    }
#endif
  }
#ifdef SVR4
  sigset(SIGFPE,old_handler);
#else
  signal(SIGFPE, (sighandler_t)old_handler);
#endif
  return(1); /* successful... done */
}