/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 }

Generated on Mon Jan 29 17:09:13 2007 for Diffit by  doxygen 1.4.6