!------------------------------------------------------------------------------
!PHY2063 Assignment: Modified Euler Method Approximator for coupled ODEs
!URN 6484960
!October 9 2018       Numerical Physics Lab Group B (TUESDAY 1000-1300)
!------------------------------------------------------------------------------
PROGRAM assignment1 !https://www.desmos.com/calculator/4lrp1jh2b7
    IMPLICIT NONE

    INTEGER, PARAMETER                       :: dp = KIND(1.d0)
    ! Helps with setting variables as double-precision reals

    REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: x, y, vx, vy, t
    ! There are 3 one-dimensional arrays: used to store values of x(t),
    ! v(t) & t

    REAL(KIND=dp), PARAMETER                 :: h = 0.005
    ! The distance between each consecutive value of t at which x(t) and
    ! v(t) will be evaluated

    REAL(KIND=dp), PARAMETER                 :: finalt = 50.0
    ! Tells the program the maximum t-value to reach

    REAL(KIND=dp) :: dxdt, dydt, dvxdt, dvydt, xmid, ymid, vxmid, vymid,&
                     dxdtmid, dvxdtmid, dydtmid, dvydtmid
    ! Declares re-usable variables for time, dx/dt and dv/dt, as well as
    ! the midpoints between two consecutive x and v points, and their
    ! time derivatives

    INTEGER                                  :: i
    !A simple index to use in loops

    WRITE(6,*) "This program will solve the Assignment 1 coupled ODEs, "
    WRITE(6,*) "using the modified Euler method, and a t-step of 5e-3"

    ALLOCATE( x(0:(NINT(finalt/h))))
    ALLOCATE( y(0:(NINT(finalt/h))))
    ALLOCATE(vx(0:(NINT(finalt/h))))
    ALLOCATE(vy(0:(NINT(finalt/h))))
    ALLOCATE( t(0:(NINT(finalt/h))))
    ! Looks at the finalt and h parameters and initializes all 1D arrays
    ! to the correct size

     x(0) = -2.5
     y(0) =  0.0
    vx(0) =  2.0
    vy(0) =  2.0
     t(0) =  0.0
    !Sets the boundary conditions as given

    OPEN (15, FILE='data.dat')

    WRITE(15,*) " #  t                         x(t)                      y(t)  &
                         &                    vx(t)                     vy(t)"
    WRITE(15,*) t(0), x(0), y(0), vx(0), vy(0)
    !WRITE(15,'(1x, f6.3, 3x, es12.5E2, 3x, es12.5E2)') t(0),x(0),v(0)
    !Prints the column headings as well as the values of x(t) and v(t)

    DO i = 1, (NINT(finalt/h))
    !Main loop

        t(i)     = t(i-1) + h
        ! Saves the current value for time to the appropriate index in
        ! the array

        dxdt     = vx(i-1)
        dvxdt    = -5*x(i-1) + (0.5*COS(2*t(i-1)))

        dydt     = vy(i-1)
        dvydt    = -5*y(i-1) + (0.5*COS(2*t(i-1)))
        ! As dictated by the problem sheet

        xmid     = x(i-1) + 0.5*h*dxdt
        vxmid    = vx(i-1) + 0.5*h*dvxdt

        ymid     = y(i-1) + 0.5*h*dxdt
        vymid    = vy(i-1) + 0.5*h*dvydt
        ! As we must do whenever we use the modified Euler method

        dxdtmid    = vxmid
        dvxdtmid   = -5*xmid + 0.5*COS(2*(t(i-1)+(0.5*h)))

        dydtmid    = vymid
        dvydtmid   = -5*ymid + 0.5*COS(2*(t(i-1)+(0.5*h)))
        ! Similar to what we did for dx and dv

        x(i)     = x(i-1)  + h*dxdtmid
        y(i)     = y(i-1)  + h*dydtmid
        vx(i)    = vx(i-1) + h*dvxdtmid
        vy(i)    = vy(i-1) + h*dvydtmid
        ! Saves the final values calculated for x and v into their
        ! respective arrays

        WRITE(15,*) t(i), x(i), y(i), vx(i), vy(i)
        ! WRITE(15,'(1x,f6.3,3x,es12.5E2,3x,es12.5E2)') t(i),x(i),v(i)
        ! Prints the final values that were just calculated, along with
        ! the time
    END DO
    CLOSE(15)
    DEALLOCATE(x,y,vx,vy,t)
    ! The loop has ended, the arrays are removed from memory, and the
    ! file will no longer be written to

    WRITE(6,*) ''
    WRITE(6,*) 'The results can be found in data.dat in the same folder&
                & as this program.'
    WRITE(6,*) 'You may plot them by running gnuplot and then using the&
                & command:'
    WRITE(6,*) "plot 'data.dat' u 1:2 w l, 'data.dat' u 1:3 w l"

END PROGRAM assignment1