! This is a general purpose module which holds information on planetary
! bodies. In addition it contains a functions to compute the gravitational
! force between two bodies, and to check for collisions.
module gravity

  use constants ! use the first module in this one.
  implicit none ! Everything in the module will be given a specific type
                ! i.e. integer, real, logical or character.
  
  private ! specifies that everything in this module is hidden to the outside
          ! world unless it is explicitly declared with a public statement.

  ! Define a 2D array which we will use to store positions of our bodies.
  ! We make this an `allocatable` array, i.e. we don't specify how many
  ! elements it contains until later.
  real(kind=dp),allocatable,dimension(:,:) :: gravity_pos

  ! Another array for the velocities
  real(kind=dp),allocatable,dimension(:,:) :: gravity_vel

  ! A 1D array to store the masses
  real(kind=dp),allocatable,dimension(:) :: gravity_mass

  ! And another one to store the radii
  real(kind=dp),allocatable,dimension(:) :: gravity_rad

  ! A single variable to hold the elapsed time
  real(kind=dp) :: time

  ! Some integer variables which we will set later
  integer :: n  ! The number of masses
  integer :: d  ! The number of dimensions

  ! The gravitational constant - from Google is G = 6.67300E-11_dp.
  ! However we cannot work in S.I. units. If we carry around
  ! very big or very small numbers we will get rounding errors.
  ! Instead choose units such that

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

  ! In this system G  = 8.992293173e6_dp
  ! We label this a parameter so that it can't be changed
  ! anywhere else in the code. We append _dp to the value to
  ! specify it as being of kind dp, which we defined above.
  real(kind=dp),parameter ::  G =  8.992293173e6_dp

  ! state the stuff we need to be read/set from outside of this module
  public :: gravity_pos
  public :: gravity_vel
  public :: gravity_mass
  public :: gravity_rad

  ! state the functions subroutines that need to be useable from outside
  ! of this module
  public :: gravity_allocate
  public :: gravity_deallocate
  public :: gravity_force
  public :: gravity_check
  
  contains ! This module contains sub programs!
       
    ! Such as this function. Given the numbers associated with two bodies
    ! this returns the force between them.
    function gravity_force(i,j)

      implicit none  ! need this in every subprogram

      ! i and j are inputs to this routine
      integer,intent(in) :: i,j

      ! the return value of this function is an array holding the 
      ! vector force between the two bodies.
      real(kind=dp),dimension(d) :: gravity_force

      ! a variable for storing the distance between the two centres 
      ! of mass
      real(kind=dp) :: distance_squared  

      ! a variable for storing the ditance between the two centres
      real(kind=dp) :: distance

      ! an array of dimension d to hold the vector pointing from 
      ! the c.o.m. of mass i to the c.o.m. of mass j
      real(kind=dp),dimension(d) :: vector_rij

      ! Now we calculate stuff - first find the vector that we 
      ! have just made room for in memory.
      vector_rij(:) = gravity_pos(:,j) - gravity_pos(:,i)

      ! N.B. (:,j) means use the whole row j from the 2D array
      ! and (:) means use the whole 1D array. i.e. we can do
      ! vector algebra!

      ! find the magnitude-squared of this vector
      distance_squared = dot_product(vector_rij(:),vector_rij(:))
      ! and the magitude
      distance = sqrt(distance_squared)

      ! Note that the functions dot_product and sqrt are intrinsic to
      ! Fortran, we don't have to write them ourselves.

      ! force on j due to i - this sets the return value of the
      ! function using Newtons law of Gravitation
      gravity_force(:) = -G * vector_rij(:)*gravity_mass(i)*gravity_mass(j)/(distance*distance_squared) 

      return !finished!

    end function gravity_force

    ! Now lets have a subroutine. A subroutine is like a function
    ! but with no return value. i.e. thing = gravity_check()
    ! will return an error. Subroutines must be invoked by calling
    ! them, i.e. call gravity_check

    subroutine gravity_check(i,j)

      ! subroutine checks if bodies i and j have collided.
      implicit none

      ! i and j are inputs to this routine
      integer,intent(in) :: i,j

      ! a variable for storing the distance between the two centres 
      ! of mass
      real(kind=dp) :: distance_squared  

      ! an array of dimension d to hold the vector pointing from 
      ! the c.o.m. of mass i to the c.o.m. of mass j
      real(kind=dp),dimension(d) :: vector_rij

      ! set up a pretty format

10 format('Body ',I2,' has collided with body ',I2)

      ! Now we calculate stuff - first find the vector that we 
      ! have just made room for in memory.
      vector_rij(:) = gravity_pos(:,j) - gravity_pos(:,i)

      ! N.B. (:,j) means use the whole row j from the 2D array
      ! and (:) means use the whole 1D array. i.e. we can do
      ! vector algebra!

      ! find the magnitude-squared of this vector
      distance_squared = dot_product(vector_rij(:),vector_rij(:))

      ! compate this to the sum of the radii squared - note the 
      ! use of brackets
      if ( distance_squared < (gravity_rad(i)+gravity_rad(j))**2 ) then
         ! the two bodies have collided - scream to screen
         write(*,*)'Error in gravity module!'
         write(*,10)i,j
         write(*,*)'KABLAMO!'
         call gravity_deallocate() ! tidy up memory
         stop ! exit the program
      end if

      ! note we have avoided using sqrt() - taking the square root
      ! costs a lot of cpu cycles.

      return ! done

      ! note that some of the code in this subroutine is identical
      ! to that in the above function. We could save on coding
      ! by writing another function which computed distance_squared.

    end subroutine
    
    
    subroutine gravity_allocate(my_d,my_n)

      implicit none  ! need this in every subprogram

      ! intent(in) indicates that my_d and my_n are inputs to this
      ! subroutine.
      integer,intent(in) :: my_d,my_n

      ! we need to declare an integer to use as an error flag
      integer :: ierr

      ! set the values of n and d stored in this module to be
      ! equal to these inputs. Note we don't need ; on the end
      ! like we do in Maple.
      d = my_d
      n = my_n

      ! we can now allocate memory space in main memory for the arrays
      ! we declared at the top of this module. In each case we check
      ! that the allocation worked. If it didn't then the integer ierr
      ! will not be 0. ( /= means 'not equal to' )
      allocate(gravity_pos(1:d,1:n),stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error allocating gravity_pos'
         stop
      end if

      allocate(gravity_vel(1:d,1:n),stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error allocating gravity_vel'
         stop
      end if

      allocate(gravity_rad(1:n),stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error allocating gravity_rad'
         stop
      end if

      allocate(gravity_mass(1:n),stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error allocating gravity_mass'
         stop
      end if

      return ! done

    end subroutine gravity_allocate
    
    subroutine gravity_deallocate()

      ! release the memory used by this module

      ! again, need an error flag
      integer :: ierr

      deallocate(gravity_pos,stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error deallocating gravity_pos'
         stop
      end if
      deallocate(gravity_vel,stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error deallocating gravity_vel'
         stop
      end if

      deallocate(gravity_rad,stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error deallocating gravity_rad'
         stop
      end if

      deallocate(gravity_mass,stat=ierr)
      if (ierr/=0) then
         write(*,*)'Error deallocating gravity_mass'
         stop
      end if

      return ! done

    end subroutine gravity_deallocate
    

end module gravity
