00001 #include<iostream>
00002 #include<fstream>
00003 #include<stdlib.h>
00004
00005 #include "def.h"
00006 #include "model.h"
00007 #include "nr.h"
00008
00009 using namespace std;
00010
00011
00012
00013 GlobExp *parseopts(int argc, char *argv[],Glob *globs,char *outstr);
00014
00015 void initialise(GlobExp ex[],Glob *globs,int simit);
00016
00017 void tabulateValues (Glob *globs,GlobExp *ex,double t0, double t1, double *t, long n, double *state,
00018 double **y, double ***dmds, double ***dmdp, double **dyds,double **dydp);
00019
00020 void outSimit(GlobExp ex[],Glob *globs,double *t,long n,double **y);
00021
00022 void freeMem(GlobExp *ex,Glob *globs,int simit);
00023
00024
00025
00026
00027
00028
00029 ofstream *dbg;
00030
00031 main(int argc, char *argv[])
00032 {
00033 long k,i,j;
00034 char outstr[100];
00035 Glob globs;
00036 GlobExp *ex;
00037 double *t,time,**y,*state;
00038 long n;
00039
00040 ex=parseopts(argc,argv,&globs,outstr);
00041
00042
00043 dbg=new ofstream("simit.dbg");
00044 if (!dbg)
00045 {
00046 cerr << "Error opening DEBUG file.\n";
00047 exit(1);
00048 }
00049 dbg->precision(4);
00050 dbg->unsetf(ios::floatfield);
00051
00052
00053 *dbg << DefModelDescription << "\n";
00054 *dbg << NPARAMS << " parameter(s):";
00055 for (k=1; k<=NPARAMS; ++k)
00056 {
00057 *dbg << " " << ex[1].par[k];
00058 }
00059 *dbg << "\n";
00060 dbg->flush();
00061
00062
00063
00064 *dbg << "\n\n";
00065 if (ex[1].y0)
00066 {
00067 *dbg << NEQNS << " starting value(s):";
00068 for (k=1; k<=NEQNS; ++k)
00069 *dbg << ' ' << ex[1].y0[k];
00070 *dbg << "\n";
00071 }
00072 else
00073 {
00074 *dbg << "No starting values specified\n";
00075 ex[1].y0=dvector(1,ex[1].nvar);
00076 for (k=1; k<=NEQNS; ++k)
00077 ex[1].y0[k]=DefYValues[k-1];
00078 }
00079 *dbg << "\n";
00080 dbg->flush();
00081
00082
00083 if(globs.dt <= 0.)
00084 {
00085 cerr << "illegal integration step: " << globs.dt << " \n";
00086 exit(1);
00087 }
00088 *dbg << "Integration step: " << globs.dt << endl;
00089
00090 time=ex[1].fitstart;
00091 n=0;
00092
00093 while(time < ex[1].fitend-globs.eps)
00094 {
00095 time+=globs.dt;
00096 n++;
00097 }
00098 n++;
00099 *dbg << "Number of time points : " << n << endl;
00100
00101 t=dvector(1,n);
00102 t[1]=ex[1].fitstart;
00103 t[n]=ex[1].fitend;
00104
00105 for(k=2;k<=n-1;k++)
00106 t[k]=t[k-1]+globs.dt;
00107 #ifdef PLOTTIMEPOINTS
00108 *dbg << "Time points : ";
00109 for(k=1;k<=n;k++)
00110 {
00111 *dbg << t[k];
00112 if(k!=n)
00113 *dbg << ";";
00114 }
00115 *dbg << endl;
00116 #endif
00117 dbg->flush();
00118
00119 y=dmatrix(1,n,1,ex[1].nvar);
00120 state=dvector(1,ex[1].nvar);
00121 for(k=1;k<=ex[1].nvar;k++)
00122 state[k]=ex[1].y0[k];
00123
00124 initialise(ex,&globs,TRUE);
00125
00126 try
00127 {
00128 tabulateValues(&globs,&ex[1],t[1],t[n],t,n,state,y,NULL,NULL,NULL,NULL);
00129 }
00130 catch (int i)
00131 {
00132 exit(1);
00133 }
00134
00135 outSimit(ex,&globs,t,n,y);
00136 cout << "creating " << ex[1].fileName << endl;
00137
00138 freeMem(ex,&globs,TRUE);
00139 }