!-------------------------------------------------------------------------
! ellisse.f90
! legge le posizioni della traiettoria di un pianet orbitante
! dalla 3a e 4a colonna del file posizioni.txt prodotto da "Keplero" in Java
! e, nel caso semplice in cui si parta dall'afelio o dal perielio
! (e uno dei fuochi, com'e' per costruzione in "Keplero", e' all'origine),
! se si forniscono le coordinate di afelio e perielio calcola 
! le deviazioni relative dall'ellisse ideale
!-------------------------------------------------------------------------
program ellisse
implicit none

integer,parameter::realkind=selected_real_kind(14)  !parametro per lavorare
                                                    !con reali con almeno 14
                                                    !cifre dec. significative  
real(kind=realkind)            :: x_A,x_B,x_F1,x_F2=0,y_F1=0,y_F2=0,a
real(kind=realkind)            :: dummy,x,y,dist,error
integer                        :: i, j, npoints

open(unit=1,file='traiettoria.txt')
open(unit=2,file='deviazioni-ellisse.txt')

print*," valido per ellisse con asse maggiore lungo x e un fuoco all'origine"
print*," A e B : intercette ellisse con asse x"

read(1,*)j,dummy,x,y
print*,' coordinate del primo punto letto nel file :',x,y
rewind(1)

print*, 'x_A >'
read*, x_A
print*, 'x_B >'
read*, x_B
print*, 'npoints >'
read*,npoints

a = abs(x_A - x_B) / 2 
x_F1 = x_A + x_B

print*,' F1 = (',x_F1,',',y_F1,')'
print*,' F2 = (',x_F2,',',y_F2,')'
print*,' a  = ',a

do i = 0,npoints       

   read(1,*)j,dummy,x,y
   dist = sqrt((x-x_F1)**2+(y-y_F1)**2)+sqrt((x-x_F2)**2+(y-y_F2)**2)
   error = (dist-2*a)/(2*a)

   write(2,*)error

end do

close(1)
close(2)

end program ellisse