%% Траектория спутников GPS (GPS Satellite Ground Tracks)
%%
PRN=10; % Номер спутника
filename='zimm1500.07n'; % Navigation message file's name
[week sec nav]=prnnav(PRN,filename); % Получение навигационного сообщения от спутника
%% Рассчет эфемериды для полного дня
mu=3.986005*10^14; % WGS-84 value of the earths universal gravitational parameter
OMe=7.2921151467*10^-5; % WGS-84 value of the earths roration rate (rad/s)
i=1;
TO=linspace(1,86400,100);

for j=1:length(TO)
    
Crs=nav(i,5); % Amplitude of the cosine harmonic correction term to orbit radius(meter)
dn=nav(i,6); % Correction for mean motion
M0=nav(i,7); % Mean Anomaly at refernce time
Cuc=nav(i,8); % Amplitude of the sine harmonic correction term to the argument of lattitude(Rad)
e=nav(i,9); % Orbit eccentricity
Cus=nav(i,10); % Amplitude of the cosine harmonic correction term to the argument of lattitude(Rad)
toe=nav(i,12); % Reference time of ephemeris parameters
Cic=nav(i,13); % Amplitude of the sine harmonic correction term to the angle of inclination(Rad)
OM0=nav(i,14); % Longitude of Ascending Node
Cis=nav(i,15); % Amplitude of the cosine harmonic correction term to the angle of inclination(Rad)
i0=nav(i,16); % Inlination angle at reference time
Crc=nav(i,17); % Amplitude of the sine harmonic correction term to orbit radius(meter)
om=nav(i,18); % Argument of perigee (semicicrcle)
OM=nav(i,19); % Rate of right ascension
IDOT=nav(i,20); % Rate of inlination angle

A=(nav(i,11))^2; % Semi-major orbit axis (meter)
n0=sqrt(mu/(A^3)); % Compute mean motion (rad/s)
tk=TO(j)-toe; % Time from ephemeris reference epoch
n=n0+dn; % Corrected mean motion
Mk=M0+n*tk; % Mean anomaly
Ek=fzero(@(x) Kepler_Eq(x,Mk,e),Mk); % Eccentric anomaly Ek (rad)
nuk=atan((sqrt(1-e^2)*sin(Ek))/((cos(Ek)-e))); % True anomaly as a function of the eccentric anomaly
if nuk<0
    nuk=nuk+2*pi;
end;
Fk=nuk+om; % Argument of lattitude

duk=Cus*sin(2*Fk)+Cuc*cos(2*Fk); % Argument of lattitude correction
drk=Crs*sin(2*Fk)+Crc*cos(2*Fk); % Argument of radius correction
dik=Cis*sin(2*Fk)+Cic*cos(2*Fk); % Argument of inclination correction

uk=Fk+duk; % Corrcted argument of lattitude
rk=A*(1-e*cos(Ek))+drk; % Corrected radius
ik=i0+dik+IDOT*tk; % Corrected inclination

xkshtr=rk*cos(uk); % Satellite position in orbital plane
ykshtr=rk*sin(uk);

OMk=OM0+(OM-OMe)*tk-OMe*toe; % Corrected longitude of ascending node (rad)

xk(j)=xkshtr*cos(OMk)-ykshtr*cos(ik)*sin(OMk); % Satellite Position in ECEF
yk(j)=xkshtr*sin(OMk)+ykshtr*cos(ik)*cos(OMk);
zk(j)=ykshtr*sin(ik);

[longitudetek latitudetek R] = cart2sph(xk(j),yk(j),zk(j)); % Преобразование x, y, z (errh-fixed) в широту и долготу
longitude(j)=longitudetek*180/pi; % Долгота в градусах
latitude(j)=latitudetek*180/pi; % Широта в градусах
end;

% Построение карты мира
load coast0;
plot(lam,phi,'k');
axis([-180,180,-90,90]);
hold on;

% Построение наземной проекции траектории спутника
plot(longitude,latitude);