#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[])) |
|
Definition at line 8 of file integrateRK.cc. |
|
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 }
|