www.pudn.com > final-GPS.rar > SV_Ephemeris_Model.m, change:2008-04-04,size:9679b

```%This Function use Ephemeris Data and Calculate satellite Position
%CopyRight By Moein Mehrtash
%**************************************************************************
% Written by Moein Mehrtash, Concordia University, 3/28/2008              *
% Email: moeinmehrtash@yahoo.com                                          *
%**************************************************************************
%**************************************************************************
% Satellite Position By Ephemeris Model
%Function's Inputs:
% GPS_RVC_T:(Sec):
% A:semi-major orbit axis (meter)
% n0: Compute mean motion (rad/s)
% tk:time from ephemeris reference epoch
% n: corrected mean motion
% Dealta_n:correction for mean motion
% M0:mean Anomaly at refernce time
% Mk:Mean anomaly
% Ek: Keplers equation for the eccentric anomaly Ek (rad)
% e: orbit eccentricity
% Toe: Refernce time of ephemeris parameters
% i0: inlination angle at reference time
% Phik:Argument of lattitude
% Omega0: Longitude of Ascending Node
% omega: Argument of perigee (semicicrcle)
% Omegadot: rate of right ascension
% IDOT: rate of inlination angle
% Delta_uk:Argument of lattitude correction
% Delta_rk:Argument of radius correction
% Delta_ik:Argument of inclination correction
% uk:Corrcted argument of lattitude
% rk: corrected radius
% ik:corrected inclination
% xpk:satellite position in orbital plane
% ypk:satellite position in orbital plane
% SEcond Harmonic Perturbation coefiicient
%C_us:Amplitude of the cosine harmonic correction term to the argument
%     of lattitude(Rad)
%C_uc:Amplitude of the sine harmonic correction term to the argument
%     of lattitude(Rad)
%C_rs:Amplitude of the cosine harmonic correction term to orbit
%radius(meter)
%C_rc:Amplitude of the sine harmonic correction term to orbit
%radius(meter)
%C_is:Amplitude of the cosine harmonic correction term to the angle of
%     inclination(Rad)
%C_ic:Amplitude of the sine harmonic correction term to the angle of
%     inclination(Rad)

%Function's Outputs:
% xk,yk,zk: Satellite position in ECEF
% Ek:Eccentric Anomaly
% A:semi-major orbit axis
% e: orbit eccentricity
% ik: Corrected Inlination
% Omegak: Corrected longitude of ascending node
% Nuk:True anomaly as a function of the eccentric anomaly
%Function's Constants:
% Mu:WGS-84 value of the earths universal gravitational
%    parameter(3.986005*10^14 m4/s2)
% Omega_dote:WGS-84 value of the earths roration rate
%           (7.2921151467*10^(-5)rad/s)
%Input:(GPS_RVC_T,[C_rs Dealta_n M0 C_uc e C_us sqrt(A) Toe C_ic Omega0 C_is
%                  i0 C_rc omega Omegadot IDOT]

%**************************************************************************
function [Pos_xyz_Mat,Orbit_parameter]=SV_Ephemeris_Model(GPS_RVC_T,data);

[Parameter,ST_NO]=size(data);
% Constant Parameter of this function
Mu=3.986005*10^14;
Omega_dote=7.2921151467*10^(-5);

for ST_No=1:ST_NO

%Load Data
A=data(7,ST_No)*data(7,ST_No);
t=GPS_RVC_T;
toe=data(8,ST_No);
Delta_n=data(2,ST_No);
M0=data(3,ST_No);
omega=data(14,ST_No);
C_us=data(6,ST_No);
C_uc=data(4,ST_No);
C_rs=data(1,ST_No);
C_rc=data(13,ST_No);
C_is=data(11,ST_No);
C_ic=data(9,ST_No);
i0=data(12,ST_No);
IDOT=data(16,ST_No);
Omega0=data(10,ST_No);
Omegadot=data(15,ST_No);
e=data(5,ST_No);

%**************************************************************************
%**************************************************************************
%Compute Time From Ephemeris Refernce Epoch
%Input:
%t(sec):GPS Sys. Time
%toe(sec):Reference Time Ephemeris
%Output:
%tk(sec):Time from ephemeris reference epoch
tk=t-toe;
%**************************************************************************
%**************************************************************************
%Compute corrected mean motion
%Input:
%A(meter): Semi-major axis
%Mu(meter^3/sec^2): value of Earth's universal gravitational parameters
%tk(sec):Time from Ephemeris Reference Epoch
%Delta_n(semi-circle/sec):Mean Motion Diference From Computed Value
%M0(semi-circles):
%Output:
%n0(Rad/sec):Computed mean motion
%n(semi-circle/sec):Corrected mean motion
%Mk(semi-circle):Mean anomaly
n0=sqrt(Mu/(A^3));
n=n0/pi+Delta_n;
Mk=M0+n*tk;

%**************************************************************************
%**************************************************************************
%Solve Kepler Eq. with Initial value Ek0=Mk
%Use Fzero Function to solve Kepler Eq.
%The fzero command is an M-file. The algorithm, which was originated by
%T. Dekker, uses a combination of bisection, secant, and inverse quadratic
%interpolation methods. An Algol 60 version, with some improvements, is
%given in [Brent, R., Algorithms for Minimization Without Derivatives,
%Prentice-Hall, 1973.]
%Input:
%Mk(semicircle): Mean Anomaly
%e:  Eccentricity
%Output:
%Ek(Rad):  Eccentric Anomaly
Ek=fzero(@(x) Kepler_Eq(x,Mk*pi,e),Mk*pi);
%**************************************************************************
%**************************************************************************
%Compute True Anomaly as a Function of Eccentric Anomaly
%Input:
%Ek(Rad):  Eccentric Anomaly
%e      :  Eccentricity
%Output:
%nuk(Rad): True Anomaly
Sin_nuk=sqrt(1-e^2)*sin(Ek)/(1-e*cos(Ek));
Cos_nuk=(cos(Ek)-e)/(1-e*cos(Ek));
nuk=atan2(Sin_nuk,Cos_nuk);
if nuk<0
nuk=nuk+2*pi;
end
%**************************************************************************
%**************************************************************************
%Compute Argumet of Lattitude and Correction
%Input:
%nuk(Rad): True Anomaly
%omega(semicircle):  Argument of perigee
%Output:
%Phik(Rad):  Argument of lattitude
Phik=nuk+omega*pi;

%**************************************************************************
%**************************************************************************
%Compute Argumet of Lattitude,Radious and Inclination Correction
%SEcond Harmonic Perturbation
%Input:
%Phik(Rad):  Argument of lattitude
%C_us:Amplitude of the cosine harmonic correction term to the argument
%     of lattitude(Rad)
%C_uc:Amplitude of the sine harmonic correction term to the argument
%     of lattitude(Rad)
%C_rs:Amplitude of the cosine harmonic correction term to orbit
%radius(meter)
%C_rc:Amplitude of the sine harmonic correction term to orbit
%radius(meter)
%C_is:Amplitude of the cosine harmonic correction term to the angle of
%     inclination(Rad)
%C_ic:Amplitude of the sine harmonic correction term to the angle of
%     inclination(Rad)
%Output:
% Delta_uk(Rad):Argument of lattitude correction
% Delta_rk(meter):Argument of radius correction
% Delta_ik(Rad):Argument of inclination correction

Delta_uk=C_us*sin(2*Phik)+C_uc*cos(2*Phik);
Delta_rk=C_rs*sin(2*Phik)+C_rc*cos(2*Phik);
Delta_ik=C_is*sin(2*Phik)+C_ic*cos(2*Phik);
%**************************************************************************
%**************************************************************************
%Compute Corrected Value of Lattitude,Radious,Inclination and Ascendong
%node
%Input:
% Phik(Rad):  Argument of lattitude
% A(meter): Semi-major axis
% Ek(Rad):  Eccentric Anomaly
% e:  Eccentricity
% i0(semi-circle): inlination angle at reference time
% Omega0(semi-circle): Refernce Longitude of Ascending Node
% Omegadot(semi-circle/sec): rate of right ascension
% Omega_dote(rad/sec):WGS 84 value of the earth's rotation rate
% IDOT(semi-circle/sec): rate of inlination angle
% Delta_uk(Rad):Argument of lattitude correction
% Delta_rk(meter):Argument of radius correction
% Delta_ik(Rad):Argument of inclination correction
%Output:
% uk(Rad):Corrcted argument of lattitude
% rk(meter): corrected radius
% ik(Rad):corrected inclination
% Omegak(Rad):Corrected longitude of ascending node.
uk=Phik+Delta_uk;              %Latitude
rk=A*(1-e*cos(Ek))+Delta_rk;   %Radious
ik=i0*pi+Delta_ik+IDOT*tk*pi;  %Inclination
Omegak=Omega0*pi+(Omegadot*pi-Omega_dote)*tk-Omega_dote*toe;

%**************************************************************************
%**************************************************************************
%Compute satellite vehicle position
%Satellite position in orbital plane
x_Perk=rk*cos(uk);
y_Perk=rk*sin(uk);
%Satellite Position in ECEF
xk=x_Perk*cos(Omegak)-y_Perk*cos(ik)*sin(Omegak);
yk=x_Perk*sin(Omegak)+y_Perk*cos(ik)*cos(Omegak);
zk=y_Perk*sin(ik);
%Output:
% xk,yk,zk(meter): Satellite Position in ECEF
% Mk(semicircle): Mean Anomaly
% Ek(Rad):  Eccentric Anomaly
% nuk(Rad): True Anomaly
% Phik(Rad):  Argument of lattitude
% uk(Rad):Corrcted argument of lattitude
% rk(meter): corrected radius
% ik(Rad):corrected inclination
% Omegak(Rad):Corrected longitude of ascending node.
% A(meter):semi-major orbit axis
% e: orbit eccentricity

Pos_xyz=[xk yk zk Mk Ek nuk Phik uk rk ik Omegak A e];

Pos_xyz_Mat(ST_No,1)=Pos_xyz(1);Pos_xyz_Mat(ST_No,2)=Pos_xyz(2);Pos_xyz_Mat(ST_No,3)=Pos_xyz(3);
for j=1:10
Orbit_parameter(j,ST_No)=Pos_xyz(j+3);
end
end

```