! -*- mode: F90 ; mode: font-lock ; column-number-mode: true
!=============================================================================!
!                            E A R T H - M O O N                              !
!=============================================================================!
! Simple Fortran 90 program. Uses the modules contained in the files          !
! constants.f90 and gravity.f90.                                              !
!                                                                             !
! Written to illustrate the level of fortran knowledge requried for the hpc   !
! course. If you can understand most of this then you will be more than set!  !
!                                                                             !
!-----------------------------------------------------------------------------!
! Written by David Quigley, June 2005                                         !
!-----------------------------------------------------------------------------!
program earthmoon

  use constants   ! Use both the modules that we created
  use gravity
  implicit none

  ! an array to store accelerations. Note that in this program unit we
  ! are coding for the two body case only. Whereas in the general 
  ! purpose gravity module we allowed for any number of bodies.
  real(kind=dp),dimension(3,2) :: accel

  ! intermediate K values needed by the Runge-Kutta
  real(kind=dp),dimension(3,2) :: k0_pos,k1_pos,k2_pos,k3_pos,k4_pos
  real(kind=dp),dimension(3,2) :: k0_vel,k1_vel,k2_vel,k3_vel,k4_vel

  ! an error flag ( a variable to store error information )
  integer :: ierr

  ! some loop counters
  integer :: istep

  ! the calculation timestep and the evolved time
  real(kind=dp) :: delta_t , time

  ! some numbers we will need
  real(kind=dp),parameter :: third = 1.0_dp/3.0_dp
  real(kind=dp),parameter :: sixth = 1.0_dp/6.0_dp

  ! number of steps to perform
  integer(kind=dp) :: nsteps

  ! set up some pretty format statements 
  ! three columns or 15 characters with 8 after the decimal
