! PROGRAM Planet
!
! Bewegung in einem Zentralpotential V(r)=c/r: 4 gekoppelte
! (Hamilton) Differentialgleichungen 1. Ordnung.
! Die Gleichungen beschreiben zum Beispiel die Bewegung eines Planeten
! um eine viel schwerere Sonne, also in guter Naeherung auch unser
! Sonnensystem. Variieren Sie die Anfangswerte. 

MODULE planet_parameters
   integer, parameter :: nvar=4
   real, parameter :: pi=3.14259265358
   real :: mass,cpot
END MODULE planet_parameters

PROGRAM Planet
   use planet_parameters
   implicit none
   
   integer::i,Nstep
   real :: time,v10,v20,h,x1max,x1min,x2max,x2min,Viewtime,x,y
   real, dimension(1:nvar)::yvar,yder,outyvar
   character (len=1)::antwort
   character (len=256)::gnucmd,gnudat
   
   write(*,*) 'Masse  (z.B.  1)     = '
   read(*,*) mass
   write(*,*) 'Potential = c/r:'
   write(*,*) 'c      (z.B.  -1)    = '
   read(*,*) cpot
   write(*,*) 'max(x) (z.B.  10)    = '
   read(*,*) x1max
   x2max=x1max
   x1min=-x1max
   x2min=-x2max
   
   write(*,*) 'Anfangswerte:'
   write(*,*) 'x1     (z.B.  -5)    = '
   read(*,*) yvar(1)
   write(*,*) 'v1     (z.B.   0)    = '
   read(*,*) v10
   yvar(2)=v10*mass
   write(*,*) 'x2     (z.B.   0)    = '
   read(*,*) yvar(3)
   write(*,*) 'v2     (z.B.-0.3)    = '
   read(*,*) v20
   yvar(4)=v20*mass
   write(*,*) 'Beobachtungszeitraum (z.B.100) = '
   read(*,*) Viewtime
   write(*,*) 'Delta t (z.B. 0.1)   = '
   read(*,*) h
   
   gnudat='planet_plot.dat'
   gnucmd='planet_plot.cmd'
   
   print *,'Ich erzeuge zwei Files:'
   print *,trim(gnudat)
   print *,trim(gnucmd)
   print *,'Falls diese Files schon existieren, erfolgt eine Fehlermeldung.'
   open(unit=10,file=trim(gnudat),status='new')
   open(unit=11,file=trim(gnucmd),status='new')

   write(11,'(a)') 'set title "Planetenbewegung"'
   write(11,'(a,e15.5,a,e15.5,a)') 'set xrange [',x1min,':',x1max,']'
   write(11,'(a,e15.5,a,e15.5,a)') 'set yrange [',x2min,':',x2max,']'
   write(11,*)'set size ratio ',abs((x2max-x2min)/(x1max-x1min))
   write(11,'(a)') 'plot "'//trim(gnudat)//'" with points pt 7 '
   write(11,'(a)') 'pause -1'
   write(11,'(a)') 'quit'
   close(11)
   
   
!  Kraftzentrum
     write(10,*)0., 0.
   do i=1,21
     write(10,*)sin(i*pi/10)*x1max/40,cos(i*pi/10)*x1max/40
     write(10,*)sin(i*pi/10)*x1max/60,cos(i*pi/10)*x1max/60
   end do
   Nstep=int(Viewtime/h)
   do i=0,Nstep-1
     time=(i+1)*h
     call derivs(time,yvar,yder)
     call rk4(yvar,yder,nvar,time,h,outyvar)
     if ((outyvar(1) >= x1min).and.(outyvar(1) <= x1max).and.&
       & (outyvar(3) >= x2min).and.(outyvar(3) <= x2max)) then
       write(10,*)outyvar(1),outyvar(3)
     end if
     yvar=outyvar
   end do
   close(10)

!  Wenn das Unterprogramm "system" nicht verfuegbar ist, so kommentieren
!  Sie die folgende Zeile aus. Dann muessen Sie die Grafik selbst mittels
!  "gnuplot planet_plot.cmd" aufreufen.

   call system(trim('gnuplot '//trim(gnucmd)))
   print *,'Berechnung und Darstellung sind abgeschlossen.'

end

SUBROUTINE derivs(x,y,dydx)
!  Hier werden die Differentialgleichungen definiert
   use planet_parameters
   implicit none
  
   real, intent(in)::x
   real, dimension(1:nvar), intent(in)::y
   real, dimension(1:nvar), intent(out)::dydx
   real:: r,rd

   r=Sqrt(y(1)*y(1)+y(3)*y(3))
   rd=-cpot/(r*r*r)
   dydx(1)=y(2)/mass
   dydx(2)=-rd*y(1)
   dydx(3)=y(4)/mass
   dydx(4)=-rd*y(3)

END SUBROUTINE derivs

SUBROUTINE rk4(y,dydx,n, x,h,yout)

!  Programme, die RK4 verwenden, muessen eine
!  PROCEDURE derivs(x:real y:glnarray VAR dydx:glnarray)
!  bereitstellen, die die Ableitungen dy/dx (x) liefert, wenn x und y
!  gegeben ist.´
!  vereinbaren, wobei nvar die Anzahl (Dimension) der Variablen y ist.
!  RK4 ist modifiziert nach der urspruenglichen Vorlage von Press et. al. (PrFl99)}

   implicit none
  
   integer, parameter::nvarmax=200
   real, dimension(1:n), intent(in)::y,dydx
   real, dimension(1:n), intent(out)::yout
   integer, intent(in):: n
   real, intent(in)::x,h
   real :: xh,hh,h6
   real, dimension(1:nvarmax)::dym,dyt,yt
  
   if(n.gt.nvarmax) then
     write(*,*) 'SUBROUTINE rk4: internal maximum number of parameters exceeded. Stopped.'
     stop
   end if

   hh = h*0.5
   h6 = h/6.0
   xh = x+hh

   yt(1:n) = y(1:n)+hh*dydx(1:n)
   call derivs(xh,yt,dyt)
   yt(1:n) = y(1:n)+hh*dyt(1:n)
   call derivs(xh,yt,dym)
   yt(1:n) = y(1:n)+h*dym(1:n)
   dym(1:n) = dyt(1:n)+dym(1:n)
   call derivs(x+h,yt,dyt)
   yout(1:n) = y(1:n)+h6*(dydx(1:n)+dyt(1:n)+2.0*dym(1:n))
    
END SUBROUTINE rk4


