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