Project final Code Gabriel, Mathilde, Alexis

module constantes
implicit none
save
real*8, parameter :: pi=3.14159274 ! pi
real*8, parameter :: mu=3.98600E5 ! GMearth (km^3/s^2)
real*8, parameter :: Re=6378.14 ! Earth radius (km)
real*8, parameter :: deg2rad=2.0*pi/360.0 ! deg -> rad
real*8, parameter :: rad2deg=180.0/pi ! rad -> deg
end module constantes

program main
use constantes
implicit none 
real*8 :: sma,ecc,inc,Omega,om,M,epoch,E,Ma,M0,Mm,t,x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3,AoP,Lati,Longi,P
call readTLE(sma,ecc,inc,Omega,M0,epoch,Mm,P)
Mm = Mm * 86400
t=epoch

do while ((t-epoch)<P)    

Ma=M0+Mm*(t-epoch)
call eccano(ecc,Ma,E)
call coordinate(x,y,z,sma,E,ecc)
call rotation(x,y,z,x3,y3,z3,inc,Omega,AoP)
call latilongi(x3,y3,z3,epoch,Lati,Longi,t)
write(*,*)Longi*rad2deg, Lati*rad2deg,(t-epoch)
t=t+0.00138888888
enddo

end program main

!==============================================================================================================
subroutine readTLE(sma,ecc,inc,Omega,M0,epoch,Mm,P)
use constantes
implicit none
integer ::y,e1
real*8 ::d,sma,ecc,inc,Omega,epoch,AoP,M0,Mm,Mm1,Ma,P

read(*,"(18x,i2,f12.8)")y,d
epoch=d
if (y==20) then
    epoch=d
endif
if (y==21) then 
    epoch=d+366.0
endif
if (y==22) then 
    epoch=d+366.0+365.0
endif

!print*, "The epoch is :",epoch

read(*,"(9x,f7.4,2x,f6.4,5x,i4,2x,f7.4,2x,f7.4,x,f17.14)")inc,Omega,e1,AoP,M0,Mm1
inc = inc*deg2rad
!write(*,"(a25,f7.4)")"inclination (radiant) : ",inc
Omega = Omega*deg2rad
!write(*,"(a49,f6.4)")"Longitude of the ascending node Ω (radiant) : ",Omega
ecc=e1*1E-7
!write(*,"(a50,f9.7)") "Eccentricity (radiant) (decimal point assumed) :  ",ecc
AoP = AoP*deg2rad
!write(*,"(a31,f7.4)")"Argument of Perigee (radiant) :",AoP
M0=M0*deg2rad
!write(*,"(a30,f7.4)")"Mean anomaly M (in radiant) :",M0
!print*, "The revolution per day is :",Mm1
P=1/Mm1
!print*,"The period is",P,"day"
Mm=((Mm1*2*pi)/86400)
!write(*,"(a38,f18.14)")"Mean motion (radiant per second) is :",Mm
sma= (mu**(1.0/3.0))/(Mm**(2.0/3.0))
!print*, "The semi major axis is:", sma

end subroutine readTLE

!=============================================================================================================================
subroutine eccano(ecc,Ma,E)
use constantes
implicit none
real*8, intent(in) :: Ma,ecc
real*8, intent(out) :: E
integer :: n
n=0
E=Ma
n=1

do while (n<=3)

E=(E-(((E-ecc*sin(E)-Ma)/(1-ecc*cos(E)))))
n=n+1

enddo

!print*,"The eccentric anomaly is :",E

end subroutine
!==============================================================================================================

subroutine coordinate(x,y,z,sma,E,ecc)
use constantes
implicit none
real*8, intent(in) :: sma,E,ecc 
real*8, intent(out) :: x,y,z
x = sma*(cos(E)-ecc)
y = sma*sqrt(1-ecc**2)*sin(E)
z = 0

!print*,"The x coordinate is",x 
!print*,"The coordinate y is",y

end subroutine

!======================================================================================
subroutine rotation(x,y,z,x3,y3,z3,inc,Omega,Aop)
use constantes
implicit none
real*8, intent(in) :: x,y,z,inc,Omega,Aop
real*8, intent(out) :: x3,y3,z3
real*8 :: x1,y1,z1,x2,y2,z2,norm
x1 = cos(AoP)*x-sin(AoP)*y
y1 = sin(AoP)*x+cos(AoP)*y
z1 = 0 

x2 = x1
y2 = cos(inc)*y1-sin(inc)*z1
z2 = sin(inc)*y1+cos(inc)*z1

x3 = cos(Omega)*x2-sin(Omega)*y2
y3 = sin(Omega)*x2+cos(Omega)*y2
z3 = z2

!print*,"The coordinate x3 is",x3
!print*,"The coordinate y3 is",y3
!print*,"The coordinate z3 is",z3


end subroutine

!========================================================================================================================
subroutine latilongi(x3,y3,z3,epoch,Lati,Longi,t)
use constantes
implicit none
real*8, intent(in) :: x3,y3,z3,epoch,t
real*8, intent(out) :: Lati,Longi
real*8 :: geoctrd,Lati2,alpha,w,GMST

geoctrd = sqrt(x3**2+y3**2+z3**2)
!print*,"Geo ctr distance",geoctrd
Lati2 = asin(z3/geoctrd)
!print*,"Lati prime",Lati2
alpha = atan2(y3,x3)
!print*,"alpha",alpha
w = tan(Lati2)/(1-(1/298.257)**2)
!print*,"w is",w
Lati = atan(w)
!print*,"Latitude is",Lati
GMST = 0.177834888 + 6.300387958*(t - 275) 
!print*,"Greenwich meridian sideral time",GMST
Longi = (alpha - GMST)

do while (Lati<-pi)

    Lati = Lati + 2*pi
enddo
do while (Lati>pi)

    Lati=Lati - 2*pi
enddo

do while (Longi<-pi)
    
    Longi = Longi + 2*pi
enddo
do while (Longi>pi)

    Longi = Longi - 2*pi
enddo

!print*,"The longitude coordinate of the satellite is",Longi*rad2deg
!print*,"The Latitude coordinate of the satellite is", Lati*rad2deg

end subroutine

Comments