program main
implicit none

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

integer(4),parameter:: n=2       ! Anzahl der Himmelskörper
real(8)             :: r0(3,n)   ! aktuelle Position
real(8)             :: rm(3,n)   ! vorherige Position
real(8)             :: rp(3,n)   ! naechste Position
real(8)             :: rij(3)    ! Differenzvektor
real(8)             :: v(3,n)    ! Geschwindigkeit
real(8)             :: f(3)      ! Kraft
real(8)             :: G         ! all. GravKonstante
real(8)             :: m(n)      ! Massen der H.
real(8)             :: radius(n) ! Radien der H.
real(8)             :: dt        ! Zeitschritt
integer(4)          :: niter     ! Anzahl der Zeitschritte
integer(4)          :: iter      ! Index des Zeitschritts
real(8)             :: Etot      ! Gesamtenergie
real(8)             :: Ekin      ! kinetische Energie
real(8)             :: Epot      ! potentielle Energie

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

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

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

read(1000,*)dt,niter
read(1000,*)G

read(1000,*)m(1),radius(1)
read(1000,*)r0(:,1)
read(1000,*)v(:,1)

read(1000,*)m(2),radius(2)
read(1000,*)r0(:,2)
read(1000,*)v(:,2)

!== bereite Anfangsbedingungen vor ======================

rij(:)=r0(:,2)-r0(:,1)
f=G*m(1)*m(2)*(rij(:))/((sqrt(dot_product(rij(:),rij(:))))**3)
rm(:,1)=r0(:,1)-v(:,1)*dt+f*(dt**2)/m(1)

rij(:)=r0(:,1)-r0(:,2)
f=G*m(2)*m(1)*(rij(:))/((sqrt(dot_product(rij(:),rij(:))))**3)
rm(:,2)=r0(:,2)-v(:,2)*dt+f*(dt**2)/m(1)

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

do iter=1,niter

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

  Ekin=0.d0  

  rij(:)=r0(:,2)-r0(:,1)
  call calcepot(G,m(1),m(2),rij(:),Epot)
  rp(:,1)=r0(:,1)*2.d0 &
        - rm(:,1) &
        + G*(dt**2)*m(2)*(rij(:)) &
          /((sqrt(dot_product(rij(:),rij(:))))**3)
  
  v(:,1)=(rp(:,1)-rm(:,1))/(2*dt)
  Ekin=Ekin+(0.5d0*m(1)*dot_product(v(:,1),v(:,1)))

  rij(:)=r0(:,1)-r0(:,2)
  rp(:,2)=r0(:,2)*2.d0 &
        - rm(:,2) &
        + G*(dt**2)*m(1)*(rij(:)) &
          /((sqrt(dot_product(rij(:),rij(:))))**3)

  v(:,2)=(rp(:,2)-rm(:,2))/(2*dt)
  Ekin=Ekin+(0.5d0*m(2)*dot_product(v(:,2),v(:,2)))

  Etot=Ekin+Epot

  !== schreibe Energie und Koordinaten ==================
  write(*,fmt='("iter ",i10," von ",i10," berechnet!" )')iter,niter
  write(1001,'(I15)') n
  write(1001,'(I8,A6,I8)') iter,".Frame von ",niter

  write(1001,'(A8,3ES23.15)') "SONNE   ", r0(:,1)
  write(1001,'(A8,3ES23.15)') "ERDE    ", r0(:,2)

  write(1002,'(I8,3ES23.15)')iter,Ekin,Epot,Etot

  !== 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(G,m1,m2,rij,epot)
implicit none
real(8)   ,intent(in)  :: G
real(8)   ,intent(in)  :: m1
real(8)   ,intent(in)  :: m2
real(8)   ,intent(in)  :: rij(3)
real(8)   ,intent(out) :: epot

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

epot=-G*m1*m2*(1.0d0/sqrt(dot_product(rij,rij)))
return
end subroutine calcepot
