program main
implicit none

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

integer(4),parameter:: n=10      ! Anzahl der Himmelskoerper
integer(4)          :: i,j       ! Zaehlvariablen
character(8)        :: namen(n)  ! Liste mit Namen der H.
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)             :: A(3)      ! Hilfsvariable
real(8)             :: Etot      ! Gesamtenergie
real(8)             :: Ekin      ! kinetische Energie
real(8)             :: Epot      ! potentielle Energie
real(8)             :: Ehilf     ! Hilfsvariable pot 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

do i=1,n
  read(1000,*)namen(i)
  read(1000,*)m(i),radius(i)
  read(1000,*)r0(:,i)
  read(1000,*)v(:,i)
end do

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

  do i=1,n
    A(:)=(0.d0)
    do j=1,n
      if (i/=j) then
        rij(:)=r0(:,j)-r0(:,i)
        A(:)=A(:) + m(j)*(rij(:))/((sqrt(dot_product(rij(:),rij(:))))**3)
      end if
    end do
    rm(:,i)=r0(:,i) - dt*v(:,i) - (1/2)*G*(dt**2)*A
  end do
!== Das ist die Schleife ueber die Zeitschritte =========

do iter=1,niter

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

  Epot=0.d0
  Ekin=0.d0

  do i=1,n
    A(:)=(0.d0)
    do j=1,n
      if (i/=j) then
        rij(:)=r0(:,j)-r0(:,i)
        call calcepot(G,m(i),m(j),rij(:),Ehilf)
        Epot=Epot+Ehilf
        A=A(:) + m(j)*(rij(:))/((sqrt(dot_product(rij(:),rij(:))))**3)
      end if
    end do
    rp(:,i)=r0(:,i)*2.d0 - rm(:,i) + G*(dt**2)*A
    v(:,i)=(rp(:,i)-rm(:,i))/(2*dt)
    Ekin=Ekin+(0.5d0*m(i)*dot_product(v(:,i),v(:,i)))
  end do
  Etot=Epot+Ekin

  !== schreibe Energie und Koordinaten ==================

  if (mod(iter,10000)==0) then
  write(*,fmt='("iter ",i10," von ",i10," berechnet!" )')iter,niter
  end if

  write(1001,'(I15)') n
  write(1001,'(I8,A11,I8)') iter,".Frame von ",niter

  do i=1,n
    write(1001,'(A8,3ES23.15)') namen(i), r0(:,i)
  end do

  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,ehilf)
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) :: ehilf

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

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