- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/simit.cc

Go to the documentation of this file.

    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 //Begin definition of module prototypes
    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 //Module modules/freeMem.cc
    00022 void freeMem(GlobExp *ex,Glob *globs,int simit);
    00023
    00024
    00025 //End definition of module prototypes
    00026
    00027
    00028 // DEBUG stream
    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   //open DEBUG stream
    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);     // set output format for floating point data
    00050   dbg->unsetf(ios::floatfield);
    00051
    00052   //print debugging information
    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   //initial values
    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   //setting up time points
    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 }