!-------------------------------------------------------------------------
! keplero.f90
!
! Integra l'equazione del modo per un pianeta orbitante attorno al sole
! (considerato fisso) usando l'algoritmo di Verlet;
! calcola posizioni, velocita', energie, momento angolare 
! e area spazzata dal pianeta orbitante  
!
!-------------------------------------------------------------------------
program keplero
implicit none

integer,parameter::realkind=selected_real_kind(14)  !parametro per lavorare
                                                    !con reali con almeno 14
                                                    !cifre dec. significative  
type::vett2d                   ! def tipo derivato in fortran. Equivalente 
                               ! ad un tipo RECORD in Pascal
                               ! puo' anche essere sostituito da un ARRAY 
      real(kind=realkind)::x,y
end type

! unita' astronomiche: 1 AU = 1.496 * 10^11 m
!                      1 yr = 3.15 * 10^7 s
!                      GM = (4 pi^2 a^3)/T^2 = 4 pi^2 AU^3/yr^2

!real(kind=realkind), parameter :: pi = 3.141569265358979323846
!real(kind=realkind), parameter :: gm = 4.0*pi*pi 

real(kind=realkind), parameter :: g= 6.67E-11 ! m^3/(kg s) 
real(kind=realkind)            :: massa
real(kind=realkind)            :: dt,ekin,epot,r,r_old,k,Lz,area,deltar
type(vett2d)                   :: pos,vel,a,pos_old
integer                        :: nstep,it

open(unit=1,file='pos.out')
open(unit=2,file='area.out')
open(unit=3,file='en.out')

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")"massa che genera il campo (sole = &
&1.99E30): "
read*,massa
write(unit=*,fmt="(a)",advance="no")"posizione x (terra, afelio = 1.521E11): "
read*, pos%x
write(unit=*,fmt="(a)",advance="no")"posizione y (terra = 0): "
read*, pos%y
write(unit=*,fmt="(a)",advance="no")"velocita' x (terra = 0): "
read*, vel%x
write(unit=*,fmt="(a)",advance="no")"velocita' y (terra, vel_min = 2.929E4): "
read*, vel%y

it=0                     ! step 0 : valori iniziali

r = (pos%x**2 + pos%y**2)**0.5
a%x    = - g * massa * pos%x / r**3     ! legge di gravitazione
a%y    = - g * massa * pos%y / r**3     !       "
ekin =  0.5 * (vel%x**2+vel%y**2)       ! en. cinetica iniziale / massa pianeta
epot = - g*massa/r                      ! en. potenziale iniziale / massa pianeta

write(unit=1,fmt=*)it,it*dt,pos%x,pos%y,vel%x,vel%y
write(unit=3,fmt=*)it,it*dt,ekin,epot,ekin+epot

do it = 1,nstep

   pos_old%x = pos%x
   pos_old%y = pos%y
   r_old     = r

   pos%x = pos%x + vel%x * dt + 0.5* a%x * dt**2
   pos%y = pos%y + vel%y * dt + 0.5* a%y * dt**2

   vel%x = vel%x + 0.5 * dt * a%x
   vel%y = vel%y + 0.5 * dt * a%y

! qui le posizioni sono gia' quelle allo step t+dt mentre per le velocita'
! manca ancora il termine con l' accelerazione a t+dt

   r = (pos%x**2 + pos%y**2)**0.5

   deltar=((pos%x-pos_old%x)**2+(pos%y-pos_old%y)**2)**0.5
   k     = (r+r_old+deltar) * 0.5  ! semiperimetro del triangolino
         ! spazzato dal vettore posizione del pianeta in [t,t+dt]
   area  = (k*(k-r_old)*(k-r)*(k-deltar))**0.5   ! area del triangolino
         ! suddetto calcolata con la formula di Erone; 
         ! moltiplicata per 2 e divisa per dt, 
         ! dev'essere uguale al momento angolare 

   a%x    = - g * massa * pos%x / r**3
   a%y    = - g * massa * pos%y / r**3
   epot =  -g*massa/r

! ora posso completare l'aggiornamento delle velocita' 
! con il termine con l' accelerazione a t+dt

   vel%x = vel%x + 0.5 * dt * a%x
   vel%y = vel%y + 0.5 * dt * a%y

!adesso anche le velocita' sono state aggiornate al valore a t+dt
!e calcolo en. cinetica e componente z del momento angolare,
!a meno della massa del pianeta orbitante
  
   ekin =  0.5* (vel%x**2+vel%y**2)
   Lz   =  pos%x*vel%y-pos%y*vel%x   !L_z = (r x mv)_z [qui non c'e' m]

   write(unit=1,fmt=*)it,it*dt,pos%x,pos%y
   write(unit=2,fmt=*)it,it*dt,area,Lz
   write(unit=3,fmt=*)it,it*dt,ekin,epot,ekin+epot

end do
end program keplero