! Runge Kutta 2nd Order Solver (implicit and explicit)
!   Uses RK solvers to estimate solution to dy/dt = -y
!     with exact solution y(t) = exp(-t)
!
!   This is a good code to show numerical methods newcomers 
!     the difference between iterative residual convergence and 
!     solution convergence.
!
!   Warning: This code is optimized for readability 
!            *not* function evaluations.
!
! Copyright (c) 2008 Charles O'Neill
! 
! Permission is hereby granted, free of charge, to any person
! obtaining a copy of this software and associated documentation
! files (the "Software"), to deal in the Software without
! restriction, including without limitation the rights to use,
! copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following
! conditions:
! 
! The above copyright notice and this permission notice shall be
! included in all copies or substantial portions of the Software.
! 
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
! OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
! HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
! WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
! OTHER DEALINGS IN THE SOFTWARE.

module Rkiteration
  implicit none

  type RK2
    real :: A(2,2)     !  In Butcher Table
    real :: B(2)       !                       C A
    real :: C(2)       !                         B
  end type

contains 

  subroutine ERK2(A,B,C)
    ! Setup Butcher Table for Explicit 2nd order Runge Kutta
    real :: A(2,2)
    real :: B(2)
    real :: C(2)  
    A(1,1) = 0.0;   A(1,2) = 0.0
    A(2,1) = 1.0;   A(2,2) = 0.0
    
    B(1) = 0.5
    B(2) = 0.5
    
    C(1) = 0.0
    C(2) = 1.0
    
  end subroutine

  subroutine IRK2(A,B,C)
    ! Setup Butcher Table for Implicit 2nd order Runge Kutta
    real :: A(2,2)
    real :: B(2)
    real :: C(2)  
    A(1,1) = 0.25;   A(1,2) = -0.25
    A(2,1) = 0.25;   A(2,2) = 5.0/12.0
    
    B(1) = 0.25;  B(2) = 0.75
    
    C(1) = 0.0
    C(2) = 2.0/3.0
    
  end subroutine
  

  subroutine testTimeIntegration
    integer :: iter, IMax
    real :: A(2,2)
    real :: B(2)
    real :: C(2)
    
    real :: zeta1, zeta2, yn1, yn
    real :: tn,tn1,t1,t2
    real :: h
    
    h = 0.1    ! Timestep (\Delta T)
    yn = 1.0   ! Initial Condition (Y_n)
    tn = 0.0   ! Initial Time (T_n)
    tn1= tn+h  ! New Time (T_{n+1})
    Imax = 20  ! Maximum number of iterations
 
    ! Implicit Setup
    !call IRK2(A,B,C)
    
    ! Explicit Setup
    call ERK2(A,B,C)

    ! Time at subpoint 1 and 2
    t1 = tn + C(1)*h
    t2 = tn + C(2)*h

    ! Init subpoint states
    zeta1 = yn 
    zeta2 = yn 
    
    ! Output:          Time, z1, z2, y_{n+1}, y(t_{n+1})
    write(*,*) 0,zeta1,zeta2, yn + h*( B(1)*f(t1,zeta1) + B(2)*f(t2,zeta2)), exp(-tn1)

    do iter=1,IMax
      ! Update state at subpoint 1
      zeta1 = yn + h*( A(1,1)*f(t1,zeta1) + A(1,2)*f(t2,zeta2) )  
      ! Update state at subpoint 2
      zeta2 = yn + h*( A(2,1)*f(t1,zeta1) + A(2,2)*f(t2,zeta2) )
      ! Output:  Time, z1, z2, y_{n+1}, y(t_{n+1})      
      write(*,*) iter,zeta1,zeta2, yn + h*( B(1)*f(t1,zeta1) + B(2)*f(t2,zeta2)), exp(-tn1)
    enddo
    
    ! Final State at new timestep
    yn1   = yn + h*( B(1)*f(t1,zeta1) + B(2)*f(t2,zeta2))
    
    ! Compare Old and New states with Actual Function
    write(*,*) tn, yn, exp(-tn)
    write(*,*) tn1, yn1, exp(-tn1)
  
  end subroutine

  ! Example da/dt function
  function f(t,a)
    real :: f, a,t
    f = -a +0.0*t
  end function

end module