!--
!-- Fortran program to reproduce table 1.1, textbook 
!-- 
!-- Finite Difference Methods for Ordinary and Partial
!-- Differential Equations by Randall J. LeVeque
!--
!-- By Nasser M. Abbasi   sept 25, 2010
!-- Math 228A, UC DAVIS, Fall 2010
!-- 
!-- to help me learn Fortran
!--
!-- OUTPUT IS

!gfortran -std=f2003 -O2 -Wextra -Wall -pedantic -fcheck=all -fwhole-file  -funroll-loops 
! -ftree- vectorize -march=native -Wsurprising -Wconversion file.f90
!
!$ ./a.exe
!       h            D+             D_             D0             D3
!   1.0E-001    -4.2939E-02     4.1138E-02    -9.0005E-04     6.8207E-05
!   5.0E-002    -2.1257E-02     2.0807E-02    -2.2510E-04     8.6491E-06
!   1.0E-002    -4.2163E-03     4.1983E-03    -9.0050E-06     6.9941E-08
!   5.0E-003    -2.1059E-03     2.1014E-03    -2.2513E-06     8.7540E-09
!   1.0E-003    -4.2083E-04     4.2065E-04    -9.0050E-08     6.9979E-11


program main
    implicit none
    DOUBLE PRECISION, parameter :: exact_value = COS(1.0D0)
    DOUBLE PRECISION, parameter :: h(5)=(/0.1D0,0.05D0,0.01D0,0.005D0,0.001D0/)
    DOUBLE PRECISION  :: result(4) = 0.0D0
    INTEGER :: I
    DOUBLE PRECISION :: derivative_left,derivative_right,derivative_center,derivative_3rd_order
    
    write (*,'(5((7x),A))')  'h','     D+','      D_','      D0','      D3' 

    DO I = 1, 5 ,1              
       result(1) = derivative_right(h(I))     -  exact_value 
       result(2) = derivative_left(h(I))      -  exact_value 
       result(3) = derivative_center(h(I))    -  exact_value 
       result(4) = derivative_3rd_order(h(I)) -  exact_value 

       write (*,'(ES11.1E3,(4(4x,ES11.4E2)))') h(I), result 
    ENDDO

  
end program main

!--------------------
DOUBLE PRECISION FUNCTION derivative_left(h)
 implicit none
 DOUBLE PRECISION, intent(in) :: h
 derivative_left = ( sin(1.0D0) - sin(1.0D0 - h)  ) / h
end

!--------------------
DOUBLE PRECISION FUNCTION derivative_right(h)
 implicit none
 DOUBLE PRECISION, intent(in) :: h
 derivative_right = ( sin(1.0D0+h) - sin(1.0D0) ) / h
end

!--------------------
DOUBLE PRECISION FUNCTION derivative_center(h)
 implicit none
 DOUBLE PRECISION, intent(in) :: h
 derivative_center = ( sin(1.0D0+h) - sin(1.0D0-h) ) / (2.D0*h)
end

!--------------------
DOUBLE PRECISION FUNCTION derivative_3rd_order(h)
 implicit none
 DOUBLE PRECISION, intent(in) :: h
 derivative_3rd_order = ( 2.D0*sin(1.0D0+h) + 3.D0*sin(1.0D0) - 6.D0*sin(1.D0-h) + sin(1.D0-2.0D0*h) ) / (6.D0*h)
end