program main
  !
  ! Rozwiązywanie zagadnienia Arenstorfa 3 ciał metodą RK
  ! -------------------------------------------------------
  ! 
  ! /user> gfortran drode.dr ode.f90
  ! /user> ./a.out > out
  ! /user> gnuplot
  !        gnuplot>  plot 'out' u 2:3 with lines
  !
  ! (Otrzymany rysunek jest orbitą satelity)
  !
  ! -------------------------------------------------------
  ! Inny wariant obliczeń:
  ! Zamienić warunki początkowe 1 <--> 2 podane poniżej
  ! powtórzyć obliczenia
  ! -------------------------------------------------------
  ! ZADANIA
  ! 1. Wykonać rysunki (t,y1), (t,y2), (t,y3), (t,y4), (y1,y2)
  ! 2. Sprawdzić wyniki dla różnej liczby zadanych kroków
  !    całkowania
  ! 3. Powtórzyć obliczenia dla różnych dokładności:
  !    parametry abser, relerr
  !    Przeczytać komentarze w pliku ode.f90
  !
  implicit none

  integer, parameter :: neqn = 4

  real ( kind = 8 ) abserr
  external faren
  integer i
  integer iflag
  integer iwork(5)
  real ( kind = 8 ) relerr
  integer, parameter :: step_num = 100
  real ( kind = 8 ) t
  real ( kind = 8 ) tout
  real ( kind = 8 ) work(100+21*neqn)
  real ( kind = 8 ) y(neqn)
  real (kind=8) tend

  write ( *, '(a)' ) '# Arenstorf orbits'
  write ( *, '(a)' ) '#    T           Y(1)         Y(2)        Y(3)       Y(4)'

  abserr = 0.00001D+00
  relerr = 0.00001D+00

  iflag = 1

  t = 0.0D+00

  ! Warunki początkowe 1
  !
  y(1)=.87978d0
  y(2)=0.0d0
  y(3)=0.0d0
  y(4)=-0.3797d0
  tend=19.14045706162071d0

  ! Warunki początkowe 2
  !
  y(1)=0.994d0
  y(2)=0.0d0
  y(3)=0.0d0
  y(4)=-2.00158510637908252240537862224d0
  tend=17.0652165601579625588917206249d0

  write ( *, '(2x,f8.4,2x,4g14.6)' ) t, y(1:neqn)

  do i = 1, step_num

     tout = real ( i, kind = 8 ) * tend / real ( step_num, kind = 8 )

     call ode ( faren, neqn, y, t, tout, relerr, abserr, iflag, work, iwork )

     if ( iflag /= 2 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' )    '* Fatal error!'
        write ( *, '(a,i8)' ) '  DR returned IFLAG = ', iflag
        exit
     end if

     write ( *, '(2x,f8.4,2x,4g14.6)' ) t, y(1:neqn)

  end do

  return
end program main
!
!
subroutine faren(x,y,f)

  ! arenstorf orbit

  implicit real*8 (a-h,o-z)
  dimension y(4),f(4)

  amu=0.012277471d0
  amup=1.d0-amu

  f(1)=y(3)
  f(2)=y(4)
  r1=(y(1)+amu)**2+y(2)**2
  r1=r1*sqrt(r1)
  r2=(y(1)-amup)**2+y(2)**2
  r2=r2*sqrt(r2)
  f(3)=y(1)+2*y(4)-amup*(y(1)+amu)/r1-amu*(y(1)-amup)/r2
  f(4)=y(2)-2*y(3)-amup*y(2)/r1-amu*y(2)/r2

  return
end subroutine faren
