!
!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
!
MODULE Precision
      PUBLIC
      INTEGER, PARAMETER :: sp = KIND (0.0E0), dp = KIND (0.0D0)
      INTEGER, PARAMETER :: wp = sp ! Use single precision in this run
END MODULE Precision
!
MODULE IVPSolve
      USE Precision
      IMPLICIT NONE
!
CONTAINS
      SUBROUTINE Euler (f, y, n_unknowns, n_points, dx, x0, y0)
      ! The function f is passed as an argument:
         INTERFACE
            SUBROUTINE f (x, y, y_prime, n_unknowns)
               USE Precision
               INTEGER n_unknowns
               REAL (KIND=wp), INTENT(IN) :: x
               REAL (KIND=wp), DIMENSION (n_unknowns), INTENT(IN) :: y
               REAL (KIND=wp), DIMENSION (n_unknowns), INTENT(OUT) :: y_prime
            END SUBROUTINE f
         END INTERFACE
         INTEGER, INTENT (IN) :: n_unknowns, n_points
         ! Upon exit, this array should contain the solution
         REAL (KIND=wp), DIMENSION (n_unknowns, n_points), INTENT (OUT) :: y
         REAL (KIND=wp), INTENT (IN) :: dx, x0 ! Step size and initial x
         REAL (KIND=wp), DIMENSION (n_unknowns), INTENT (IN) :: y0 ! Initial y
!
         REAL (KIND=wp) :: x
         REAL (KIND=wp), DIMENSION (n_unknowns) :: y_prime
         INTEGER :: i
!
         x = x0
         y (:, 1) = y0
         DO i = 1, n_points - 1
            CALL f (x, y(:, i), y_prime, n_unknowns)
            y (:, i+1) = y (:, i) + y_prime * dx ! Euler's method
            x = x + dx
         END DO
!
      END SUBROUTINE Euler
END MODULE IVPSolve
!
MODULE Harmonic
      USE Precision
      IMPLICIT NONE
      PUBLIC
      REAL (KIND=wp) :: k = 1.0_wp
CONTAINS
      SUBROUTINE kx (x, y, dy, n_eq)
         USE Precision
         INTEGER, INTENT (IN) :: n_eq
         REAL (KIND=wp), INTENT (IN) :: x
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (IN) :: y
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (OUT) :: dy
         dy = (/ - k * y (2), y (1) /)
      END SUBROUTINE kx
END MODULE Harmonic
!
PROGRAM ODE_test
      USE Precision
      USE Harmonic
      USE IVPSolve
      USE SimpleGraphics
      IMPLICIT NONE      
!
      INTEGER, PARAMETER :: n_eq = 2, n_pts = 1000
      REAL (KIND=wp), DIMENSION (n_eq, n_pts) :: x
      REAL (KIND=wp), DIMENSION (n_pts) :: t
      REAL (KIND=wp) :: dt, t_min = 0.0_wp, t_max = 10.0_wp
      INTEGER :: i
      CHARACTER (50), DIMENSION (2) :: title ! The title of the plot
!
      dt = (t_max-t_min) / REAL (n_pts-1, wp)
      t = (/ (t_min+REAL(i-1, wp)*dt, i=1, n_pts) /)
      CALL Euler (f=kx, y=x, n_unknowns=n_eq, n_points=n_pts, dx=dt, x0=t_min, y0= (/ 0.0_wp, 1.0_wp /))
!
     ! This is an example of using Plot2D plot:
      title (1) = "x(t) for an anharmonic oscillator" ! First line of the title
      title (2) = "F(x)=-k*x+k/2*x^3" ! Second line of title
      CALL InitGraphics (file="Harmonic.png", file_type="PNG", plot_title=title, x_label="t", y_label="x(t)")
      CALL Plot2D (x=REAL(t, sp), y=REAL(x(1, :), sp), plot_spec="L-R", new_plot=.TRUE.)! x=f(t)
      CALL EndGraphics ()
!
END PROGRAM ODE_test
!