!=========================================================
!---------------------------------------------------------
program         euler1d
 implicit none
 integer, parameter :: P = 1000000
 real, dimension(P) :: t , x , v
 real               :: tf, x0, v0
 integer            :: Nt, i

 print *,'# Enter x0,v0,tf,Nt (t_i=0):'
 read  *,         x0,v0,tf,Nt
 print *,'# ',    x0,v0,tf,Nt

 if(Nt > P) stop 'Error: Nt > P'

 !--------------------------------------------------------
 call euler
 open(unit=20,file='euler.dat')
 do i = 1, Nt
  write(20,*) t(i), x(i), v(i)
 end do
 close(20)
 !--------------------------------------------------------
 call euler_cromer
 open(unit=20,file='euler_cromer.dat')
 do i = 1, Nt
  write(20,*) t(i), x(i), v(i)
 end do
 close(20)
 !--------------------------------------------------------
 call euler_verlet
 open(unit=20,file='euler_verlet.dat')
 do i = 1, Nt
  write(20,*) t(i), x(i), v(i)
 end do
 close(20)

!---------------------------------------------------------
contains
!---------------------------------------------------------

 function       accel(x)
  real :: x, accel
  accel = -10.0*sin(x)
 end function   accel
 !--------------------------------------------------------
 subroutine     euler
  integer :: i
  real    :: dt

  t(1)  = 0.0
  x(1)  = x0
  v(1)  = v0
  dt    = tf/(Nt-1)

  do i  = 2, Nt
   t(i) =                 (i-1)  * dt
   x(i) = x(i-1) +       v(i-1)  * dt
   v(i) = v(i-1) + accel(x(i-1)) * dt
  end do

 end subroutine euler
 !--------------------------------------------------------
 subroutine     euler_cromer
  integer :: i
  real    :: dt

  t(1)  = 0.0
  x(1)  = x0
  v(1)  = v0
  dt    = tf/(Nt-1)

  do i  = 2, Nt
   t(i) =                 (i-1)  * dt
   v(i) = v(i-1) + accel(x(i-1)) * dt ! Statement goes first
   x(i) = x(i-1) +       v(i  )  * dt ! Notice the difference
  end do


 end subroutine euler_cromer
 !--------------------------------------------------------
 subroutine     euler_verlet
  integer :: i
  real    :: dt, dt2,idt2
  real    :: xx0
  
  t(1)    = 0.0
  x(1)    = x0
  v(1)    = v0
  dt      = tf  / (Nt-1)

  dt2     = dt  * dt
  idt2    = 0.5 / dt ! = 1 / ( 2 dt )

  !One step backwards:
  xx0     = x(1) - v(1) * dt  + accel(x(1)) * dt2 / 2.0
  t(2)    = dt
  x(2)    = 2.0  * x(1) - xx0 + accel(x(1)) * dt2

  !Iterations: start from i = 3
  do  i   = 3, Nt
   t(i  ) = (i-1) * dt
   x(i  ) = 2.0 * x(i-1) - x(i-2) + accel(x(i-1)) * dt2
   v(i-1) = (x(i) - x(i-2)) * idt2
  end do

  !Last step:
  v(Nt) = (x(Nt) - x(Nt-1))/dt

 end subroutine euler_verlet
 !--------------------------------------------------------

end program     euler1d
!---------------------------------------------------------
!=========================================================

!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2022)
!  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/>.
!  -----------------------------------------------------------------------
