!**********************************************************************!
!
! File: cromosc.f90 (06-Apr-2001)
!                   (27-May-2005)
!                   (27-Mar-2008)
!
! Harmonischer Oszillator mit dem Cromer-Verfahren
!
!   dv/dt = -x
!   dx/dt = v
!
! dt..........Zeitschritt
! itpr........Ausgabe alle itpr Zeitschritte
! itmax.......maximale Anzahl Zeitschritte
! x(0),v(0)...Startwerte
!
! Ausgabe auf File "cromosc.dat"
!
!**********************************************************************!

program cromosc

!**********************************************************************!

  implicit none

  integer,parameter::iout=2

  integer::it,itmax,itpr
  real::dt,t,v,vn,x,xn
  
  write(unit=*,fmt="(a)",advance="no") " dt,itpr,itmax="
  read(unit=*,fmt=*) dt,itpr,itmax
  write(unit=*,fmt="(a)",advance="no") " x(0),v(0)="
  read(unit=*,fmt=*) x,v
  
  open(unit=iout,file="cromosc.dat",status="replace",action="write", &
       form="formatted",position="rewind")
  
  t=0.0
  write(unit=iout,fmt="(tr1,f11.6,2(tr1,es14.7))") t,x,v
  
  do it=1,itmax

    vn=v-dt*x      ! Cromer
    xn=x+dt*vn     !

    t=it*dt
    x=xn
    v=vn

    if(modulo(it,itpr)==0) then
      write(unit=iout,fmt="(tr1,f11.6,2(tr1,es14.7))") t,x,v
    end if
  end do
  
  close(unit=iout)
  
end program cromosc

!**********************************************************************!