10 format("Vector position of Earth = ",3E15.5,' at t = '3E15.5)
20 format("Vector position of Moon  = ",3E15.5,' at t = '3E15.5)

  ! Set up the gravity module, i.e. tell it we are using
  ! two bodies in three dimensions so it can allocate
  ! memory.
  call gravity_allocate(3,2)

  ! Reminder of units system

  ! 1 unit of length = radius of earth
  ! 1 unit of mass   = mass of earth
  ! 1 unit of time   = 28 days

  ! Set the masses - note how we set the value of an array element
  gravity_mass(1) =  1.00_dp        ! Mass of the Earth is 1 unit
  gravity_mass(2) =  0.012319641_dp ! Moon in Earth mass units

  ! Set the radii
  gravity_rad(1) = 1.00_dp        ! Radius of the Earth is 1 length unit
  gravity_rad(2) = 0.272400871_dp ! Radius of the Moon 

  ! set up some initial conditions
  gravity_pos(:,1) = 0.0_dp  ! Earth is at the origin

  ! put the moon at 405,500 km from the Earth along the x-axis
  gravity_pos(1,2)   = 63.5769273_dp ! 405,500 km in earth radii
  gravity_pos(2:3,2) = 0.0_dp        ! y and z components

  ! let's have the Earth initially at rest
  gravity_vel(:,1)   = 0.0_dp

  ! to create an orbit the moon should be moving in say the z
  ! direction only. 
  gravity_vel(1:2,2) = 0.0_dp ! x and y components

  ! If one unit of time is 28 days the moon velocity should be
  ! two * pi * radius or orbit / 28 days
  gravity_vel(3,2) = 2.0_dp * pi * 63.5769273_dp 

  ! A timestep of about 4 hours should be small enough for accuracy
  delta_t = 0.006_dp

  ! simulate for 500 steps
  nsteps = 500

  ! Main do loop. This says loop over values of istep from one to nsteps.
  ! All code between here and end do is executed on each pass through the
  ! loop.
  do istep = 1 , nsteps 
 
     !====================================================================
     !  Start of Runge-Kutta stuff - there's no new fortran in here.     !
     !  This is simple integration of two coupled 1st order differential !
     !  equations, just like in the MAPLE assignment. Ignore it if you   !
     !  want.                                                            !
     !====================================================================

     ! keep a copy of the current positions and velocities
     k0_pos(:,:) = gravity_pos(:,:)
     k0_vel(:,:) = gravity_vel(:,:)

     ! Acceleration of Earth due to moon - using the force function in 
     ! the useful_stuff module
     accel(:,1) =  gravity_force(2,1) / gravity_mass(1)
     ! acceleration of moon due to Earth
     accel(:,2) =  gravity_force(1,2) / gravity_mass(2)

     ! calculate k1 positions and velocities
     k1_pos(:,:) = delta_t * gravity_vel(:,:)
     k1_vel(:,:) = delta_t * accel(:,:)
     
     ! update 1
     gravity_pos(:,:) = k0_pos(:,:) + 0.5_dp * k1_pos(:,:)
     gravity_vel(:,:) = k0_vel(:,:) + 0.5_dp * k1_vel(:,:)

     ! Acceleration of Earth due to moon - using the force function in 
     ! the useful_stuff module
     accel(:,1) =  gravity_force(2,1) / gravity_mass(1)
     ! acceleration of moon due to Earth
     accel(:,2) =  gravity_force(1,2) / gravity_mass(2)
     
     ! calculate k2 positions and velocities
     k2_pos(:,:) = delta_t * gravity_vel(:,:)
     k2_vel(:,:) = delta_t * accel(:,:)

     ! update 2
     gravity_pos(:,:) = k0_pos(:,:) + 0.5_dp * k2_pos(:,:)
     gravity_vel(:,:) = k0_vel(:,:) + 0.5_dp * k2_vel(:,:)

     ! Acceleration of Earth due to moon - using the force function in 
     ! the useful_stuff module
     accel(:,1) =  gravity_force(2,1) / gravity_mass(1)
     ! acceleration of moon due to Earth
     accel(:,2) =  gravity_force(1,2) / gravity_mass(2)

     ! calculate k3 positions and velocities
     k3_pos(:,:) = delta_t * gravity_vel(:,:)
     k3_vel(:,:) = delta_t * accel(:,:)

     ! update 3
     gravity_pos(:,:) = k0_pos(:,:) + k3_pos(:,:)
     gravity_vel(:,:) = k0_vel(:,:) + k3_vel(:,:)

     ! Acceleration of Earth due to moon - using the gravity_force function in 
     ! the useful_stuff module
     accel(:,1) =  gravity_force(2,1) / gravity_mass(1)
     ! acceleration of moon due to Earth
     accel(:,2) =  gravity_force(1,2) / gravity_mass(2)

     ! calculate k4 positions and velocities
     k4_pos(:,:) = delta_t * gravity_vel(:,:)
     k4_vel(:,:) = delta_t * accel(:,:)

     ! final update - note the & to continue onto a second line
     gravity_pos(:,:) = k0_pos(:,:) + sixth*k1_pos(:,:) + third*k2_pos(:,:) &
                                    + third*k3_pos(:,:) + sixth*k4_pos(:,:) 
     gravity_vel(:,:) = k0_vel(:,:) + sixth*k1_vel(:,:) + third*k2_vel(:,:) &
                                    + third*k3_vel(:,:) + sixth*k4_vel(:,:) 

     !=========================================================
     !  End of Runge-Kutta stuff                            
     !=========================================================

     ! increment time
     time = time + delta_t
     
     ! check for the end of the world
     call gravity_check(1,2)

     ! write some information to the screen. This can be redirected to a file
     ! and `grepped'. These write commands use the pretty formats we set up earlier
     write(*,10)gravity_pos(:,1),time
     write(*,20)gravity_pos(:,2),time

  end do

  ! release the memory used by the gravity module
  call gravity_deallocate()



end program earthmoon
