/home/peifer/diffit/modules/integrateRK.cc File Reference

#include <math.h>
#include <iostream>
#include "../nr.h"
#include "../def.h"

Go to the source code of this file.

Defines

#define TINY   1.0e-30

Functions

void integrateRK (Glob *globs, double ystart[], double p[], long nvar, double x1, double x2, double eps, double h1, double hmin, long *nok, long *nbad, void(*derivs)(double, double[], double[], double[]))


Define Documentation

#define TINY   1.0e-30
 

Definition at line 8 of file integrateRK.cc.


Function Documentation

void integrateRK Glob globs,
double  ystart[],
double  p[],
long  nvar,
double  x1,
double  x2,
double  eps,
double  h1,
double  hmin,
long *  nok,
long *  nbad,
void(*)(double, double[], double[], double[])  derivs
 

Definition at line 19 of file integrateRK.cc.

References dvector(), Glob::maxstp, nrerror(), and SIGN.

Referenced by tabulateValues().

00022 {
00023   long ode_fkmax = 0, ode_fkount = 0, ode_kmax = 0, ode_kount = 0;
00024   double ode_dxsav = 0, *ode_fxp = 0, *ode_xp = 0, **ode_yp = 0;
00025   long rkqs_ign = 0;
00026   long nstp,i;
00027   double xsav,x,hnext,hdid,h;
00028   double *yscal,*y,*dydx;
00029   long odeint_MAXSTP=globs->maxstp;
00030 
00031   yscal=dvector(1,nvar);
00032   y=dvector(1,nvar);
00033   dydx=dvector(1,nvar);
00034   x=x1;
00035   h=SIGN(h1,x2-x1);
00036   *nok = (*nbad) = ode_kount = 0;
00037   for (i=1;i<=nvar;i++) y[i]=ystart[i];
00038   if (ode_dxsav<0) nrerror("odeint:dxsav<0");
00039   if (ode_kmax>0 && ode_kount<ode_kmax-1) {
00040      ode_xp[++ode_kount]=x;
00041      for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
00042      xsav=x;
00043   } /* save first point, w.h. 19.03.99 */
00044   for (nstp=1;nstp<=odeint_MAXSTP;nstp++) {
00045     (*derivs)(x,y,dydx,p);
00046     for (i=1;i<=nvar;i++)
00047       yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
00048     if (ode_fkmax) {
00049       if (ode_fkount>=ode_fkmax) nrerror("ode_fkount exceeds ode_fkmax");
00050       ode_fxp[++ode_fkount]=x;
00051     }
00052     if (ode_kmax>0 && ode_kount<ode_kmax-1 && fabs(x-xsav)>ode_dxsav) {
00053       ode_xp[++ode_kount]=x;
00054       for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
00055       xsav=x;
00056     }
00057     if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
00058     rkqs(y,dydx,p,nvar,&x,h,eps,yscal,&hdid,&hnext,globs->rkqs_ign,derivs);
00059     if (hdid == h) ++(*nok); 
00060     else ++(*nbad);
00061     if ((x-x2)*(x2-x1) >= 0.0) {
00062       for (i=1;i<=nvar;i++) ystart[i]=y[i];
00063       if (ode_kmax>0 && ode_kount<ode_kmax-1) {
00064         ode_xp[++ode_kount]=x;
00065         for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
00066       }
00067       free_dvector(dydx,1,nvar);
00068       free_dvector(y,1,nvar);
00069       free_dvector(yscal,1,nvar);
00070       return;
00071     }
00072     if (fabs(hnext) <= hmin)
00073       {
00074         cerr << "Step size too small in odeint\n";
00075         throw 1;
00076       }
00077     h=hnext;
00078   }
00079   cerr << "Too many steps in routine odeint; try -maxstp <n> greater n\n";
00080   throw 1;
00081 }


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