/*   This handles the delay stuff    */

#include <stdio.h>
#include <math.h>
#include "xpplim.h"

double AlphaMax=100,OmegaMax=100;

double *DelayWork;
int LatestDelay;
int MaxDelay;
int DelayFlag=0;

int NDelay,del_stab_flag,WhichDelay,DelayGrid=1000;
double variable_shift[2][MAXODE];
double delay_list[MAXDELAY];

double evaluate();
extern double DELTA_T,T0,DELAY;
extern int NODE,NCON,NSYM,NSYM_START,NCON_START,NMarkov;
 
extern char delay_string[MAXODE][80];
extern double variables[];



double delay_stab_eval(delay,var)  
/* this returns appropriate values for delay jacobian */ 
     double delay;
     int var;
{
  int i;

  if(del_stab_flag==0) /* search for all delays  */
    {
      for(i=0;i<NDelay;i++){
	if(delay==delay_list[i])
	  return(variables[var]);
      }
      delay_list[NDelay]=delay;
      NDelay++;
      return(variables[var]);
    }
 /*  now we must determine the value to return  */
 /*  del_stab_flag =-1 */    
     for(i=0;i<NDelay;i++){
       if(delay==delay_list[i])
	if(i==WhichDelay)
	  return variable_shift[1][var-1];
     }
   return variable_shift[0][var-1];
}



alloc_delay(big)
double big;
{
 int n,i;
 double zip;
 
 n=(int)(big/fabs(DELTA_T))+1;

 MaxDelay=n;
 LatestDelay=1;
 DelayFlag=0;
 DelayWork=(double *)calloc(n*(NODE ),sizeof(double));
 if(DelayWork==NULL){
  err_msg("Could not allocate memory for Delay");
  return(0);
 }
 DelayFlag=1;
 NDelay=0;
 WhichDelay=-1;
 del_stab_flag=1;
 for(i=0;i<n*(NODE );i++)DelayWork[i]=0.0;
 return(1);
}

free_delay()
{
 if(DelayFlag)free(DelayWork);
 DelayFlag=0;
}

stor_delay(y)
double *y;
{
 int i,in;
 if(DelayFlag==0)return;
 --LatestDelay;
 if(LatestDelay<0)LatestDelay+=MaxDelay;
 in=LatestDelay*(NODE );
 for(i=0;i<(NODE );i++)DelayWork[i+in]=y[i];
}

double get_delay(in,tau)
int in;
double tau;
{
 double x=tau/fabs(DELTA_T);
 int n1=(int)x;
 int n2=n1+1;
 int i1,i2;
 double x1,x2;
 if(tau<0.0||tau>DELAY){
			 err_msg("Delay negative or too large");
			stop_integration();
			return(0.0);
			}
  i1=(n1+LatestDelay)%MaxDelay;
  i2=(n2+LatestDelay)%MaxDelay;
  x1=DelayWork[in+(NODE )*i1];
  x2=DelayWork[in+(NODE )*i2];
  return(x1+(x-n1)*(x2-x1));
 }

/*  Handling of the initial data  */
do_init_delay(big)
double big;
{
 double t=T0,old_t,y[MAXODE];
 int i,nt,j;
 int len;
 
 int *del_form[MAXODE];
 nt=(int)(big/fabs(DELTA_T));
 NCON=NCON_START;
 NSYM=NSYM_START;
 for(i=0;i<(NODE );i++){
	del_form[i]=(int *)calloc(200,sizeof(int));
	if(del_form[i]==NULL){
		err_msg("Failed to allocate delay formula ...");
		for(j=0;j<i;j++)free(del_form[j]);
                NCON=NCON_START;
		NSYM=NSYM_START;
		return(0);
		}
	 if(add_expr(delay_string[i],del_form[i],&len)){
		err_msg("Illegal delay expression");
                for(j=0;j<=i;j++)free(del_form[j]);
		 NCON=NCON_START;
		NSYM=NSYM_START;
		return(0);
		}
	 }        /*  Okay all formulas are cool... */
  LatestDelay=1;
	
  get_val("t",&old_t);
 
  for(i=nt;i>=0;i--){
	t=T0-fabs(DELTA_T)*i;
	set_val("t",t);
	for(j=0;j<(NODE );j++)
		y[j]=evaluate(del_form[j]);
	stor_delay(y);
  }
   for(j=0;j<(NODE );j++)free(del_form[j]);
   NCON=NCON_START;
   NSYM=NSYM_START;
  set_val("t",old_t);
   return(1);
 }
		
	 
	
        
 
 
  

 
 
 

 






