!===========================================================
!
! File: sch_util.f90, module schutil
!
! Module that contains common procedures:
! integrate and rkstep, as well as the functions f1,f2
! Also defines common variables energy and dp
module          schutil
 use, intrinsic        :: iso_fortran_env
 use                   :: schV
 implicit none 
 integer , parameter   :: dp = real64
! ----------------------------------------------------------
 private
 public                :: integrate,rkstep
! ----------------------------------------------------------
contains
! ----------------------------------------------------------
!===========================================================
!----- trivial function: derivative of psi      
 function       f1(x,psi,psip) 
  real(dp) ::   f1
  real(dp) ::   x,psi,psip
  f1=psip
 end function   f1
!===========================================================
!----- the second derivative of wavefunction:
!psip(x)' = psi(x)'' = -(E-V) psi(x)
 function       f2(x,psi,psip)
  real(dp) ::   f2
  real(dp) ::   x,psi,psip
!----- Schroedinger eq: RHS
  f2 = (V(x)-energy)*psi
 end function   f2
!===========================================================
!Simpson's rule to integrate psi(x)*psi(x) for proper
!normalization. For n intervals of width dx (n even)
!Simpson's rule is:
!int(f(x)dx) = 
! (dx/3)*(f(x_0)+4 f(x_1)+2 f(x_2)+...+4 f(x_{n-1})+f(x_n))
!
!Input:   Discrete values of function psi(Nx)
!         Integration step dx
!Returns: Integral(psi(x)psi(x) dx)
!===========================================================
 function       integrate(psi,dx,Nx)
  real(dp) ::   integrate
  integer  ::   Nx
  real(dp) ::   psi(Nx),dx
!-------------
  real(dp) ::   int
  integer  ::   i
!----- zeroth order point:
  i     = 1
  int   = psi(i)*psi(i)
!----- odd  order points (i=k+1 is even):
  do i  = 2,Nx-1,2
   int  = int + 4.0_dp*psi(i)*psi(i)
  enddo
!----- even order points:
  do i  = 3,Nx-2,2
   int  = int + 2.0_dp*psi(i)*psi(i)
  enddo
!----- last point:
  i     = Nx
  int   = int + psi(i)*psi(i)
!----- measure normalization:
  int   = int*dx/3.0_dp
!----- final result:
  integrate = int
 end function   integrate
!===========================================================
!Subroutine rkstep(t,x1,x2,dt)
!Runge-Kutta Integration routine of ODE
!dx1/dt=f1(t,x1,x2) dx2/dt=f2(t,x1,x2)
!User must supply derivative functions:
!real(8) function f1(t,x1,x2)
!real(8) function f2(t,x1,x2) 
!Given initial point (t,x1,x2) the routine advnaces it
!by time dt.
!Input : Inital time t    and function values x1,x2
!Output: Final  time t+dt and function values x1,x2
!Careful!: values of t,x1,x2 are overwritten...
!===========================================================
 subroutine     rkstep(t,x1,x2,dt)
  real(dp) ::   t,x1,x2,dt
  real(dp) ::   k11,k12,k13,k14,k21,k22,k23,k24
  real(dp) ::   h,h2,h6

  h =dt                           !h =dt, integration step
  h2=0.5_dp*h                     !h2=h/2
  h6=h/6.0_dp                     !h6=h/6

  k11=f1(t,x1,x2)
  k21=f2(t,x1,x2)
  k12=f1(t+h2,x1+h2*k11,x2+h2*k21)
  k22=f2(t+h2,x1+h2*k11,x2+h2*k21)
  k13=f1(t+h2,x1+h2*k12,x2+h2*k22)
  k23=f2(t+h2,x1+h2*k12,x2+h2*k22)
  k14=f1(t+h ,x1+h *k13,x2+h *k23)
  k24=f2(t+h ,x1+h *k13,x2+h *k23)

  t =t+h
  x1=x1+h6*(k11+2.0_dp*(k12+k13)+k14)
  x2=x2+h6*(k21+2.0_dp*(k22+k23)+k24)

 end subroutine rkstep

end  module     schutil
!===========================================================

!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2021)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------