C
C <<<<<<<<<<<<<<<<<<<< PROMOEQN - PROMOEQN - PROMOEQN >>>>>>>>>>>>>>>>>>>>
C
SUBROUTINE PROMOEQN(T, Y, DYDT, NEQN, ERRMSG)
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C Calculate the derivative of Y at T for NEQN equations and store them in
C DYDT. Note that this subroutine has been set up for use with the pendulum
C ODE problem.
C
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION Y(NEQN), DYDT(NEQN)
CHARACTER ERRMSG*80
C
COMMON /PROMO/ GRAV, DRAGK, VEL, N
C
C Check input.
C
IF (NEQN .NE. 4) THEN
ERRMSG = 'NEQN must equal 4 in PROMOEQN.'
RETURN
ENDIF
C
C Note that our 4 differential equations (DYDT) are
C DYDT(1) = dx/dt = vx
C DYDT(2) = dvx/dt = -k * vx^n * vx/abs(v)
C DYDT(3) = dy/dt = vy
C DYDT(4) = dvy/dt = -g - k * vy^n * vy/abs(v)
C
C Note that our 4 corresponding dependent variables (Y) are:
C Y(1) = x(i)
C Y(2) = vx(i)
C Y(3) = y(i)
C Y(4) = vy(i)
C
C Calculate total velocity at time T (note that this is always positive):
C
VEL = DSQRT(Y(2)*Y(2) + Y(4)*Y(4))
IF (VEL .LE. 0.0D0) THEN
ERRMSG = 'Total velocity is less than or equal to zero.'
RETURN
ENDIF
C
C Calculate the derivatives at time T:
C
DYDT(1) = Y(2)
DYDT(2) = -DRAGK * Y(2)**N * Y(2) / VEL
DYDT(3) = Y(4)
DYDT(4) = -GRAV - DRAGK * Y(4)**N * Y(4) / VEL
C DYDT(1) = Y(2)
C DYDT(2) = -DRAGK * Y(2) * VEL
C DYDT(3) = Y(4)
C DYDT(4) = -GRAV - DRAGK * Y(4) * VEL
C
RETURN
END
C
C <<<<<<<<<< RK4_PROJMOTION - RK4_PROJMOTION - RK4_PROJMOTION >>>>>>>>>>>>
C
SUBROUTINE RK4_PROJMOTION(Y, DYDT, NEQN, T, H, YOUT, PROMOEQN,
+ ERRMSG)
C
C Solve up to 10 ordinary differential equations using 4th-order Runge-Kutta.
C
C Variable definitions:
C
C Y: Dependent variable array of size NEQN. [x(i), vx(i), y(i), vy(i)]
C DYDT: Value of time derivatives of our dependent variables.
C [dx(i)/dt, dvx(i)/dt, dy(i)/dt, dvy(i)/dt]
C NEQN: The number of equations being solve which is also the number
C of dependent variables stored in Y.
C T: The current value of the independent variable. [TIME(I)]
C H: The step size. [TSTEP]
C YOUT: The new values of Y (NEQN element array).
C ERMSG: A string containing any errors that are discovered in the
C calculations.
C
C The user supplies the subroutine PROMOEQN(T, Y, DYDT, NEQN, SUBR, ERRMSG)
C which returns derivatives DYDT at T.
C
IMPLICIT REAL*8 (A-H, O-Z)
C
EXTERNAL PROMOEQN
C
C Set the maximum number of functions.
C
PARAMETER (NMAX=10)
C
DIMENSION Y(NEQN), DYDT(NEQN), YOUT(NEQN)
DIMENSION YT(NMAX), DYT(NMAX), DYM(NMAX)
CHARACTER ERRMSG*80
C
C Check number of functions that were passed.
C
IF (NEQN .GT. NMAX) THEN
ERRMSG = 'Number of functions passed exceeds NMAX in '
1 // 'RK4_PROJMOTION.'
RETURN
ENDIF
C
HH = H * 0.5D0
H6 = H / 6.D0
TH = T + HH
C
C First step.
C
DO 11 I = 1, NEQN
YT(I) = Y(I) + HH*DYDT(I)
11 CONTINUE
C
C Second step.
C
CALL PROMOEQN(TH, YT, DYT, NEQN, ERRMSG)
IF (ERRMSG(1:1) .NE. ' ') RETURN
DO 12 I = 1, NEQN
YT(I) = Y(I) + HH*DYT(I)
12 CONTINUE
C
C Third step.
C
CALL PROMOEQN(TH, YT, DYM, NEQN, ERRMSG)
IF (ERRMSG(1:1) .NE. ' ') RETURN
DO 13 I = 1, NEQN
YT(I) = Y(I) + H*DYM(I)
DYM(I) = DYT(I) + DYM(I)
13 CONTINUE
C
C Fourth step.
C
CALL PROMOEQN(T+H, YT, DYT, NEQN, ERRMSG)
IF (ERRMSG(1:1) .NE. ' ') RETURN
DO 14 I = 1, NEQN
YOUT(I) = Y(I) + H6*(DYDT(I) + DYT(I) + 2.D0*DYM(I))
14 CONTINUE
C
RETURN
END