ground track satellite

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,epoch,Ma,E,Mm1,Mm,M,P,t,x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3,AoP
real*8 :: GMST,longitude,r,alpha,phiprime,w,latitude
call readTLE(sma,ecc,inc,Omega,Ma,epoch,Mm,P)

Mm = Mm*86400
t=epoch
do while ((t-epoch)<P)
    M = Ma+Mm*(t-epoch)
    call eccano(ecc,M,E)
    call XYZ(ecc,E,sma,x,y,z)
    call rota(x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3,inc,Omega,AoP)
    call coorGPS(x3,y3,z3,epoch,longitude,latitude)
    t=t+0.00208333
    write(*,*)longitude*rad2deg, latitude*rad2deg,(t-epoch)*1440,"minutes"
    
enddo

end program main
!---------------------------------------------------------------------------!
subroutine readTLE(sma,ecc,inc,Omega,Ma,epoch,Mm,P)
use constantes
implicit none
integer :: y,e1

real*8 :: d,sma,ecc,inc,Omega,Ma,epoch,AoP,Mm,Mm1,P
read(*,'(18x,i2,f12.8)') y,d
epoch=d
if(y==20) then 
    epoch=d
endif
if(y==21) then 
    epoch=d+366
endif
print*, "epoch (days): ", epoch

read(*,"(9x,f7.4,2x,f7.4,5x,i3,1x,f8.4,1x,f8.4,x,f17.14)")inc,Omega,e1,AoP,Ma,Mm1

inc = inc*deg2rad
write(*,"(a21,f7.4)")"inclinaison (deg) : ",inc*rad2deg
Omega = Omega*deg2rad
write(*,"(a41,f7.4)")"Longitude of the ascending node (deg) : ",Omega*rad2deg
ecc=e1*1E-7
write(*,"(a16,f9.7)")"Eccentricity : ",ecc
AoP = AoP*deg2rad
write(*,"(a29,f8.4)")"Argument of Perigee (deg) : ",AoP*rad2deg
Ma = Ma*deg2rad
write(*,"(a21,f8.4)")"Mean anomaly (deg): ",Ma*rad2deg
Mm=((Mm1*2*pi)/(3600*24))
P=1/Mm1
print*,"period (min):",P*1440
sma= (mu**(1.0/3.0))/(Mm**(2.0/3.0))
print*, "semi major axis (km):",sma

end subroutine readTLE
!---------------------------------------------------------------------------------!
subroutine eccano(ecc,M,E)
use constantes
implicit none
real*8, intent(in) :: M,ecc
real*8, intent(out) :: E
integer :: n

E=M
E=(E-((E-ecc*sin(E)-M)/(1-ecc*cos(E))))

end subroutine


!---------------------------------------------------------------------------------!

subroutine XYZ(ecc,E,sma,x,y,z)
implicit none
real*8, intent(in) :: ecc,E,sma
real*8, intent(out):: x,y,z

x= sma*(cos(E)-ecc)
y= sma*sqrt(1-ecc**2)*sin(E)
Z=0

end subroutine

!---------------------------------------------------------------------!

subroutine rota(x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3,inc,Omega,AoP)
use constantes
implicit none
real*8, intent(in)::x,y,z,inc,Omega,AoP
real*8, intent(out)::x1,x2,x3,y1,y2,y3,z1,z2,z3

x1 = cos(AoP)*x - sin(AoP)*y
y1 = sin(AoP)*x + cos(AoP)*y
z1 = z


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


end subroutine

!-----------------------------------------------------------------------------!

subroutine coorGPS(x3,y3,z3,epoch,longitude,latitude)
use constantes
implicit none
real*8, intent(in) :: epoch,x3,y3,z3
real*8, intent(out) :: longitude,latitude
real*8:: GMST,phiprime,alpha,r,w

GMST= 0.177834888 + (6.300387958*(epoch - 275))
alpha = atan2(y3,x3)
longitude = (alpha-GMST)

r=sqrt(x3**2+y3**2+z3**2)
phiprime = asin(z3/r)
w=(tan(phiprime))/(1-0.0033528**2)
latitude=atan(w)


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

do while(latitude>pi)
    latitude=latitude-2*pi
enddo

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

do while(longitude>pi)
    longitude=longitude-2*pi
enddo

end subroutine

Comments