program main
implicit none

!== Definition der Variablen ============================

integer(4),parameter:: n=3    ! Dimension der Koordinaten
real(8)             :: r0(n)  ! aktuelle Position
real(8)             :: rm(n)  ! vorherige Position
real(8)             :: rp(n)  ! naechste Position
real(8)             :: v(n)   ! Geschwindigkeit
real(8)             :: f(n)   ! Kraft
real(8)             :: c(n,n) ! Kraftkonstanten
real(8)             :: m      ! Masse
real(8)             :: dt     ! Zeitschritt
real(8)             :: alpha  ! Reibungsfaktor
real(8)             :: Etot   ! Gesamtenergie
real(8)             :: Ekin   ! kinetische Energie
real(8)             :: Epot   ! potentielle Energie
integer(4)          :: niter  ! Anzahl der Zeitschritte
integer(4)          :: iter   ! Index des Zeitschritts
real(8)             :: a      ! Hilfsvariable

!********************************************************
!== oeffne Dateien fuer Ein- und Ausgabe ================

open(unit=1000,file='in.dat',form='formatted')
open(unit=1001,file='energies.dat',form='formatted')
open(unit=1002,file='r.dat',form='formatted')

!== lese Eingabedaten ===================================

read(1000,*)m,alpha,c
read(1000,*)dt,niter
read(1000,*)r0
read(1000,*)v
niter=1000/dt
!== bereite Anfangsbedingungen vor ======================

rm(:)=r0(:)-v(:)*dt*(1+a)-0.5d0/m*matmul(c(:,:),r0(:))*dt**2

!== Das ist die Schleife ueber die Zeitschritte =========

do iter=1,niter

  !== berechne potentielle Energie und Kraefte ==========

  call calcepot(n,c,r0,epot,f)

  !== propagiere ========================================

  a=0.5d0*alpha*dt/m
  rp(:)=r0(:)*2.d0/(1+a) &
        - rm(:)*(1-a)/(1+a) &
        + f(:)*dt**2/m/(1+a)

  !== schreibe Energie und Koordinaten ==================
  v(:)=(rp(:)-rm(:))/(2.d0*dt)
  ekin=0.5d0*m*dot_product(v,v)
  etot=ekin+epot
  write(*,fmt='("iter ",i10," etot ",e10.5 ," ekin "&
          &,e10.5 ," epot ",e10.5)')iter,etot,ekin,epot
  write(1001,*)real(iter,kind=8)*dt,etot,ekin,epot
  write(1002,*)real(iter,kind=8)*dt,r0(:)

  !== Vertausche Koordinaten fuer naechsten Zeitschritt =

  rm(:)=r0(:)
  r0(:)=rp(:)
end do

!== nun schliesse die Dateien ===========================

close(1000)
close(1001)
close(1002)
stop
end program main

!........................................................

subroutine calcepot(n,c,r0,epot,f)
implicit none
integer(4),intent(in)  :: n
real(8)   ,intent(in)  :: c(n,n)
real(8)   ,intent(in)  :: r0(n)
real(8)   ,intent(out) :: epot
real(8)   ,intent(out) :: f(n)

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

epot=0.5d0*dot_product(r0(:),matmul(c(:,:),r0(:)))
f(:)=-matmul(c(:,:),r0(:))
return
end subroutine calcepot
