/* Template file for coding d.d.e. models for solv95 - Example is a switched dde
   model

*/



#include <math.h>
//#include <windows.h>
#include "solv95.h"

/***************************************************************************/
/*  Put global variables  here. These should never be written to from      */
/* grad() or switchfunctions(), directly or indirectly.                  	*/
/***************************************************************************/



/***************************************************************************/
/*             Problem specific routines           								*/
/***************************************************************************/

void initcons(no_vars,no_cons) int *no_cons,*no_vars;
/* this routine specifies the number of constants and number of variables in
	the model.
*/
{ *no_cons=6;
  *no_vars=2;
}

void switchfunctions(sw,s,c,t)
double *sw,*s,*c,t;


/* This routine sets the values of the switch functions. When the switch
	functions pass through zero from positive to negative the state variables
	may be reset in function map(). The switch functions should pass smoothly
	through 0 and should never have both value and first derivative zero. The
	same switch must not pass through zero from positive to negative more than
	once in an integration timestep. An example of a switch function is:
						sw[0]=sin(pi*t/30.0)
	which passes through zero every 60 time units. Switches may include state
	variables provided the above conditions are met. Note that to get 'Solver'
	style switches define twice as many switches and let e.g. sw[1]=-sw[0] */


{ static double pi;
  static int first=1;
  if (first) {first=0;pi = 2*asin(1.0);}
  sw[0]=sin(2*pi*t/c[1]); // add resource
  sw[1]=sin(2*pi*(t-c[4])/c[1]); // stop integrator as discontinuity in lagged variable occurs
}



void map(s,c,t,swno)
double *s,*c,t;int swno;

/* This routine is called whenever one of the switch functions passes through
	zero. 'swno' is the number of the switch function. The state variables
	can be changed discontinuously within this routine. eg:
   if (swno==1)
	  { s[0]=coeff[1]*(s[0]);}
	time and the coefficients should not be changed.
*/

{ if (swno==0)
  { s[0]+=c[2]; // adding resources
  }             // otherwise nothing required
}




void grad(g,s,c,t,pastvalue)
double *g,*s,*c,t;
double (*pastvalue)();

/* This routine must provide the gradients g for the state variables s.
	So ds[i]/dt=g[i]=fi(s,c,t) where c is the coefficient vector. lagged
   variables may be accessed here using pastvalue(i,x,j) which returns the
   ith (starting at zero) lagged variable at time x, using lag pointer k

   (lag pointers are used by pastvalue to store the history buffer location
    corresponding to a lag in order to save exectution time. For example if
    your code requires lagged varaible 0 at lags T and 2T for each gradient
    calculation then it is efficient to obtain these values using:
    pastvalue(0,t-T,0) and pastvalue(0,t-2*T,1) rather than
    pastvalue(0,t-T,0) and pastvalue(0,t-2*T,0). The latter works, it's just
    slower because more time is spent searching for lagged values)
*/

{ g[0]=-c[0]*s[0]*s[1]; // resource
  if (t>c[4])
  g[1]=c[3]*pastvalue(0,t-c[4],0)*pastvalue(1,t-c[4],0)-c[5]*s[1]; // consumer
  else
  g[1]=-c[5]*s[1];
}

void storehistory(his,ghis,g,s,c,t)
double *his,*ghis,*g,*s,*c,t;

/* This is the routine in which the values of the history variables at time
	t are calculated and put in the array his, along with gradients in ghis,
	using state variables s, gradients of s, g, and coefficients c
   e.g. if the state variable 2 is history variable 0, you would need the line:
   his[0]=s[2];ghis[0]=g[2];
*/

{ his[0]=s[0];ghis[0]=g[0];
  his[1]=s[1];ghis[1]=g[1];
}



void statescale(double *scale)

/* In this routine you can set scale factors for error control. For each
   state variable the maximum permisable error will be bounded below by the
   tolerance multiplied by scale[i]. If you don't supply values then zero will
   be used.
	Non-zero scale values are useful for variables that start at zero and
   leave zero without 3rd order continuity. */

{
}


void initst(s,c,t)
double *s,*c,t;

/* initialise state variables and any global constants here, you can use c */

{ s[0]=0.0;s[1]=1.0;

}

void initout(out) usercontrol *out;

/* This routine is where the output windows are first set up. You have to
   fill in the global structure "out", which tells solv95 which
   variables to plot and in which windows */




{ int i;
  /* how many windows */
  out->no_windows=2;

  /* how many lines in each window e.g. for two lines in window zero set
     out->lines[0]=2 ..... */
  out->lines[0]=1;
  out->lines[1]=1;

  /* which state variables go in which windows and what are they called?
     e.g. to put state variable 4 in window 3 as curve 0, and label it `trash':
     out->index[4].win=3;out->index[4].cur=0;out->label[4]="trash";
     If you don't provide this information for a variable then it isn't plotted
  */
  out->label[0]="Resource";
  out->index[0].win=0;out->index[0].cur=0;

  out->label[1]="Consumer";
  out->index[1].win=1;out->index[1].cur=0;

  /* labels for windows */
  out->wname[0]="Output";
  out->wname[1]="More Output";
  /* Now set up initial range of y axis variables for each window */
  out->range[0].y0=0.0;out->range[0].y1=100.0;
  out->range[1].y0=0.0;out->range[1].y1=200.0;

  out->xlabel="Time";
  /* initialise file output. Set out->fileno to the number of columns to output
     (in addition to time). Set out->fout[i] to the number of the state variable
     that is to go in column i. Set out->fileno=0 for no output. */
  out->fileno=2;
  out->fout[0]=0;
  out->fout[1]=1;

  /* integration details*/
  out->t0=0.0;        /* default start time */
  out->t1=300.0;     /* default stop time */
  out->dt=1.0;        /* initial timestep */
  out->tol=0.000005;  /* integration tolerance */
  out->dout = 1.0;    /* approximate average output timestep */
  out->hbsize=1000L;  /* how many past values to store for each history variable */
  out->nhv=2;         /* Number of history (lagged) variables */
  out->nlag=1;        /* Number of lag markers per history variable (set to 1 if unsure)*/
  out->nsw=2;         /* number of switch varaibles */

  /* Initial Values and names for constants to be prompted for at run time */

  out->c[0]=0.1; out->cname[0]="alpha";
  out->c[1]=10.0; out->cname[1]="T";
  out->c[2]=50.0;out->cname[2]="R_a";
  out->c[3]=0.05;out->cname[3]="gamma";
  out->c[4]=5.0;out->cname[4]="tau";
  out->c[5]=0.02;out->cname[5]="delta";
  /* Optional information text on parameters - split text between lines with `\'*/

  out->cinfo[0]="Consumption rate per consumer per unit resource.";
  out->cinfo[1]="Time between resource top up.";
  out->cinfo[2]="Amount by which resources are topped up.";
  out->cinfo[3]="Consumer per capita fecundity per unit resource";
  out->cinfo[4]="Consumer development time.";
  out->cinfo[5]="Consumer per capita death rate.";
  /* Text and title for initial welcome box (optional) */
  out->initialtitle="Switched DDE Model";
  out->initialtext=" This model is an example that came with solv95 ";

}