PROGRAM Morse IMPLICIT NONE REAL*8 :: Y(2), DT, Time, X0, P0, Energy, E0 INTEGER :: I, Icnt = 0 CHARACTER (20) :: FileName PRINT *, "What should the ouptput file be called?" READ *, FileName PRINT *, "X0 = ?" READ *, X0 PRINT *, "Energy = ?" READ *, E0 OPEN(20, File = FileName) P0 = exp(-X0) * SQRT(2. * ((E0 - 1.) * exp(2. * X0) + 2. * exp(X0) - 1.)) PRINT *, P0 ! Initial conditions Y(1) = X0 Y(2) = P0 ! Time step DT = 1.0E-2 DO I = 1, 10001 Icnt = Icnt + 1 Energy = Y(2) * Y(2) / 2. + (1. - exp(-Y(1))) * (1. - exp(-Y(1))) Time = (I - 1.) * DT IF (Icnt == 5) THEN WRITE (*,*) Time, Y(1), Y(2), Energy WRITE (20,*) Time, Y(1), Y(2), Energy CALL RK4 (Y, Time, DT) Icnt = 0 END IF END DO STOP CONTAINS ! ------------------------------------------------------------------------------ SUBROUTINE FDT(Y, Time, F) REAL*8 :: Y(2), F(2), Time F(1) = Y(2) F(2) = -2. * (exp(-Y(1)) - exp(-2. * Y(1))) RETURN END SUBROUTINE FDT ! ------------------------------------------------------------------------------ SUBROUTINE RK4 (Y, Time, DT) REAL*8 :: YH(2), Y(2), DT, Y1DT(2), Y2DT(2), Y3DT(2), Y4DT(2), Time INTEGER :: L ! Calculate k1 = FDT(t, y(t)) CALL FDT(Y, Time, Y1DT) ! ------------------------------------------------------------------------------ ! Calculate k2 = FDT(t + DT, y(t)+ 0.5*DT K1) DO L = 1, 2 YH(L) = Y(L) + 0.5 * DT * Y1DT(L) END DO CALL FDT(YH, Time + DT / 2., Y2DT) ! ------------------------------------------------------------------------------ ! Calculate k3 = FDT(t + DT/2, y(t)+ 0.5*DT K2) DO L = 1, 2 YH(L) = Y(L) + 0.5 * DT * Y2DT(L) END DO CALL FDT(YH, Time + DT / 2., Y3DT) ! ------------------------------------------------------------------------------ ! Calculate k4 = FDT(t + DT, y(t) + dt K3) DO L = 1, 2 YH(L) = Y(L) + DT * Y3DT(L) END DO CALL FDT(YH, Time + DT, Y4DT) ! ------------------------------------------------------------------------------ ! RK4 step DO L = 1, 2 Y(L) = Y(L) + (Y1DT(L) + 2. * Y2DT(L) + 2. * Y3DT(L) + Y4DT(L)) * DT / 6. END DO RETURN END SUBROUTINE RK4 END PROGRAM