module gravitazione
implicit none

interface operator (.x.)
    module procedure vec
end interface operator (.x.)

integer,parameter :: kr=selected_real_kind(15), nbody=3
real(kind=kr),parameter :: G=6.673E-11_kr  ! costante di gravitazione universale in unita' SI
real(kind=kr),dimension(nbody) :: mass

contains


 function vec(x,y) result(r)
   real(kind=kr),intent(in),dimension(3) :: x,y
   real(kind=kr),dimension(3) :: r
   r(1) = x(2)*y(3)-x(3)*y(2)
   r(2) = x(3)*y(1)-x(1)*y(3)
   r(3) = x(1)*y(2)-x(2)*y(1)
 end function vec


  subroutine interazione(pos,f,epot)
   real(kind=kr), intent(in), dimension(:,:) :: pos
   real(kind=kr), intent(out) :: epot
   real(kind=kr), intent(out), dimension(:,:) :: f
   real(kind=kr), dimension(size(pos,1)) :: posij
   real(kind=kr) :: rij
   integer ::i,j
   epot = 0
   f    = 0
   do i=1,nbody
      do j=1,nbody
         if( i==j ) cycle
         posij(:)  = pos(:,j)-pos(:,i)
         rij=sqrt( dot_product(posij,posij) )
         epot   = epot   -(G*mass(i)*mass(j))/rij
         f(:,i) = f(:,i) +(G*mass(i)*mass(j)) * posij(:)/rij**3
      end do
   end do
   
   epot = epot/2    ! divisione per 2 perche' tutti i contributi 
                    ! sono stati contati due volte: (i,j) e (j,i) 
end subroutine interazione
end module gravitazione


program keplero
use gravitazione, massa=> mass
implicit none

real(kind=kr)    :: dt,ekin,epot,dist,ang,ang_prec
real(kind=kr),dimension(3):: L,vcm, posTS
real(kind=kr),dimension(3,nbody):: pos,vel,f
integer :: nstep,it,i, igiro


write(unit=*,fmt="(a)",advance="no")"delta t : " 
read*,dt                                             
write(unit=*,fmt="(a)",advance="no")"n.step: "
read*,nstep                                       
write(unit=*,fmt="(a)",advance="no")"masse: "   
read*,massa
write(unit=*,fmt="(a)",advance="no")"array posizione iniziale: "
read*,pos
write(unit=*,fmt="(a)",advance="no")"array velocità iniziale: "
read*,vel

it=0 
vcm = matmul(vel,massa)/sum(massa)  !vel. centro di massa
! vcm_i = sum_{j=1,nbody} m_j vel_{i,j} / sum_{j=1,nbody} m_j
vel = vel - spread(vcm,2,nbody)     !va sottratta all'inizio da tutte le vel. dei 3 corpi per mettersi nel sistema del CM
! l'array vcm(1:3) viene replicato nbody volte lungo la dimensione (2) di vel(1:3,1:nbody)

write(unit=1,fmt=*)"# it,it*dt,pos,vel"
write(unit=1,fmt=*)it,it*dt,pos,vel
dist = sqrt(sum((pos(:,2)-pos(:,1))**2))
write(unit=3,fmt=*)"# it,it*dt,pos(:,2)-pos(:,1),dist"
write(unit=3,fmt=*)it,it*dt,pos(:,2)-pos(:,1),dist

call interazione(pos,f,epot)

ekin=0
do i=1, nbody
   ekin=ekin + 0.5*massa(i)*(dot_product(vel(:,i),vel(:,i)))
end do

L=0
do i=1, nbody
   L=L+(pos(:,i).x.(massa(i)*vel(:,i)))
end do

write(unit=2,fmt=*)"# it,dt*it,ekin,epot,ekin+epot,L"
write(unit=2,fmt=*)it,dt*it,ekin,epot,ekin+epot,L
igiro = 0
do it = 1,nstep
    pos = pos + vel * dt + 0.5* f/spread(massa,1,3) * dt**2
    vel = vel + 0.5 * dt * f/spread(massa,1,3) 
     call interazione(pos,f,epot)
    vel = vel + 0.5 * dt * f/spread(massa,1,3)
    
    write(unit=1,fmt=*)it,it*dt,pos,vel
    dist = sqrt(sum((pos(:,2)-pos(:,1))**2))
    write(unit=3,fmt=*)it,it*dt,pos(:,2)-pos(:,1),dist


    ekin=0
    do i=1, nbody
    ekin=ekin + 0.5*massa(i)*(dot_product(vel(:,i),vel(:,i)))
    end do
    
    L=0
    do i=1, nbody
    L=L+(pos(:,i).x.(massa(i)*vel(:,i)))
    end do
    
    write(unit=2,fmt=*)it,it*dt,ekin,epot,ekin+epot,L
end do


end program keplero

!dati di input:
!46200
!70000
!1.989e30 5.89e24 1898.6e24
!0 0 0
!152.1e9 0 0
!816.62e9 0 0
!0 0 0
!0 26.29e3 0
!0 12.44e3  0