- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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