\section{/home/peifer/diffit/modules/sim\-Init.cc File Reference}
\label{simInit_8cc}\index{/home/peifer/diffit/modules/simInit.cc@{/home/peifer/diffit/modules/simInit.cc}}
{\tt \#include $<$iostream$>$}\par
{\tt \#include $<$fstream$>$}\par
{\tt \#include $<$math.h$>$}\par
{\tt \#include $<$string.h$>$}\par
{\tt \#include $<$stdio.h$>$}\par
{\tt \#include \char`\"{}../def.h\char`\"{}}\par
{\tt \#include \char`\"{}../model.h\char`\"{}}\par
{\tt \#include \char`\"{}../nr.h\char`\"{}}\par
\subsection*{Functions}
\begin{CompactItemize}
\item 
void {\bf tabulate\-Values} ({\bf Glob} $\ast$globs, {\bf Glob\-Exp} $\ast$ex, double t0, double t1, double $\ast$t, long n, double $\ast$state, double $\ast$$\ast$y, double $\ast$$\ast$$\ast$dmds, double $\ast$$\ast$$\ast$dmdp, double $\ast$$\ast$dyds, double $\ast$$\ast$dydp)
\item 
void {\bf sim\-Init} ({\bf Glob\-Exp} $\ast$ex, {\bf Glob} $\ast$globs)
\end{CompactItemize}


\subsection{Function Documentation}
\index{simInit.cc@{sim\-Init.cc}!simInit@{simInit}}
\index{simInit@{simInit}!simInit.cc@{sim\-Init.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void sim\-Init ({\bf Glob\-Exp} $\ast$ {\em ex}, {\bf Glob} $\ast$ {\em globs})}\label{simInit_8cc_1c8ff8f8413349876df41981e468ea78}




Definition at line 21 of file sim\-Init.cc.

References dmatrix(), dvector(), Glob\-Exp::n\-Points, Glob\-Exp::nvar, Glob\-Exp::y0, and Glob\-Exp::y\-Try.

\footnotesize\begin{verbatim}22 {
23   long i,j;
24   long nvar=ex->nvar;
25   long nPoints=ex->nPoints;
26   long idum=-42;
27   double *t=dvector(1,2);
28   double **y=dmatrix(1,2,1,nvar);
29   
30    if (ex->y0) 
31     {
32       for (i=1; i<=nvar; ++i) 
33         ex->yTry[1][i]=ex->y0[i];
34     } 
35    else
36      {
37         for (i=1; i<=nvar; ++i) 
38           {
39             ex->yTry[1][i]=DefYValues[i-1];
40           }
41      }
42 
43    //integrating 
44 
45    for(i=2;i<=nPoints;i++)
46      {
47        for(j=1;j<=nvar;j++)
48          ex->yTry[i][j]=ex->yTry[i-1][j];
49        t[1]=ex->mesh[i-1];
50        t[2]=ex->mesh[i];
51        tabulateValues(globs,ex,t[1],t[2],t,2,ex->yTry[i],y,NULL,NULL,NULL,NULL);
52       for(j=1;j<=nvar;j++)
53         ex->yTry[i][j]+=globs->pert*ex->yTry[i][j]*gasdev(&idum);
54      }
55    
56 
57 
58    free_dvector(t,1,2);
59    free_dmatrix(y,1,2,1,nvar);
60 }
\end{verbatim}\normalsize 


\index{simInit.cc@{sim\-Init.cc}!tabulateValues@{tabulateValues}}
\index{tabulateValues@{tabulateValues}!simInit.cc@{sim\-Init.cc}}
\subsubsection{\setlength{\rightskip}{0pt plus 5cm}void tabulate\-Values ({\bf Glob} $\ast$ {\em globs}, {\bf Glob\-Exp} $\ast$ {\em ex}, double {\em t0}, double {\em t1}, double $\ast$ {\em t}, long {\em n}, double $\ast$ {\em state}, double $\ast$$\ast$ {\em y}, double $\ast$$\ast$$\ast$ {\em dmds}, double $\ast$$\ast$$\ast$ {\em dmdp}, double $\ast$$\ast$ {\em dyds}, double $\ast$$\ast$ {\em dydp})}\label{simInit_8cc_49291cd052f6e6ad7da733e5175315f5}




Definition at line 108 of file int\-ODE.cc.

References call\_\-odessa(), dfill0(), dmatrix(), Glob::do\-P, dvector(), Glob::eps, init\-Yt(), integrate\-RK(), Glob::integrator, max(), Glob::maxstp, Glob::npar, Glob\-Exp::nvar, Glob\-Exp::par, Glob::rkqs\_\-ign, sens\-Derivs(), Glob::stiff, and TRUE.

\footnotesize\begin{verbatim}110 {
111   // integrate from t0 to t1, storing values at t[1..n]
112   // t[n] may be ==t1 or <t1
113   // in the latter case we will integrate from t[n] to t1 in the
114   // final loop k=n+1 without storing any data points
115   // If dmds!=NULL, compute derivatives as well
116   // and store them in dmds,dmdp,dyds,dydp
117   // with respect to observation, see odin.doc: Observation function
118   
119   //BEGIN von Felix veraendert
120   //   double hmin=max(fabs(t0),fabs(t1))*1e-8;
121   double hmin = max (fabs (t0), fabs (t1)) * 1e-15;
122   // Das legt die minimale Schrittweite, die ODESSA erlaubt wird fest.
123   // hmin sollte so gross sein, dass es noch nicht zu Abschneidefehlern
124   // beim Aufaddieren zum Zeitpunkt kommt, deshalb die max(...)-Konstruktion.
125   // Der Faktor 1e-15 gewaehrleistet, dass noch richtig kleine Schritte gemacht
126   // werden duerfen, aber man trotzdem von der rel. Maschinengenauigkeit (im Bereich
127   // von 1e-18 noch einigermassen weit entfernt bleibt. 
128   // Man koennte das in Zukunft noch in einen externen Parameter stecken.
129   //END von Felix veraendert
130 
131   long nPExp=globs->npar;
132   long nVar=ex->nvar;
133   long nobs=NOBS;
134   int generic;
135   double *p=ex->par;
136   double eps=globs->eps;  
137   double h1 = (t1 - t0) * .01;
138   long nok,nbad;
139 
140   initInt(globs,ex);
141   if (n > 0 && (t0 > t[1] || t1 < t[n]))
142     {
143       cerr << "tabVal: t[1..n] exceeds [t0;t1]" << endl;
144       throw 1;
145     }
146 
147   if (globs->integrator==3)
148     {                           // said -odessa
149 #ifdef ODESSA
150       double &rtol = eps;
151 //BEGIN von Felix veraendert
152 //       double atol=eps*1e-6;  
153       double atol = eps;
154 // Fehlerkontrolle fuer ODESSA: Das Spiel mit relativem und absolutem LOKALEN(!)
155 // Fehler sollte man unbedingt noch in externe Parameter stecken. Doku dazu
156 // siehe odessa.f: ATOL und RTOL
157 //END von Felix veraendert
158       // odeint uses TINY=1e-30 in a similar way as odessa's atol, 
159       // but atol<1e-16 or so results in step size underflow
160       
161       int odeint_MAXSTP=globs->maxstp;
162       int d_flag=1;
163       int stif=globs->stiff;
164       char doPStr[globs->npar+1];
165       
166       for(long i=1;i<=globs->npar;i++)
167         doPStr[i]=(char)globs->doP[i];
168       
169       call_odessa(globs,ex,nVar, nPExp, doPStr,
170          d_flag, odeint_MAXSTP, stif,
171          TRUE, n, t, t0, t1, state, y,
172          dmds, dmdp, dyds, dydp, p, rtol, atol, hmin, t1 - t0, h1);
173 
174       return;
175 #else
176       cerr << "-odessa used but not linked" << endl;
177       throw 1;
178 #endif // ODESSA
179     }                           // if odessa
180       double *gy=dvector(1,nobs);
181       long i, j, k;
182       double xx1 = t0, xx2;
183       
184       if (dmds)
185         {                           // compute sensitivities
186           // create extended state vector Y
187           double **Yt = dmatrix (0, nPExp + nVar, 1, nVar);
188           initYt (Yt[0], state);
189           
190           // chain rule for extra observables
191           double **dgdy, **dgdp;
192           
193           dgdy = dmatrix (1, nobs, 1, nVar);
194           dgdp = dmatrix (1, nobs, 1, nPExp);
195           dfill0 (dgdy[1], nobs * nVar);
196           dfill0 (dgdp[1], nobs * nPExp);
197           
198           
199           for (k = 1; k <= n + 1; k++)
200             {
201               xx2 = (k > n) ? t1 :  // end point
202                 t[k];               // data point
203               
204               globs->rkqs_ign = (nPExp + nVar) * nVar;
205               if (xx2 > xx1 && nVar > 0)
206                 {                 
207                   if(globs->integrator==1)
208                     {
209                       //Runge-Kutta integation
210                       integrateRK(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,
211                                   eps, h1, hmin, &nok, &nbad, sensDerivs);
212                     }
213                   else if(globs->integrator==2)
214                     {
215                       cerr << "CVODES not optimized skip integrator" << endl;
216                       throw 1;
217             
218                       //integrateCVODES(globs,Yt[0], p, (1 + nPExp + nVar) * nVar, xx1, xx2,eps,TRUE);
219                     }
220                 }
221               // write state to y and sensitivities to dm/dp and dm/ds
222               if (k <= n)
223                 {
224                   dcopy (y[k], Yt[0], nVar);
225                   generic=observation(globs,ex,xx2, y[k],gy, p, dgdy, dgdp);
226                   for(i=1;i<=nobs;i++)
227                     y[k][i]=gy[i];
228                   
229                   for (i = 1; i <= nobs; i++)
230                     {
231                       for (j = 1; j <= nPExp; j++)
232                         {
233                           double &dest = dmdp[k][i][j];
234                           if (!generic)
235                             {
236                               dest = dgdp[i][j];
237                               for (long l = 1; l <= nVar; l++)
238                                 dest += dgdy[i][l] * Yt[j][l];
239                             }
240                           else
241                             dest = Yt[j][i];
242                         }           // for j
243                       for (j = 1; j <= nVar; j++)
244                         {
245                           double &dest = dmds[k][i][j];
246                           if (!generic)
247                             {
248                               dest = 0;
249                               for (long l = 1; l <= nVar; l++)
250                                 dest += dgdy[i][l] * Yt[nPExp + j][l];
251                             }
252                           else
253                             dest = Yt[nPExp + j][i];
254                         }           // for j
255                     }               // for i
256                 }                   // if k<=n
257               xx1 = xx2;
258             }                       // for k
259           
260           // copy state and its derivatives to state, dydp and dyds
261           dcopy (state, Yt[0], nVar);       // original state vector
262           // write sensitivities to dy/dp and dy/ds
263           for (i = 1; i <= nVar; i++)
264             {
265               for (j = 1; j <= nPExp; j++)
266                 dydp[i][j] = Yt[j][i];
267               for (j = 1; j <= nVar; j++)
268                 dyds[i][j] = Yt[nPExp + j][i];
269             }
270           
271           free_dmatrix (dgdp,1 ,nobs , 1, nPExp);
272           free_dmatrix (dgdy,1 ,nobs, 1, nVar);
273           free_dmatrix (Yt, 0, nPExp + nVar, 1, nVar);
274         }
275       else
276         {                           // !dmds
277           for (k = 1; k <= n + 1; k++)
278             {
279               xx2 = (k <= n) ? t[k] : t1;
280               
281               globs->rkqs_ign = 0;
282               if (xx2 > xx1 && nVar > 0)
283                 { 
284                   if(globs->integrator==1)
285                     {
286                       //Runge-Kutta integration
287                       integrateRK(globs,state, p, nVar, xx1, xx2,
288                                   eps, h1, hmin, &nok, &nbad, ode);
289                     }
290                   else if(globs->integrator==2)
291                     {
292                       cerr << "CVODES not optimized skip integrator" << endl;
293                       throw 1;
294 
295                       //integrateCVODES(globs,state, p, nVar, xx1, xx2,eps,FALSE);
296                     }
297                 }
298               if (k <= n)
299                 {
300                   dcopy (y[k], state, nVar);
301                   generic=observation (globs,ex,xx2, y[k],gy,p, NULL, NULL);
302                   for(i=1;i<=nobs;i++)
303                     y[k][i]=gy[i];
304                 }
305               xx1 = xx2;
306             }                       // for k
307         }           
308       // if dmds else
309       free_dvector(gy,1,nobs); 
310 }
\end{verbatim}\normalsize 


