code

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,X3,Y3,Z3,latitude,X2,Y2,Z2,longitude,X0,Y0,Z0,X1,Y1,Z1,AoP
call readTLE(sma,ecc,inc,Omega,M0,epoch,Mm) 
t=326
Ma=M0+Mm*(t-epoch)
call eccano(ecc,Ma,E)
call xyz(X0,Y0,Z0,sma,E,ecc)
call xyz1(X0,Y0,Z0,X1,Y1,Z1,AoP)
call xyz3(inc,X1,Y1,Z1,X2,Y2,Z2)
call xyz2(Omega,X2,Y2,Z2,X3,Y3,Z3)
call latitudelongitude(X3,Y3,Z3,epoch,latitude,longitude,sma)


end program main
!==============================================================================================================
subroutine readTLE(sma,ecc,inc,Omega,M0,epoch,Mm)
use constantes
implicit none
integer ::y,e1
real*8 ::d,sma,ecc,inc,Omega,epoch,AoP,M0,Mm,Mm1,Ma
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
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 
write(*,"(a25,f7.4)")"inclination (degrees) : ",inc
write(*,"(a49,f7.4)")"Longitude of the ascending node Ω (degrees) : ",Omega
ecc=e1*1E-7
write(*,"(a41,f9.7)") "Eccentricity (decimal point assumed) : ",ecc
write(*,"(a22,f7.4)")"Argument of Perigee :",AoP
M0=M0*deg2rad
write(*,"(a30,f7.4)")"Mean anomaly M (in radiant) :",M0
print*, "The revolution per day is :",Mm1
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 
!=============================================================================================================================
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
print*,"E number 0 =",E
n=1
do while (n<=3)
E=(E-(((E-ecc*sin(E)-Ma)/(1-ecc*cos(E)))))
print*, "E number", n, "=", E
n=n+1
enddo
print*,"The eccentric anomaly is :",E

end subroutine
!=============================================================================================================================
subroutine xyz(X0,Y0,Z0,sma,E,ecc)
implicit none 
real*8, intent(in) :: sma,E,ecc
real*8, intent(out) :: X0,Y0,Z0
X0=sma*(cos(E)-ecc)
Y0=sma*sqrt(1-ecc**2)*sin(E)
Z0=0 
print*,"x is",X0
print*,"y is",Y0
print*,"z is",Z0

end subroutine
!=============================================================================================================================
subroutine xyz1(X0,Y0,Z0,X1,Y1,Z1,AoP)
implicit none 
real*8,intent(in) :: X0,Y0,Z0,AoP
real*8,intent(out) :: X1,Y1,Z1

X1=cos(AoP)*X0-sin(AoP)*Y0
Y1=sin(Aop)*X0+cos(AoP)*Y0
Z1=0
print*,"The position X1, Y1 and Z1 are",X1,Y1,Z1

end subroutine 
!=============================================================================================================================
subroutine xyz3(inc,X1,Y1,Z1,X2,Y2,Z2)
implicit none
real*8,intent(in) :: X1,Y1,Z1,inc
real*8,intent(out) :: X2,Y2,Z2
X2=X1
Y2=cos(inc)*Y1-sin(inc)*Z1
Z2=sin(inc)*Y1-cos(inc)*Z1

print*,"The position X2, Y2 and Z2 are", X2,Y2,Z2


end subroutine

!=============================================================================================================================
subroutine xyz2(Omega,X2,Y2,Z2,X3,Y3,Z3)
implicit none
real*8,intent(in) :: Omega,X2,Y2,Z2
real*8,intent(out) :: X3,Y3,Z3
X3=cos(Omega)*X2-sin(Omega)*Y2
Y3=sin(Omega)*X2+cos(Omega)*Y2
Z3=Z2
print*,"The position X3, Y3 and Z3 are",X3,Y3,Z3

end subroutine
!=============================================================================================================================
subroutine latitudelongitude(X3,Y3,Z3,epoch,latitude,longitude,sma)
use constantes
implicit none 
real*8, intent(in) :: X3,Y3,Z3,epoch
real*8, intent(out) :: latitude,longitude 
real*8 :: r,geocentriclat,alpha,GMST,w,t,sma
r=sma
print*,"r is",r 
geocentriclat=asin(Z3/r)
print*,"geocentric latitude is",geocentriclat
alpha=atan2(Y3,X3)
print*,"alpha is",alpha
w=tan(geocentriclat)
print*,"w",w 
latitude=atan(w)
GMST=0.177834888+6.300387958*(t-275)
print*,"GMST",GMST
longitude=(alpha-GMST)
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
    

print*,"Longitude",longitude*deg2rad
print*,"Latitude",latitude*deg2rad
end do
end subroutine

Comments