************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE ODEINT(NVAR,X1,X2,EPS,H1,HMIN,STEPR,DERIVS, C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C in C $ YINOUT, C ^^^^^^ C inout C $ NOK,NBAD) C ^^^^^^^^ C out C driver routine for integration of differential equations. C NVAR starting values are integrated from X1 to X2 with accuracy C EPS. H1 is a guessed first step size, HMIN the minimum allowed C step size, which can be zero. STEPR is the name of the stepper C routine to be used. DERIVS is a user-supplied routine for C calculating the right hand sides. C C The starting values YINOUT are replaced with the results. C C NOK and NBAD are the numer of good and bad steps taken. C C Based on ideas from NUMERICAL RECIPES. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C named constants as array bounds INTEGER KMAXX, NMAX PARAMETER (KMAXX=200, NMAX=50) ************************************************************************ C input variables in the parameter list C diemension of vector INTEGER NVAR C accuracy, guessed first step size, minimum step size, starting C and end point of integration, dependent variable DOUBLE PRECISION EPS, H1, HMIN, X1, X2, YINOUT(NVAR) ************************************************************************ C output variables in the parameter list C number of good and bad steps for bookkeeping INTEGER NOK, NBAD ************************************************************************ C output to called routines C all called routines C stepper routine, evaluation of right-hand sides EXTERNAL STEPR, DERIVS C variables C x an y for evaluation of right-hand sides, this result, step, C scaling vector for stepper routine DOUBLE PRECISION X, Y(NMAX), DYDX(NMAX), H, YSCAL(NMAX) ************************************************************************ C output from called routines C number of successful steps and next step size DOUBLE PRECISION HDID, HNEXT ************************************************************************ C variables in common blocks C bookkeeping (not used in this version, but put in in case someone C wants to debug or whatever) INTEGER KMAX, KOUNT DOUBLE PRECISION DXSAV, XP(KMAXX), YP(NMAX,KMAXX) COMMON /PATH/ DXSAV, XP, YP, KMAX, KOUNT C this is to keep track of (very rare) errors here and in related C routines INTEGER MSTEP COMMON /MERROR/ MSTEP C save until the next call SAVE /MERROR/ ************************************************************************ C internally used quantities C named constants C maximum number of steps INTEGER MAXSTP C just for numerical stability; see below DOUBLE PRECISION TINY PARAMETER (MAXSTP=10000, TINY=1.0D-30) C variables C DO loops INTEGER I, NSTP C for bookkeeping DOUBLE PRECISION XSAV ************************************************************************ ************************************************************************ C run-time initialisation C reset stepsize error to zero whenever routine is called MSTEP = 0 C KMAX not initialised in DATA since it's in a common block; C BLOCK DATA not worth the effort KMAX = 0 C DXSAV isn't needed if KMAX = 0 * DXSAV = 0 X = X1 H = SIGN(H1,X2-X1) NOK = 0 NBAD = 0 KOUNT = 0 ************************************************************************ ************************************************************************ C calculations C fill up work vector DO 100 I = 1, NVAR Y(I) = YINOUT(I) 100 CONTINUE C uncomment the next line if KMAX isn't 0; for KMAX = 0 neither C XSAV nor DXSAV is needed, so we don't need the next line C (This assures storage of the first step.) C IF (KMAX .GT. 0) XSAV = X - DXSAV*2.0 C at most NMAX steps DO 600 NSTP = 1, MAXSTP C get the right-hand sides CALL DERIVS(X,Y, DYDX) DO 200 I = 1, NVAR C scaling to monitor accuracy YSCAL(I) = ABS(Y(I)) + ABS(H*DYDX(I)) + TINY 200 CONTINUE C storage for internal bookkeeping IF (KMAX .GT. 0) THEN IF ( (ABS(X - XSAV)) .GT. (ABS(DXSAV)) ) THEN IF( KOUNT .LT. (KMAX - 1) ) THEN KOUNT = KOUNT + 1 XP(KOUNT) = X DO 300 I = 1, NVAR YP(I,KOUNT) = Y(I) 300 CONTINUE XSAV = X ENDIF ENDIF ENDIF C make sure the step size isn't too large IF ( ((X + H - X2)*(X + H - X1)) .GT. 0.0 ) H = X2 - X ************************************************************************ C abnormal return; set a flag and get out IF (MSTEP .EQ. 1) RETURN C This happens if there has been a failure in the stepper routine C on the previous call. ************************************************************************ C call the stepper if everything's OK CALL STEPR(NVAR,DYDX,H,EPS,YSCAL,DERIVS, X,Y, HDID,HNEXT) IF (HDID .EQ. H)THEN NOK = NOK + 1 ELSE NBAD = NBAD + 1 ENDIF C finished? IF ( ((X - X2)*(X2 - X1)) .GE. 0.0 ) THEN C replace input values with output values DO 400 I = 1, NVAR YINOUT(I) = Y(I) 400 CONTINUE IF (KMAX .NE. 0) THEN C save final step KOUNT = KOUNT + 1 XP(KOUNT) = X DO 500 I = 1, NVAR YP(I,KOUNT) = Y(I) 500 CONTINUE ENDIF ************************************************************************ C normal return; we have to get out now, since the error C conditions which follow this one are irrelevant if we're C successful, which means we've reached the next line RETURN ENDIF ************************************************************************ C abnormal return (stepsize smaller than minimum); C set a flag and get out IF( (ABS(HNEXT)) .LT. HMIN ) THEN MSTEP = 2 RETURN END IF ************************************************************************ C set the next H value H = HNEXT 600 CONTINUE ************************************************************************ C abnormal return (too many steps); set a flag and get out MSTEP = 3 RETURN ************************************************************************ ************************************************************************ END