#include <stdlib.h>
/* this is the main integrator routine
for phase-plane
It takes the steps looks at the interrupts, plots and stores the data
It also loads the delay stuff if required
*/
/*
New stuff for 9/96 -- cvode added
cvode(command,y,t,n,tout,kflag,atol,rtol)
command =0 continue, 1 is start 2 finish
kflag is error < 0 is bad -- call cvode_err_msg(kflag)
call end_cv() to end it normally
on return y is new stuff, t is new time kflag is error if any
if kflag < 0 thats bad
NOTE: except for the structure MyGraph, it is "x-free" so it
is completely portable
*/
#include <stdio.h>
#include <X11/Xlib.h>
#include <math.h>
#include <string.h>
#include "xpplim.h"
#include "struct.h"
#include "phsplan.h"
extern GRAPH *MyGraph;
#include "menudrive.h"
#define MSWTCH(u,v) memcpy((void *)(u),(void *)(v),xpv.node*sizeof(double))
#define READEM 1
#define ESCAPE 27
#define FIRSTCOLOR 30
#define MAX_LEN_SBOX 25
#define PARAM 1
#define IC 2
#define BMAXCOL 20
#define MAXFP 400
#define NAR_IC 50
#define VOLTERRA 6
#define BACKEUL 7
#define RKQS 8
#define STIFF 9
#define GEAR 5
#define CVODE 10
#define DP5 11
#define DP83 12
#define RB23 13
extern int animation_on_the_fly;
extern double ShootIC[8][MAXODE];
extern int ShootICFlag;
extern int ShootIndex;
extern int SimulPlotFlag,current_pop,num_pops,ActiveWinList[];
typedef struct {
int index0,type;
char formula[256];
int n;
char var[20];
int j1,j2;
} ARRAY_IC;
int ar_ic_defined=0;
double ndrand48();
ARRAY_IC ar_ic[NAR_IC];
typedef struct
{
int n,flag;
double *x[MAXFP];
double *er[MAXFP];
double *em[MAXFP];
} FIXPTLIST;
FIXPTLIST fixptlist;
typedef struct
{
int n;
double tol;
double xlo[MAXODE],xhi[MAXODE];
} FIXPTGUESS;
FIXPTGUESS fixptguess;
typedef struct
{
int nvec,node;
double *x;
} XPPVEC;
XPPVEC xpv;
int SuppressBounds=0;
extern char *info_message,*ic_hint[],*sing_hint[];
extern int Xup,batch_range;
extern char batchout[256];
double atof();
extern int NMarkov,STOCH_FLAG;
extern int color_total,SCALEY,DCURY,PSFlag,PointRadius;
int DelayErr;
float **get_browser_data();
double get_ivar();
double MyData[MAXODE],MyTime;
int MyStart;
extern int DelayFlag,DCURY,NKernel;
int RANGE_FLAG;
extern int PAR_FOL,SHOOT;
extern char upar_names[MAXPAR][11];
extern double last_ic[MAXODE];
double LastTime;
extern double DELAY;
extern int R_COL;
extern int colorline[11];
extern (*rhs)();
int STOP_FLAG=0;
int PSLineStyle;
struct {
char item[30];
int steps,shoot,col,movie;
double plow,phigh;
} eq_range;
struct {
char item[30],item2[30];
int steps,steps2,reset,oldic,index,index2,cycle,type,type2,movie;
double plow,phigh,plow2,phigh2;
int rtype;
} range;
#define MAX_INTERN_SET 100
typedef struct {
char *name;
char *does;
} INTERN_SET;
extern INTERN_SET intern_set[MAX_INTERN_SET];
extern int Nintern_set;
int (*solver)();
init_ar_ic()
{
int i;
for(i=0;i<NAR_IC;i++){
ar_ic[i].index0=-1;
ar_ic[i].formula[0]=0;
ar_ic[i].n=0;
ar_ic[i].var[i]=0;
ar_ic[i].type=0;
}
}
dump_range(fp,f)
FILE *fp;
int f;
{
char bob[256];
if(f==READEM)
fgets(bob,255,fp);
else
fprintf(fp,"# Range information\n");
io_string(eq_range.item,11,fp,f);
io_int(&eq_range.col,fp,f,"eq-range stab col");
io_int(&eq_range.shoot,fp,f,"shoot flag 1=on");
io_int(&eq_range.steps,fp,f,"eq-range steps");
io_double(&eq_range.plow,fp,f,"eq_range low");
io_double(&eq_range.phigh,fp,f,"eq_range high");
io_string(range.item,11,fp,f);
io_string(range.item2,11,fp,f);
io_int(&range.steps,fp,f,"Range steps");
io_int(&range.cycle,fp,f,"Cycle color 1=on");
io_int(&range.reset,fp,f,"Reset data 1=on");
io_int(&range.oldic,fp,f,"Use old I.C.s 1=yes");
io_double(&range.plow,fp,f,"Par1 low");
io_double(&range.plow2,fp,f,"Par2 low");
io_double(&range.phigh,fp,f,"Par1 high");
io_double(&range.phigh2,fp,f,"Par2 high");
dump_shoot_range(fp,f);
if(f==READEM)range.steps2=range.steps;
}
init_range()
{
eq_range.col=-1;
eq_range.shoot=0;
eq_range.steps=10;
eq_range.plow=0.0;
eq_range.phigh=1.0;
eq_range.movie=0;
sprintf(eq_range.item,upar_names[0]);
range.type=0;
range.rtype=0;
range.index=range.index2=0;
range.steps=20;
range.steps2=20;
range.plow=range.plow2=0.0;
range.phigh=range.phigh2=1.0;
range.reset=1;
range.oldic=1;
range.cycle=0;
range.movie=0;
sprintf(range.item,uvar_names[0]);
sprintf(range.item2,uvar_names[0]);
init_shoot_range(upar_names[0]);
init_monte_carlo();
}
set_up_eq_range()
{
static char *n[]={"*2Range over","Steps","Start","End",
"Shoot (Y/N)",
"Stability col","Movie (Y/N)"};
char values[7][MAX_LEN_SBOX];
int status,i;
static char *yn[]={"N","Y"};
sprintf(values[0],"%s",eq_range.item);
sprintf(values[1],"%d",eq_range.steps);
sprintf(values[2],"%.16g",eq_range.plow);
sprintf(values[3],"%.16g",eq_range.phigh);
sprintf(values[4],"%s",yn[eq_range.shoot]);
sprintf(values[5],"%d",eq_range.col);
sprintf(values[6],"%s",yn[eq_range.movie]);
status=do_string_box(7,7,1,"Range Equilibria",n,values,45);
if(status!=0){
strcpy(eq_range.item,values[0]);
i=find_user_name(PARAM,eq_range.item);
if(i<0){
err_msg("No such parameter");
return(0);
}
eq_range.steps=atoi(values[1]);
if(eq_range.steps<=0)eq_range.steps=10;
eq_range.plow=atof(values[2]);
eq_range.phigh=atof(values[3]);
if(values[4][0]=='Y'||values[4][0]=='y')eq_range.shoot=1;
else eq_range.shoot=0;
if(values[6][0]=='Y'||values[6][0]=='y')eq_range.movie=1;
else eq_range.movie=0;
eq_range.col=atoi(values[5]);
if(eq_range.col<=1||eq_range.col>(NEQ+1))eq_range.col=-1;
return(1);
}
return(0);
}
cont_integ()
{
double tetemp;
double *x;
double dif;
if(INFLAG==0||FFT!=0||HIST!=0)return;
tetemp=TEND;
wipe_rep();
data_back();
if(new_float("Continue until:",&tetemp)==-1)return;
x=&MyData[0];
tetemp=fabs(tetemp);
if(fabs(MyTime)>=tetemp)return;
dif=tetemp-fabs(MyTime);
/* TEND=tetemp; */
MyStart=1; /* I know it is wasteful to restart, but lets be safe.... */
integrate(&MyTime,x,dif,DELTA_T,1,NJMP,&MyStart);
ping();
refresh_browser(storind);
}
range_item()
{
int i;
char bob[256];
i=find_user_name(PARAM,range.item);
if(i>-1){
range.type=PARAM;
range.index=i;
}
else {
i=find_user_name(IC,range.item);
if(i<=-1){
sprintf(bob," %s is not a parameter or variable !",range.item);
err_msg(bob);
return(0);
}
range.type=IC;
range.index=i;
}
return 1;
}
range_item2()
{
int i;
char bob[256];
i=find_user_name(PARAM,range.item2);
if(i>-1){
range.type2=PARAM;
range.index2=i;
}
else {
i=find_user_name(IC,range.item2);
if(i<=-1){
sprintf(bob," %s is not a parameter or variable !",range.item2);
err_msg(bob);
return(0);
}
range.type2=IC;
range.index2=i;
}
return 1;
}
set_up_range()
{
static char *n[]={"*3Range over","Steps","Start","End",
"Reset storage (Y/N)",
"Use old ic's (Y/N)","Cycle color (Y/N)","Movie(Y/N)"};
char values[8][MAX_LEN_SBOX];
int status,i;
static char *yn[]={"N","Y"};
if(!Xup){
return(range_item());
}
sprintf(values[0],"%s",range.item);
sprintf(values[1],"%d",range.steps);
sprintf(values[2],"%.16g",range.plow);
sprintf(values[3],"%.16g",range.phigh);
sprintf(values[4],"%s",yn[range.reset]);
sprintf(values[5],"%s",yn[range.oldic]);
sprintf(values[6],"%s",yn[range.cycle]);
sprintf(values[7],"%s",yn[range.movie]);
status=do_string_box(8,8,1,"Range Integrate",n,values,45);
if(status!=0){
strcpy(range.item,values[0]);
/* i=find_user_name(PARAM,range.item);
if(i>-1){
range.type=PARAM;
range.index=i;
}
else {
i=find_user_name(IC,range.item);
if(i<=-1){
err_msg("No such name!");
return(0);
}
range.type=IC;
range.index=i;
}
*/
if(range_item()==0)return 0;
range.steps=atoi(values[1]);
if(range.steps<=0)range.steps=10;
range.plow=atof(values[2]);
range.phigh=atof(values[3]);
if(values[4][0]=='Y'||values[4][0]=='y')range.reset=1;
else range.reset=0;
if(values[5][0]=='Y'||values[5][0]=='y')range.oldic=1;
else range.oldic=0;
if(values[6][0]=='Y'||values[6][0]=='y')range.cycle=1;
else range.cycle=0;
if(values[7][0]=='Y'||values[7][0]=='y')range.movie=1;
else range.movie=0;
/* printf("%s %d %d %d (%d %d) %f %f ",
range.item, range.steps,
range.reset,range.oldic,range.type,range.index,
range.plow,range.phigh);*/
RANGE_FLAG=1;
return(1);
}
return(0);
}
set_up_range2()
{
static char *n[]={"*3Vary1","Start1","End1",
"*3Vary2","Start2","End2","Steps",
"Reset storage (Y/N)",
"Use old ic's (Y/N)","Cycle color (Y/N)","Movie(Y/N)",
"Crv(1) Array(2)","Steps2"};
char values[13][MAX_LEN_SBOX];
int status,i;
static char *yn[]={"N","Y"};
if(!Xup){
return(range_item());
}
sprintf(values[0],"%s",range.item);
sprintf(values[1],"%.16g",range.plow);
sprintf(values[2],"%.16g",range.phigh);
sprintf(values[3],"%s",range.item2);
sprintf(values[4],"%.16g",range.plow2);
sprintf(values[5],"%.16g",range.phigh2);
sprintf(values[6],"%d",range.steps);
sprintf(values[7],"%s",yn[range.reset]);
sprintf(values[8],"%s",yn[range.oldic]);
sprintf(values[9],"%s",yn[range.cycle]);
sprintf(values[10],"%s",yn[range.movie]);
if(range.rtype==2)
sprintf(values[11],"2");
else
sprintf(values[11],"1");
sprintf(values[12],"%d",range.steps2);
status=do_string_box(13,7,2,"Double Range Integrate",n,values,45);
if(status!=0){
strcpy(range.item,values[0]);
if(range_item()==0)return 0;
strcpy(range.item2,values[3]);
if(range_item2()==0)return 0;
range.steps=atoi(values[6]);
range.steps2=atoi(values[12]);
if(range.steps<=0)range.steps=10;
if(range.steps2<=0)range.steps2=10;
range.plow=atof(values[1]);
range.phigh=atof(values[2]);
range.plow2=atof(values[4]);
range.phigh2=atof(values[5]);
if(values[7][0]=='Y'||values[7][0]=='y')range.reset=1;
else range.reset=0;
if(values[8][0]=='Y'||values[8][0]=='y')range.oldic=1;
else range.oldic=0;
if(values[9][0]=='Y'||values[9][0]=='y')range.cycle=1;
else range.cycle=0;
if(values[10][0]=='Y'||values[10][0]=='y')range.movie=1;
else range.movie=0;
range.rtype=atoi(values[11]);
RANGE_FLAG=1;
return(1);
}
return(0);
}
init_monte_carlo()
{
int i;
fixptguess.tol=.001;
fixptguess.n=100;
for(i=0;i<NODE;i++){
fixptguess.xlo[i]=-10;
fixptguess.xhi[i]=10;
}
fixptlist.flag=0;
fixptlist.n=0;
}
monte_carlo()
{
int append=0,sb=0;
int i=0,done=0;
double z;
char name[256];
new_int("Append(1/0",&append);
new_int("# Guesses:",&fixptguess.n);
new_float("Tolerance:",&fixptguess.tol);
while(1){
sprintf(name,"%s_lo :",uvar_names[i]);
z=fixptguess.xlo[i];
done=new_float(name,&z);
if(done==0)
fixptguess.xlo[i]=z;
if(done==-1)break;
sprintf(name,"%s_hi :",uvar_names[i]);
z=fixptguess.xhi[i];
done=new_float(name,&z);
if(done==0)
fixptguess.xhi[i]=z;
if(done==-1)break;
i++;
if(i>=NODE)
break;
}
do_monte_carlo_search(append, 1);
}
do_monte_carlo_search(int append, int stuffbrowse)
{
int i,j,k,m,n=fixptguess.n;
int ierr,new=1;
double x[MAXODE],sum;
double er[MAXODE],em[MAXODE];
if(append==0)
fixptlist.n=0;
if(fixptlist.flag==0){
for(i=0;i<MAXFP;i++){
fixptlist.x[i]=(double *)malloc(NODE*sizeof(double));
fixptlist.er[i]=(double *)malloc(NODE*sizeof(double));
fixptlist.em[i]=(double *)malloc(NODE*sizeof(double));
}
fixptlist.flag=1;
}
for(i=0;i<n;i++){
/* printf("Guess:\n"); */
for(j=0;j<NODE;j++){
x[j]=ndrand48()*(fixptguess.xhi[j]-fixptguess.xlo[j])+fixptguess.xlo[j];
/* printf("x[%d]=%g \n",j,x[j]); */
}
silent_fixpt(x,NEWT_ERR,EVEC_ERR,BOUND,EVEC_ITER,NODE,er,em,&ierr);
if(ierr==0){
m=fixptlist.n;
if(m==0){ /* first fixed point found */
fixptlist.n=1;
printf("Found: %d\n",m);
for(j=0;j<NODE;j++){
fixptlist.x[0][j]=x[j];
fixptlist.er[0][j]=er[j];
fixptlist.em[0][j]=em[j];
printf(" x[%d]= %g eval= %g + I %g \n",j,x[j],er[j],em[j]);
}
}
else { /* there are others better compare them */
new=1;
for(k=0;k<m;k++){
sum=0.0;
for(j=0;j<NODE;j++)
sum+=fabs(x[j]-fixptlist.x[k][j]);
if(sum<fixptguess.tol)
new=0;
}
if(new==1){
m=fixptlist.n;
fixptlist.n++;
if(m<MAXFP){
printf("Found: %d\n",m);
for(j=0;j<NODE;j++){
fixptlist.x[m][j]=x[j];
fixptlist.er[m][j]=er[j];
fixptlist.em[m][j]=em[j];
printf(" x[%d]= %g eval= %g + I %g \n",j,x[j],er[j],em[j]);
}
}
}
}
}
}
if(stuffbrowse) {
reset_browser();
storind=0;
m=fixptlist.n;
for(i=0;i<m;i++){
storage[0][storind]=(float)i;
for(j=0;j<NODE;j++)storage[j+1][storind]=(float)fixptlist.x[i][j];
storind++;
}
refresh_browser(storind);
}
}
do_eq_range(x)
double *x;
{
double parlo,parhi,dpar,temp;
int npar,stabcol,i,j,ierr;
char bob[50];
float stabinfo;
if(set_up_eq_range()==0)return;
wipe_rep();
data_back();
parlo=eq_range.plow;
parhi=eq_range.phigh;
npar=eq_range.steps;
dpar=(parhi-parlo)/(double)npar;
stabcol=eq_range.col;
storind=0;
DelayErr=0;
ENDSING=0;
PAR_FOL=1;
PAUSER=0;
SHOOT=eq_range.shoot;
reset_browser();
if(eq_range.movie)reset_film();
for(i=0;i<=npar;i++)
{
if(eq_range.movie)
clear_draw_window();
temp=parlo+dpar*(double)i;
set_val(eq_range.item,temp);
PAR_FOL=1;
sprintf(bob,"%s=%.16g",eq_range.item,temp);
bottom_msg(2,bob);
evaluate_derived();
if(DelayFlag)
do_delay_sing(x,NEWT_ERR,EVEC_ERR,BOUND,EVEC_ITER,
NODE,&ierr,&stabinfo);
else do_sing(x,NEWT_ERR,EVEC_ERR,BOUND,EVEC_ITER,
NODE,&ierr,&stabinfo);
if(eq_range.movie)
if(film_clip()==0)err_msg("Out of film");
storage[0][storind]=temp;
for(j=0;j<NODE;j++)storage[j+1][storind]=(float)x[j];
for(j=NODE;j<NODE+NMarkov;j++)storage[j+1][storind]=0.0;
if(stabcol>0)storage[stabcol-1][storind]=stabinfo;
storind++;
if(ENDSING==1)break;
}
refresh_browser(storind);
PAR_FOL=0;
}
swap_color(col,rorw)
int *col,rorw;
{
if(rorw)MyGraph->color[0]=*col;
else *col=MyGraph->color[0];
}
set_cycle(flag,icol)
int flag,*icol;
{
if(flag==0)return;
MyGraph->color[0]=*icol+1;
*icol=*icol+1;
if(*icol==10)*icol=0;
}
do_range(x,flag)
double *x;
int flag; /* 0 for 1-param 1 for 2 parameter */
{
char esc;
char bob[50];
int ivar=0,ivar2,res=0,oldic=0;
int nit=20,i,j,itype,itype2,cycle=0,icol=0,nit2,iii;
int color=MyGraph->color[0];
double t,dpar,plow=0.0,phigh=1.0,p,plow2,phigh2,p2,dpar2;
double temp,temp2;
int ierr=0;
if(flag==0){
range.rtype=0;
if(set_up_range()==0)return(-1);
}
else {
if(set_up_range2()==0)return -1;
}
MyStart=1;
itype=range.type;
ivar=range.index;
res=range.reset;
oldic=range.oldic;
nit=range.steps;
plow=range.plow;
phigh=range.phigh;
cycle=range.cycle;
dpar=(phigh-plow)/(double)nit;
get_ic(2,x);
storind=0;
STORFLAG=1;
PAUSER=0;
nit2=0;
if(range.rtype==2)nit2=range.steps2;
if(range.type==PARAM)get_val(range.item,&temp);
alloc_liap(nit); /* make space */
if(range.rtype>0){
itype2=range.type2;
ivar2=range.index2;
plow2=range.plow2;
phigh2=range.phigh2;
if(range.rtype==2)dpar2=(phigh2-plow2)/(double)nit2;
else dpar2=(phigh2-plow2)/(double)nit;
if(range.type2==PARAM)get_val(range.item2,&temp2);
}
if(range.movie)reset_film();
for(j=0;j<=nit2;j++){
for(i=0;i<=nit;i++)
{
if(range.movie)clear_draw_window();
if(cycle)MyGraph->color[0]=icol+1;
icol++;
if(icol==10)icol=0;
p=plow+dpar*(double)i;
if(range.rtype==1)
p2=plow2+dpar2*(double)i;
if(range.rtype==2)
p2=plow2+dpar2*(double)j;
if(oldic==1){
get_ic(1,x);
if(DelayFlag){
/* restart initial data */
if(do_init_delay(DELAY)==0)break;
}
}
t=T0;
MyStart=1;
POIEXT=0;
if(itype==IC)x[ivar]=p;
else {
set_val(range.item,p);
redo_all_fun_tables();
re_evaluate_kernels();
}
if(range.rtype>0){
if(itype2==IC)x[ivar2]=p2;
else {
set_val(range.item2,p2);
redo_all_fun_tables();
re_evaluate_kernels();
}
}
if(Xup){
if(range.rtype>0)
sprintf(bob,"%s=%.16g %s=%.16g",range.item,p,range.item2,p2);
else
sprintf(bob,"%s=%.16g i=%d",range.item,p,i);
bottom_msg(2,bob);
}
do_start_flags(x,&MyTime);
if(fabs(MyTime)>=TRANS&&STORFLAG==1&&POIMAP==0)
{
storage[0][0]=(float)MyTime;
extra(x,MyTime,NODE,NEQ);
for(iii=0;iii<NEQ;iii++)storage[1+iii][0]=(float)x[iii];
storind=1;
}
if(integrate(&t,x,TEND,DELTA_T,1,NJMP,&MyStart)==1){
ierr=-1;
break;
}
/* printf("storind = %d \n",storind); */
if(STOCH_FLAG)
append_stoch(i,storind);
if(range.movie)
if(film_clip()==0){err_msg("Out of film");break;}
refresh_browser(storind);
do_this_liaprun(i,p); /* sends parameter and index back */
if(storind>2)auto_freeze_it();
if(res==1||STOCH_FLAG)
{
if(batch_range==1)
write_this_run(batchout,i);
storind=0;
}
}
}
if(oldic==1)get_ic(1,x);
else get_ic(0,x);
if(range.type==PARAM)set_val(range.item,temp);
if(range.rtype>0)
if(range.type2==PARAM)set_val(range.item2,temp2);
evaluate_derived();
MyGraph->color[0]=color;
INFLAG=1;
/* refresh_browser(storind); */
ping();
if(STOCH_FLAG)
do_stats(nit,ierr);
return(ierr);
}
find_equilib_com(int com)
{
int i,ierr;
float xm,ym;
int im,jm;
int iv,jv;
float stabinfo;
double *x,oldtrans;
x=&MyData[0];
if(FFT||HIST||NKernel>0)return;
STORFLAG=0;
POIMAP=0;
oldtrans=TRANS;
TRANS=0.0;
evaluate_derived();
switch(com){
case 2:
do_eq_range(x);
return;
case 1:
/* Get mouse values */
iv=MyGraph->xv[0]-1;
jv=MyGraph->yv[0]-1;
if(iv<0||iv>=NODE||jv<0||jv>=NODE||MyGraph->grtype>=5||jv==iv){
err_msg("Not in useable 2D plane...");
return;
}
/* get mouse click x,y */
get_ic(1,x);
MessageBox("Click on guess");
if(GetMouseXY(&im,&jm)){
scale_to_real(im,jm,&xm,&ym);
x[iv]=(double)xm;
x[jv]=(double)ym;
}
KillMessageBox();
break;
case 3: monte_carlo();
return;
case 0:
default:
get_ic(2,x);
break;
}
if(DelayFlag){
do_delay_sing(x,NEWT_ERR,EVEC_ERR,BOUND,EVEC_ITER,NODE,&ierr,&stabinfo);
ping();
}
else
do_sing(x,NEWT_ERR,EVEC_ERR,BOUND,EVEC_ITER,NODE,&ierr,&stabinfo);
TRANS=oldtrans;
}
batch_integrate()
{
int i;
if(Nintern_set==0){
batch_integrate_once();
return;
}
for(i=0;i<Nintern_set;i++){
sprintf(batchout,"%s.dat",intern_set[i].name);
printf("out=%s\n",batchout);
extract_internset(i);
chk_delay();
printf(" Ok integrating now \n");
batch_integrate_once();
}
}
batch_integrate_once()
{
FILE *fp;
double *x;
int oldstart=MyStart,i;
MyStart=1;
x=&MyData[0];
RANGE_FLAG=0;
DelayErr=0;
MyTime=T0;
STORFLAG=1;
POIEXT=0;
storind=0;
reset_browser();
/* printf("batch_range=%d\n",batch_range); */
if(batch_range==1||STOCH_FLAG>0){
reset_dae();
RANGE_FLAG=1;
if(do_range(x,0)!=0)
printf(" Errors occured in range integration \n");
}
else {
get_ic(2,x);
if(DelayFlag){
/* restart initial data */
if(do_init_delay(DELAY)==0)return;
}
do_start_flags(x,&MyTime);
if(fabs(MyTime)>=TRANS&&STORFLAG==1&&POIMAP==0)
{
storage[0][0]=(float)MyTime;
extra(x,MyTime,NODE,NEQ);
for(i=0;i<NEQ;i++)storage[1+i][0]=(float)x[i];
storind=1;
}
if(integrate(&MyTime,x,TEND,DELTA_T,1,NJMP,&MyStart)!=0)
printf(" Integration not completed -- will write anyway...\n");
INFLAG=1;
refresh_browser(storind);
}
if(!batch_range || range.reset==0){
if(STOCH_FLAG==1)mean_back();
if(STOCH_FLAG==2)variance_back();
fp=fopen(batchout,"w");
if(fp==NULL){
printf(" Unable to open %s to write \n",batchout);
return;
}
write_mybrowser_data(fp);
fclose(fp);
}
printf(" Run complete ... \n");
}
write_this_run(file,i)
char *file;
int i;
{
char outfile[256];
FILE *fp;
sprintf(outfile,"%s.%d",file,i);
fp=fopen(outfile,"w");
if(fp==NULL){
printf("Couldnt open %s\n",outfile);
return -1;
}
write_mybrowser_data(fp);
fclose(fp);
}
do_init_data(int com)
{
char sr[20],ch;
int i,si;
double *x;
double old_dt=DELTA_T;
FILE *fp;
char icfile[256];
float xm,ym;
int im,jm,oldstart,iv,jv,badmouse;
oldstart=MyStart;
MyStart=1;
x=&MyData[0];
RANGE_FLAG=0;
DelayErr=0;
reset_dae();
if(FFT||HIST)return;
if(com==M_ID){ /* dont want to wipe out everything! */
get_new_guesses();
return;
}
data_back();
wipe_rep();
MyTime=T0;
STORFLAG=1;
POIEXT=0;
storind=0;
reset_browser();
switch(com){
case M_IR: /* do range */
do_range(x,0);
return;
case M_I2:
do_range(x,1);
return;
case M_IS:
case M_IL:
if(INFLAG==0){
ping();
err_msg("No prior solution");
return;
}
get_ic(0,x);
if(com==M_IS){
T0=LastTime;
MyTime=T0;
}
if(METHOD==VOLTERRA&&oldstart==0){
ch=(char)TwoChoice("No","Yes","Reset integrals?","ny");
if(ch=='n')MyStart=oldstart;
}
break;
case M_IO:
get_ic(1,x);
if(DelayFlag){
/* restart initial data */
if(do_init_delay(DELAY)==0)return;
}
set_init_guess();
break;
case M_IM:
case M_II:
iv=MyGraph->xv[0]-1;
jv=MyGraph->yv[0]-1;
if(iv<0||iv>=NODE||jv<0||jv>=NODE||MyGraph->grtype>=5||jv==iv){
err_msg("Not in useable 2D plane...");
return;
}
/* Get mouse values */
if(com==M_IM){
get_ic(1,x);
MessageBox("Click on initial data");
if(GetMouseXY(&im,&jm)){
scale_to_real(im,jm,&xm,&ym);
im=MyGraph->xv[0]-1;
jm=MyGraph->yv[0]-1;
x[iv]=(double)xm;
x[jv]=(double)ym;
last_ic[im]=x[im];
last_ic[jm]=x[jm];
KillMessageBox();
if(DelayFlag){
/* restart initial data */
if(do_init_delay(DELAY)==0)return;
}
}
else {
KillMessageBox();
return;
}
}
else {
SuppressBounds=1;
MessageBox("Click on initial data -- ESC to quit");
while(1){
get_ic(1,x);
badmouse=GetMouseXY(&im,&jm);
if(badmouse==0)break;
scale_to_real(im,jm,&xm,&ym);
im=MyGraph->xv[0]-1;
jm=MyGraph->yv[0]-1;
x[iv]=(double)xm;
x[jv]=(double)ym;
last_ic[im]=x[im];
last_ic[jm]=x[jm];
if(DelayFlag){
/* restart initial data */
if(do_init_delay(DELAY)==0)break;
}
MyStart=1;
MyTime=T0;
usual_integrate_stuff(x);
}
KillMessageBox();
SuppressBounds=0;
return;
}
break;
case M_IN:
man_ic();
get_ic(2,x);
set_init_guess();
break;
case M_IU:
if(form_ic()==0)return;
get_ic(2,x);
break;
case M_IH:
if(ShootICFlag==0){
err_msg("No shooting data available");
break;
}
sprintf(sr,"Which? (1-%d)",ShootIndex);
si=1;
new_int(sr,&si);
si--;
if(si<ShootIndex&&si>=0){
for(i=0;i<NODE;i++)
last_ic[i]=ShootIC[si][i];
get_ic(2,x);
}
else
err_msg("Out of range");
break;
case M_IF:
icfile[0]=0;
if(!file_selector("Read initial data",icfile,"*.dat"))return;
/* if(new_string("Filename: ",icfile)==0)return; */
if((fp=fopen(icfile,"r"))==NULL){
err_msg(" Cant open IC file");
return;
}
for(i=0;i<NODE;i++)
fscanf(fp,"%lg",&last_ic[i]);
fclose(fp);
get_ic(2,x);
break;
case M_IB:
DELTA_T=-fabs(DELTA_T);
get_ic(2,x);
set_init_guess();
if(DelayFlag){
/* restart initial data */
if(do_init_delay(DELAY)==0)return;
}
break;
case M_IG:
default:
set_init_guess();
get_ic(2,x);
if(DelayFlag){
/* restart initial data */
if(do_init_delay(DELAY)==0)return;
}
break;
}
usual_integrate_stuff(x);
DELTA_T=old_dt;
}
run_now()
{
int i,si;
double *x;
MyStart=1;
x=&MyData[0];
RANGE_FLAG=0;
DelayErr=0;
reset_dae();
MyTime=T0;
get_ic(2,x);
STORFLAG=1;
POIEXT=0;
storind=0;
reset_browser();
usual_integrate_stuff(x);
}
do_start_flags(double *x,double *t)
{
int iflagstart=1;
double tnew=*t;
double sss;
one_flag_step(x,x,&iflagstart,*t,&tnew,NODE,&sss);
}
usual_integrate_stuff(x)
double *x;
{
int i;
do_start_flags(x,&MyTime);
if(fabs(MyTime)>=TRANS&&STORFLAG==1&&POIMAP==0)
{
storage[0][0]=(float)MyTime;
extra(x,MyTime,NODE,NEQ);
for(i=0;i<NEQ;i++)storage[1+i][0]=(float)x[i];
storind=1;
}
integrate(&MyTime,x,TEND,DELTA_T,1,NJMP,&MyStart);
ping();
INFLAG=1;
refresh_browser(storind);
auto_freeze_it();
redraw_ics();
}
/* form_ic -- u_i(0) = F(i) where "i" is represented by "t"
or
u[5..20]=f([j])
*/
do_new_array_ic(new,j1,j2)
char *new;
int j1,j2;
{
int i;
int ihot=-1;
int ifree=-1;
/* first check to see if this is
one that has already been used and also find the first free one
*/
for(i=0;i<NAR_IC;i++){
if(ar_ic[i].index0==-1&&ifree==-1&&ar_ic[i].type==0)
ifree=i;
if(strcmp(ar_ic[i].var,new)==0&&ar_ic[i].j1==j1&&ar_ic[i].j2==j2)
ihot=i;
}
if(ihot==-1){
if(ifree==-1){
ihot=0;
}
else {
ihot=ifree;
}
/* copy relevant stuff */
strcpy(ar_ic[ihot].var,new);
ar_ic[ihot].type=2;
ar_ic[ihot].j1=j1;
ar_ic[ihot].j2=j2;
}
new_string("Formula:",ar_ic[ihot].formula);
/* now we have everything we need */
evaluate_ar_ic(ar_ic[ihot].var,ar_ic[ihot].formula,
ar_ic[ihot].j1,ar_ic[ihot].j2);
}
store_new_array_ic(new,j1,j2,formula)
char *new,*formula;
int j1,j2;
{
int i;
int ihot=-1;
int ifree=-1;
/* first check to see if this is
one that has already been used and also find the first free one
*/
for(i=0;i<NAR_IC;i++){
if(ar_ic[i].index0==-1&&ifree==-1&&ar_ic[i].type==0)
ifree=i;
if(strcmp(ar_ic[i].var,new)==0&&ar_ic[i].j1==j1&&ar_ic[i].j2==j2)
ihot=i;
}
if(ihot==-1){
if(ifree==-1){
ihot=0;
}
else {
ihot=ifree;
}
/* copy relevant stuff */
strcpy(ar_ic[ihot].var,new);
ar_ic[ihot].type=2;
ar_ic[ihot].j1=j1;
ar_ic[ihot].j2=j2;
}
strcpy(ar_ic[ihot].formula,formula);
}
evaluate_ar_ic(v,f,j1,j2)
char *v,*f;
int j1,j2;
{
int j;
int i,flag;
double z;
char vp[25],fp[256];
for(j=j1;j<=j2;j++){
i=-1;
subsk(v,vp,j,1);
find_variable(vp,&i);
if(i>0){
subsk(f,fp,j,1);
flag=do_calc(fp,&z);
if(flag!=-1)
last_ic[i-1]=z;
else
return;
}
}
}
extract_ic_data(char *big)
{
int i,n,j;
int j1,j2,flag2;
char front[40],new[50],c;
char back[256];
de_space(big);
i=0;
n=strlen(big);
while(1){
c=big[i];
if(c=='(')break;
front[i]=c;
i++;
if(i>=n){
return(-1);
}
}
front[i]=0;
/* lets find the back part */
i=i+4;
for(j=i;j<n;j++)
back[j-i]=big[j];
back[j-i]=0;
/* now fix it up */
big[0]='#';
big[1]=' ';
search_array(front,new,&j1,&j2,&flag2);
if(flag2==1){
store_new_array_ic(new,j1,j2,back);
ar_ic_defined=1;
}
}
arr_ic_start()
{
int i;
if(ar_ic_defined==0) return;
for(i=0;i<NAR_IC;i++){
if(ar_ic[i].type==2){
evaluate_ar_ic(ar_ic[i].var,ar_ic[i].formula,
ar_ic[i].j1,ar_ic[i].j2);
}
}
}
set_array_ic()
{
char junk[50];
char new[50];
int i,index0,myar=-1;
int i1,in;
int j1,j2,flag2;
double z;
int flag;
junk[0]=0;
if(new_string("Variable: ",junk)==0)return 0;
search_array(junk,new,&j1,&j2,&flag2);
if(flag2==1)
{
do_new_array_ic(new,j1,j2);
}
else {
find_variable(junk,&i);
if(i<=-1)
return 0;
index0=i;
for(i=0;i<NAR_IC;i++){
if(ar_ic[i].type==2)continue;
if(index0==ar_ic[i].index0){
myar=i;
break;
}
}
if(myar<0){
for(i=0;i<NAR_IC;i++){
if(ar_ic[i].type==2)continue;
if(ar_ic[i].index0==-1){
myar=i;
break;
}
}
}
if(myar<0)myar=0;
/* Now we have an element in the array index */
ar_ic[myar].index0=index0;
ar_ic[myar].type=0;
new_int("Number elements:",&ar_ic[myar].n);
new_string("u=F(t-i0):",ar_ic[myar].formula);
i1=index0-1;
in=i1+ar_ic[myar].n;
/* printf("i1=%d in=%d \n",i1,in); */
if(i1>NODE||in>NODE)return 0; /* out of bounds */
for(i=i1;i<in;i++){
set_val("t",(double)(i-i1));
flag=do_calc(ar_ic[myar].formula,&z);
if(flag==-1){
err_msg("Bad formula");
return 1;
}
last_ic[i]=z;
}
}
return 1;
}
form_ic()
{
int ans;
while(1){
ans=set_array_ic();
if(ans==0)break;
}
return 1;
}
get_ic(it,x)
int it;
double *x;
{
int i;
switch(it){
case 0:
for(i=0;i<NODE+NMarkov;i++)last_ic[i]=x[i];
break;
case 1:
case 2:
for(i=0;i<NODE+NMarkov;i++)x[i]=last_ic[i];
break;
}
}
ode_int(y,t,istart,ishow)
double *y,*t;
int *istart,ishow;
{
double error[MAXODE];
int kflag,i;
int nodes=xpv.node+xpv.nvec;
int nit,nout=NJMP;
double tend=TEND;
double dt=DELTA_T,tout;
if(METHOD==0){
nit=tend;
dt=dt/fabs(dt);
}
else nit=(tend+.1*fabs(dt))/fabs(dt);
if(ishow==1){
integrate(t,y,tend,dt,1,nout,istart);
return(1);
}
MSWTCH(xpv.x,y);
evaluate_derived();
if(METHOD<GEAR ||METHOD==BACKEUL){
kflag=solver(xpv.x,t,dt,nit,nodes,istart,WORK);
MSWTCH(y,xpv.x);
if(kflag<0){
ping();
if(RANGE_FLAG)return(0);
switch(kflag)
{
case -1: err_msg(" Singular Jacobian "); break;
case -2: err_msg("Too many iterates");break;
}
return(0);
}
}
else
{
tout=*t+tend*dt/fabs(dt);
switch(METHOD){
case GEAR:
if(*istart==1)*istart=0;
gear(nodes,t,tout,xpv.x,HMIN,HMAX,TOLER,2,error,
&kflag,istart,WORK,IWORK);
MSWTCH(y,xpv.x);
if(kflag<0)
{
ping();
if(RANGE_FLAG)return(0);
switch(kflag)
{
case -1: err_msg("kflag=-1: minimum step too big"); break;
case -2: err_msg("kflag=-2: required order too big");break;
case -3: err_msg("kflag=-3: minimum step too big");break;
case -4: err_msg("kflag=-4: tolerance too small");break;
}
return(0);
}
break;
#ifdef CVODE_YES
case CVODE:
cvode(istart,xpv.x,t,nodes,tout,&kflag,&TOLER,&ATOLER);
MSWTCH(y,xpv.x);
if(kflag<0){
cvode_err_msg(kflag);
return(0);
}
end_cv();
break;
#endif
case DP5:
case DP83:
dp(istart,xpv.x,t,nodes,tout,&TOLER,&ATOLER,METHOD-DP5,&kflag);
MSWTCH(y,xpv.x);
if(kflag<0){
if(RANGE_FLAG)return(0);
dp_err(kflag);
return 0;
}
break;
case RB23:
rb23(xpv.x,t,tout,istart,nodes,WORK,&kflag);
MSWTCH(y,xpv.x);
if(kflag<0){
ping();
if(RANGE_FLAG)return(0);
err_msg("Step size too small");
return 0;
}
break;
case RKQS:
case STIFF:
adaptive(xpv.x,nodes,t,tout,TOLER,&dt,
HMIN,WORK,&kflag,NEWT_ERR,METHOD,istart);
MSWTCH(y,xpv.x);
if(kflag){
ping();
if(RANGE_FLAG)return(0);
switch(kflag){
case 2: err_msg("Step size too small"); break;
case 3: err_msg("Too many steps"); break;
case -1: err_msg("singular jacobian encountered"); break;
case 1: err_msg("stepsize is close to 0"); break;
case 4: err_msg("exceeded MAXTRY in stiff"); break;
}
return(0);
}
break;
}
}
return(1);
}
integrate(t,x,tend,dt, count, nout,start)
double *t, *x, tend,dt;
int count,nout, *start;
{
float xv[MAXODE+1],xvold[MAXODE+1];
float oldperiod=0.0;
double error[MAXODE];
double xprime[MAXODE],oldxprime[MAXODE],hguess=dt;
int kflag;
int torcross[MAXODE];
int nodes=xpv.node+xpv.nvec-NMarkov;
int perflag=0;
int rval=0;
double oldx[MAXODE],oldt,dint,dxp,sect,sect1,tout,tzero=*t;
double sss,tnew=*t;
int iflagstart=1;
float tscal=tend,tstart=*t,tv;
float dp_time,xp_time;
char esc;
char error_message[50];
int ieqn,koff,i,pflag=0;
int icount=0;
int nit;
int cwidth;
/* new poincare map stuff */
double pmapf,pmapfold;
int i_nan=0; /* NaN */
MSWTCH(xpv.x,x);
if(Xup) cwidth=get_command_width();
LastTime=*t;
evaluate_derived();
if((METHOD==GEAR)&&(*start==1))*start=0;
if(METHOD==0){
nit=tend;
dt=dt/fabs(dt);
}
else nit=(tend+fabs(dt)*.1)/fabs(dt);
/* else nit=tend/fabs(dt); */
nit=(nit+nout-1)/nout;
if(nit==0)return(rval);
one_flag_step(xpv.x,xpv.x,&iflagstart,*t,&tnew,nodes,&sss);
MSWTCH(x,xpv.x);
extra(x,*t,NODE,NEQ); /* Note this takes care of initializing Markov variables */
MSWTCH(xpv.x,x);
xv[0]=(float)*t;
for(ieqn=1;ieqn<=NEQ;ieqn++)xv[ieqn]=(float)x[ieqn-1];
if(animation_on_the_fly)on_the_fly(1);
/* if(POIMAP==4)
pmapfold=get_map_value(x,*t); */
if(POIMAP)
{
oldt=*t;
for(ieqn=0;ieqn<NEQ;ieqn++)oldx[ieqn]=x[ieqn];
}
if(dt<0.0)tscal=-tend;
if(tscal==0.0)tscal=1.0;
stor_delay(x);
while(1)
{
switch(METHOD){
case GEAR:
{
tout=tzero+dt*(icount+1);
if(fabs(dt)<fabs(HMIN)){
LastTime=*t;
return(1);
}
MSWTCH(xpv.x,x);
gear(nodes,t,tout,xpv.x,HMIN,HMAX,TOLER,2,error,&kflag,start,WORK,IWORK);
MSWTCH(x,xpv.x);
stor_delay(x);
if(DelayErr){
DelayErr=0;
LastTime=*t;
err_dae();
return(1);
}
if(kflag<0)
{
ping();
if(RANGE_FLAG||SuppressBounds){
LastTime=*t;
return(1);
}
switch(kflag)
{
case -1: err_msg("kflag=-1: minimum step too big"); break;
case -2: err_msg("kflag=-2: required order too big");break;
case -3: err_msg("kflag=-3: minimum step too big");break;
case -4: err_msg("kflag=-4: tolerance too small");break;
}
LastTime=*t;
return(1);
}
}
break;
#ifdef CVODE_YES
case CVODE:
tout=tzero+dt*(icount+1);
if(fabs(dt)<fabs(HMIN)){
LastTime=*t;
end_cv();
return(1);
}
MSWTCH(xpv.x,x);
cvode(start,xpv.x,t,nodes,tout,&kflag,&TOLER,&ATOLER);
MSWTCH(x,xpv.x);
stor_delay(x);
if(DelayErr){
DelayErr=0;
err_dae();
LastTime=*t;
return(1);
}
if(kflag<0){
ping();
if(RANGE_FLAG||SuppressBounds){
LastTime=*t;
return(1);
}
cvode_err_msg(kflag);
LastTime=*t;
return(1);
}
break;
#endif
case DP5:
case DP83:
tout=tzero+dt*(icount+1);
if(fabs(dt)<fabs(HMIN)){
LastTime=*t;
return(1);
}
MSWTCH(xpv.x,x);
dp(start,xpv.x,t,nodes,tout,&TOLER,&ATOLER,METHOD-DP5,&kflag);
MSWTCH(x,xpv.x);
stor_delay(x);
if(DelayErr){
DelayErr=0;
err_dae();
LastTime=*t;
return(1);
}
if(kflag<0){
if(RANGE_FLAG||SuppressBounds){
LastTime=*t;
return(1);
}
dp_err(kflag);
LastTime=*t;
return(1);
}
break;
case RB23:
tout=tzero+dt*(icount+1);
if(fabs(dt)<fabs(HMIN)){
LastTime=*t;
return(1);
}
MSWTCH(xpv.x,x);
rb23(xpv.x,t,tout,start,nodes,WORK,&kflag);
MSWTCH(x,xpv.x);
stor_delay(x);
if(DelayErr){
DelayErr=0;
err_dae();
LastTime=*t;
return(1);
}
if(kflag<0){
if(RANGE_FLAG||SuppressBounds){
LastTime=*t;
return(1);
}
err_msg("Step size too small");
LastTime=*t;
return(1);
}
break;
case RKQS:
case STIFF:
tout=tzero+dt*(icount+1);
if(fabs(dt)<fabs(HMIN)){
LastTime=*t;
return(1);
}
MSWTCH(xpv.x,x);
adaptive(xpv.x,nodes,t,tout,TOLER,&hguess,
HMIN,WORK,&kflag,NEWT_ERR,METHOD,start);
MSWTCH(x,xpv.x);
stor_delay(x);
if(DelayErr){
DelayErr=0;
err_dae();
LastTime=*t;
return(1);
}
if(kflag){
ping();
if(RANGE_FLAG||SuppressBounds){
LastTime=*t;
return(1);
}
switch(kflag){
case 2: err_msg("Step size too small"); break;
case 3: err_msg("Too many steps"); break;
case -1: err_msg("singular jacobian encountered"); break;
case 1: err_msg("stepsize is close to 0"); break;
case 4: err_msg("exceeded MAXTRY in stiff"); break;
}
LastTime=*t;
return(1);
}
break;
default: {
MSWTCH(xpv.x,x);
kflag=solver(xpv.x,t,dt,nout,nodes,start,WORK);
MSWTCH(x,xpv.x);
if(kflag<0)
{
ping();
if(RANGE_FLAG||SuppressBounds)break;
switch(kflag)
{
case -1: err_msg("Singular Jacobian "); break;
case -2: err_msg("Too many iterates ");break;
}
LastTime=*t;
return(1);
}
}
}
/* START POST INTEGRATE STUFF */
extra(x,*t,NODE,NEQ);
if (TORUS == 1) {
for (ieqn = 0; ieqn < NEQ; ieqn++) {
torcross[ieqn]=0;
if (itor[ieqn] == 1) {
if (x[ieqn] > TOR_PERIOD) {
x[ieqn] -= TOR_PERIOD;
torcross[ieqn]=-1;
/* for (ip = 0; ip < NPlots; ip++) {
if (ieqn == IYPLT[ip]-1)
oldypl[ip] -= (float) TOR_PERIOD;
if (ieqn == IXPLT[ip]-1)
oldxpl[ip] -= (float) TOR_PERIOD;
if ((ieqn == IZPLT[ip]-1) && (PLOT_3D == 1))
oldzpl[ip] -= (float) TOR_PERIOD;
} */
}
if (x[ieqn] < 0) {
x[ieqn] += TOR_PERIOD;
torcross[ieqn]=1;
/* for (ip = 0; ip < NPlots; ip++) {
if (ieqn == IYPLT[ip]-1)
oldypl[ip] += (float) TOR_PERIOD;
if (ieqn == IXPLT[ip]-1)
oldxpl[ip] += (float) TOR_PERIOD;
if ((ieqn == IZPLT[ip]-1) &&
(MyGraph->ThreeDFlag == 1))
oldzpl[ip] += (float) TOR_PERIOD;
} */
}
}
}
}
xvold[0]=xv[0];
for(ieqn=1;ieqn<(NEQ+1);ieqn++)
{
xvold[ieqn]=xv[ieqn];
xv[ieqn]=(float)x[ieqn-1];
/* trap NaN using isnan() in math.h
modified the out of bounds message as well
print all the variables on the terminal window, haven't decide
should I store them or not.
If use with nout=1, can pinpoint the offensive variable(s)
*/
if(isnan(x[ieqn-1])!=0)
{
sprintf(error_message," %s is NaN at t = %f ",
uvar_names[ieqn-1],*t);
/* if((STORFLAG==1)&&(storind<MAXSTOR))
{ */ i_nan=0;
fprintf(stderr,"variable\tf(t-1)\tf(t) \n");
/* storage[i_nan][storind]=*t; */
for(i_nan=1;i_nan<=ieqn;i_nan++)
{/*storage[i_nan][storind]=xv[i_nan];*/
fprintf(stderr," %s\t%g\t%g\n",
uvar_names[i_nan-1],xvold[i_nan],xv[i_nan]);
}
for(;i_nan<=NEQ;i_nan++)
{
/*storage[i_nan][storind]=(float)x[i_nan-1];*/
fprintf(stderr," %s\t%g\t%g\n",
uvar_names[i_nan-1],xv[i_nan],(float)x[i_nan-1]);
}
/* storind++;
if(!(storind<MAXSTOR))
if(stor_full()==0)break;
}
*/
err_msg(error_message);
rval=1;
break;
}
/* end of NaN */
if(fabs(x[ieqn-1])>BOUND)
{
if(RANGE_FLAG||SuppressBounds)break;
sprintf(error_message," %s out of bounds at t = %f ",
uvar_names[ieqn-1],*t);
/* if((STORFLAG==1)&&(storind<MAXSTOR))
{ */ i_nan=0;
fprintf(stderr,"variable\tf(t-1)\tf(t) \n");
/* storage[i_nan][storind]=*t; */
for(i_nan=1;i_nan<=ieqn;i_nan++)
{/*storage[i_nan][storind]=xv[i_nan];*/
fprintf(stderr," %s\t%g\t%g\n",
uvar_names[i_nan-1],xvold[i_nan],xv[i_nan]);
}
for(;i_nan<=NEQ;i_nan++)
{
/*storage[i_nan][storind]=(float)x[i_nan-1];*/
fprintf(stderr," %s\t%g\t%g\n",
uvar_names[i_nan-1],xv[i_nan],(float)x[i_nan-1]);
}
/* storind++;
if(!(storind<MAXSTOR))
if(stor_full()==0)break;
}
*/
err_msg(error_message);
rval=1;
break;
}
}
/* This is where the progresser goes */
if(Xup){ plot_command(nit,icount,cwidth);
esc=my_abort();
{
if(esc==ESCAPE) break;
if(esc=='/'){rval=1;ENDSING=1;break;}
}
}
if(STOP_FLAG==1){STOP_FLAG=0;break;}
if(DelayErr){err_dae();rval=1;ENDSING=1;DelayErr=0;break;}
/* printf(" NEQ=%d ieqn = %d \n",NEQ,ieqn); */
if(ieqn<(NEQ+1))break;
tv=(float)*t;
xv[0]=tv;
if((POIMAP==2)&&!(POIVAR==0))
{
pflag=0;
if((oldx[POIVAR-1]<x[POIVAR-1])&&!(POIEXT<0))POIEXT=1;
if((oldx[POIVAR-1]>x[POIVAR-1])&&!(POIEXT>0))POIEXT=-1;
if( ( !(oldx[POIVAR-1]<x[POIVAR-1]) && (POIEXT>0) )||
( !(oldx[POIVAR-1]>x[POIVAR-1]) && (POIEXT<0) )
)
{
if(POISGN*POIEXT>=0)
{
/* We will interpolate to get a good local extremum */
rhs(*t,x,xprime,NEQ);
rhs(oldt,oldx,oldxprime,NEQ);
dxp=xprime[POIVAR-1]-oldxprime[POIVAR-1];
if(dxp==0.0){
err_msg("Cannot zero RHS for max/min - use a variable");
return(1);
}
dint=xprime[POIVAR-1]/dxp;
tv=(1-dint)**t+dint*oldt;
xv[0]=tv;
for(i=1;i<=NEQ;i++)xv[i]=dint*oldx[i-1]+(1-dint)*x[i-1];
pflag=1;
}
POIEXT=-POIEXT;
}
goto poi;
}
/* here is code for a formula type map -- F(X,t)=0
*/
if(POIMAP==4) {
/* pmapf=get_map_value(x,*t);
*/
}
if(POIMAP==1||POIMAP==3)
{
if(POIVAR==0)
{
sect1=fmod(fabs(oldt),fabs(POIPLN));
sect=fmod(fabs(*t),fabs(POIPLN));
if(sect<sect1)
{
dint=sect/(POIPLN+sect-sect1);
i=(int)(fabs(*t)/fabs(POIPLN));
tv=(float)POIPLN*i;
xv[0]=tv;
for(i=1;i<=NEQ;i++)xv[i]=(float)(dint*oldx[i-1]+(1-dint)*x[i-1]);
pflag=1;
}
else pflag=0;
}
else
{
if(!(POISGN<0))
{
if((oldx[POIVAR-1]<POIPLN)&&!(x[POIVAR-1]<POIPLN))
{
dint=(x[POIVAR-1]-POIPLN)/(x[POIVAR-1]-oldx[POIVAR-1]);
tv=(1-dint)**t+dint*oldt;
xv[0]=tv;
for(i=1;i<=NEQ;i++)xv[i]=dint*oldx[i-1]+(1-dint)*x[i-1];
pflag=1;
goto poi;
}
else pflag=0;
}
if(!(POISGN>0))
{
if((oldx[POIVAR-1]>POIPLN)&&!(x[POIVAR-1]>POIPLN))
{
dint=(x[POIVAR-1]-POIPLN)/(x[POIVAR-1]-oldx[POIVAR-1]);
tv=(1-dint)**t+dint*oldt;
xv[0]=tv;
for(i=1;i<=NEQ;i++)xv[i]=dint*oldx[i-1]+(1-dint)*x[i-1];
pflag=1;
}
else pflag=0;
}
}
poi: for(i=0;i<NEQ;i++)oldx[i]=x[i];
oldt=*t;
if(pflag==0)goto out;
}
/* Plotting and storing data */
if(POIMAP==3&&pflag==1){
if(oldperiod==0.0){
pflag=0; /* this is the first hit !! */
oldperiod=*t;
goto out;
}
xv[0]=*t-oldperiod;
oldperiod=*t;
}
if(!(fabs(*t)<TRANS)&&Xup)
{
plot_the_graphs(xv,xvold,NODE,NEQ,fabs(dt*NJMP),torcross);
}
if((STORFLAG==1)&&(count!=0)&&(storind<MAXSTOR)&&!(fabs(*t)<TRANS))
{
if(animation_on_the_fly)on_the_fly(0);
for(ieqn=0;ieqn<=NEQ;ieqn++)
storage[ieqn][storind]=xv[ieqn];
storind++;
if(!(storind<MAXSTOR))
if(stor_full()==0)break;
if((pflag==1)&&(SOS==1))break;
}
out:
icount++;
if(icount>=nit&&count!=0)break;
/* END POST INTEGRATE ANALYSIS */
}
LastTime=*t;
#ifdef CVODE_YES
if(METHOD==CVODE)
end_cv();
#endif
return(rval);
}
send_halt(double *y, double t)
{
STOP_FLAG=1;
}
send_output(double *y,double t)
{
double yy[MAXODE];
int i;
for(i=0;i<NODE;i++)
yy[i]=y[i];
extra(yy,t,NODE,NEQ);
if((STORFLAG==1)&&(storind<MAXSTOR)){
for(i=0;i<NEQ;i++)
storage[i+1][storind]=(float)yy[i];
storage[0][storind]=(float)t;
storind++;
}
}
do_plot(oldxpl,oldypl, oldzpl,xpl, ypl, zpl)
float *oldxpl, *oldypl, *oldzpl,*xpl, *ypl,*zpl;
{
int ip,np=MyGraph->nvars;
for(ip=0;ip<np;ip++){
if(MyGraph->ColorFlag==0){
set_linestyle(MyGraph->color[ip]);
}
/* if(MyGraph->line[ip]<0)
continue; */
if(MyGraph->line[ip]<=0)
{
PointRadius=-MyGraph->line[ip];
if(MyGraph->ThreeDFlag==0) point_abs(xpl[ip],ypl[ip]);
else point_3d(xpl[ip],ypl[ip],zpl[ip]);
}
else
{
if(MyGraph->ThreeDFlag==0){
line_abs(oldxpl[ip],oldypl[ip],xpl[ip],ypl[ip]);
}
else line_3d(oldxpl[ip],oldypl[ip],oldzpl[ip],
xpl[ip],ypl[ip],zpl[ip]);
}
}
}
/*
old restore is in restore.c
*/
export_data(fp)
FILE *fp;
{
int ip,np=MyGraph->nvars;
int ZSHFT,YSHFT,XSHFT;
int i,j,kxoff,kyoff,kzoff;
int iiXPLT,iiYPLT,iiZPLT;
int strind=get_maxrow_browser();
int i1=0,i2=strind;
float **data;
data=get_browser_data();
XSHFT=MyGraph->xshft;
YSHFT=MyGraph->yshft;
ZSHFT=MyGraph->zshft;
if(i1<ZSHFT)i1=ZSHFT;
if(i1<YSHFT)i1=YSHFT;
if(i1<XSHFT)i1=XSHFT;
if(strind<2)return;
kxoff=i1-XSHFT;
kzoff=i1-ZSHFT;
kyoff=i1-YSHFT;
iiXPLT=MyGraph->xv[0];
iiYPLT=MyGraph->yv[0];
if(MyGraph->ThreeDFlag>0){
iiZPLT=MyGraph->zv[0];
for(j=i1;j<strind;j++){
fprintf(fp,"%g %g %g \n",
data[iiXPLT][kxoff],
data[iiYPLT][kyoff],
data[iiZPLT][kzoff]);
kxoff++;
kyoff++;
kzoff++;
}
return;
}
/* 2D graph so we will save y from each curve */
for(j=i1;j<strind;j++){
fprintf(fp,"%g ",data[iiXPLT][kxoff]);
for(ip=0;ip<np;ip++){
iiYPLT=MyGraph->yv[ip];
fprintf(fp,"%g ",data[iiYPLT][kyoff]);
}
fprintf(fp,"\n");
kxoff++;
kyoff++;
kzoff++;
}
return;
}
plot_the_graphs(float *xv,float *xvold,int node,int neq,double ddt,int *tc)
{
int i;
int ic=current_pop;
if(SimulPlotFlag==0){
plot_one_graph(xv,xvold,node,neq,ddt,tc);
return;
}
for(i=0;i<num_pops;i++){
make_active(ActiveWinList[i]);
plot_one_graph(xv,xvold,node,neq,ddt,tc);
}
make_active(ic);
}
plot_one_graph(float *xv,float *xvold,int node,int neq,double ddt,int *tc)
{
int *IXPLT,*IYPLT,*IZPLT;
int NPlots,ip;
float oldxpl[MAXPERPLOT],oldypl[MAXPERPLOT],oldzpl[MAXPERPLOT];
float xpl[MAXPERPLOT],ypl[MAXPERPLOT],zpl[MAXPERPLOT];
NPlots=MyGraph->nvars;
IXPLT=MyGraph->xv;
IYPLT=MyGraph->yv;
IZPLT=MyGraph->zv;
for(ip=0;ip<NEQ;ip++){
if(itor[ip]==1)
xvold[ip+1]=xvold[ip+1]+tc[ip]*TOR_PERIOD;
}
for(ip=0;ip<NPlots;ip++){
oldxpl[ip]=xvold[IXPLT[ip]];
oldypl[ip]=xvold[IYPLT[ip]];
oldzpl[ip]=xvold[IZPLT[ip]];
xpl[ip]=xv[IXPLT[ip]];
ypl[ip]=xv[IYPLT[ip]];
zpl[ip]=xv[IZPLT[ip]];
}
if(MyGraph->ColorFlag)
comp_color(xv,xvold,NODE,(float)ddt);
do_plot(oldxpl,oldypl,oldzpl,xpl,ypl,zpl);
}
restore(i1,i2)
int i1,i2;
{
int ip,np=MyGraph->nvars;
int ZSHFT,YSHFT,XSHFT;
int i,j,kxoff,kyoff,kzoff;
int iiXPLT,iiYPLT,iiZPLT;
float oldxpl,oldypl,oldzpl,xpl,ypl,zpl;
float v1[MAXODE+1],v2[MAXODE+1];
float **data;
data=get_browser_data();
XSHFT=MyGraph->xshft;
YSHFT=MyGraph->yshft;
ZSHFT=MyGraph->zshft;
if(i1<ZSHFT)i1=ZSHFT;
if(i1<YSHFT)i1=YSHFT;
if(i1<XSHFT)i1=XSHFT;
if(storind<2)return;
for(ip=0;ip<np;ip++){
kxoff=i1-XSHFT;
kzoff=i1-ZSHFT;
kyoff=i1-YSHFT;
iiXPLT=MyGraph->xv[ip];
iiYPLT=MyGraph->yv[ip];
iiZPLT=MyGraph->zv[ip];
set_linestyle(MyGraph->color[ip]);
oldxpl=data[iiXPLT][kxoff];
oldypl=data[iiYPLT][kyoff];
oldzpl=data[iiZPLT][kzoff];
for(i=i1;i<i2;i++){
{
xpl=data[iiXPLT][kxoff];
ypl=data[iiYPLT][kyoff];
zpl=data[iiZPLT][kzoff];
}
if(TORUS==1)
{
if (fabs(oldxpl-xpl)>(float)(.5*TOR_PERIOD))oldxpl=xpl;
if (fabs(oldypl-ypl)>(float)(.5*TOR_PERIOD))oldypl=ypl;
if (fabs(oldzpl-zpl)>(float)(.5*TOR_PERIOD))oldzpl=zpl;
}
if(MyGraph->ColorFlag!=0&&i>i1){
for(j=0;j<=NEQ;j++){
v1[j]=data[j][i];
v2[j]=data[j][i-1];
}
comp_color(v1,v2,NODE,
(float)fabs(data[0][i]-data[0][i+1]));
} /* ignored by postscript */
/* if(MyGraph->line[ip]<0)
goto noplot; */
if(MyGraph->line[ip]<=0){
PointRadius=-MyGraph->line[ip];
if(MyGraph->ThreeDFlag==0)point_abs(xpl,ypl);
else point_3d(xpl,ypl,zpl);
}
else {
if(MyGraph->ThreeDFlag==0)
line_abs(oldxpl,oldypl,xpl,ypl);
else
line_3d(oldxpl,oldypl,oldzpl,xpl,ypl,zpl);
}
noplot:
oldxpl=xpl;
oldypl=ypl;
oldzpl=zpl;
kxoff++;
kyoff++;
kzoff++;
}
}
}
/* Sets the color according to the velocity or z-value */
comp_color( v1,v2,n, dt)
float *v1, *v2;
int n;
float dt;
{
int i,cur_color;
float sum;
float min_scale=(float)(MyGraph->min_scale);
float color_scale=(float)(MyGraph->color_scale);
if(MyGraph->ColorFlag==2){
sum=v1[MyGraph->ColorValue];
/* printf(" %d:sum=%g \n",MyGraph->zv[0],sum); */
}
else
{
for(i=0,sum=0.0;i<n;i++)sum+=(float)fabs((double)(v1[i+1]-v2[i+1]));
sum=sum/(dt);
}
cur_color=(int)((sum-min_scale)*(float)color_total/color_scale);
/* printf("min=%f max=%f col = %d val = %f \n",min_scale,color_scale,
cur_color,sum); */
if(cur_color<0)cur_color=0;
if(cur_color>color_total)cur_color=color_total-1;
cur_color+=FIRSTCOLOR;
set_color(cur_color);
ps_do_color(cur_color);
}
shoot(x,xg,evec,sgn)
double *x,*xg,*evec;
int sgn;
{
int i;
double t=0.0;
for(i=0;i<NODE;i++)
x[i]=xg[i]+sgn*evec[i]*DELTA_T*.1;
i=1;
integrate(&t,x,TEND,DELTA_T,0,NJMP,&i);
ping();
}
stop_integration()
{
/* set some global error here... */
DelayErr=1;
}
stor_full()
{
char ch;
if(!Xup){
printf(" Storage full -- increase maxstor \n");
return(0);
}
if(FOREVER)goto ov;
ping();
ch=(char)TwoChoice("YES","NO","Storage full: Overwrite?",
"yn");
if(ch=='y')
{
ov:
storind=0;
return(1);
}
return(0);
}
syntax highlighted by Code2HTML, v. 0.9.1