00001 #include <math.h>
00002 #include <iostream>
00003 #include "../nr.h"
00004 #include "../def.h"
00005
00006 using namespace std;
00007
00008 #define TINY 1.0e-30
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 void integrateRK(Glob *globs,double ystart[], double p[], long nvar, double x1, double x2,double eps,
00020 double h1, double hmin, long *nok, long *nbad,
00021 void (*derivs) (double, double[], double[], double[]))
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 }
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 }
00082 #undef TINY