groundtrack.m

衛星軌道運算可算出課譜樂六元素的一個好程式

% Orbit groundtrack plot Latitude longitude lat long
% Richard Rieber
% December 18, 2006
% [e-mail]
%
% Revision 8/21/07: Fixed small typo in help comments
%                   Added H1 line for lookfor functionality
%
% function groundtrack(Kepler,GMSTo,Tf,fig)
%
% Purpose: This function plots the groundtrack of a given orbit.
% 
% Inputs: o Kepler - A vector of length 6 containing all of the keplerian
%                    orbital elements [a,ecc,inc,Omega,w,nu] in km and radians
%         o GMSTo - The Greenwich Mean Siderial Time at the given initial position
%                   in radians.
%         o Tf - The length of time to plot the groundtrack in seconds
%         o fig - The figure number on which to plot the groundtrack [OPTIONAL]
%
% Outputs: NONE

function groundtrack(Kepler,GMSTo,Tf,fig)

if nargin < 3 || nargin > 4
    error('Incorrect number of inputs, see help groundtrack.m')
elseif nargin == 3
    figure
elseif nargin == 4
    figure(fig)
end

a = Kepler(1);
ecc = Kepler(2);
inc = Kepler(3);
O = Kepler(4);
w = Kepler(5);
nuo = Kepler(6);

r2d = 180/pi;

R_earth = 6378.1363;       %km  Radius of Earth
mu_earth = 398600.4415;    %km^3/s^2  Gravitational constant of Earth
w_earth = 7.2921158553e-5; %rad/s  %Rotation rate of Earth

n = (mu_earth/a^3)^.5;  %Mean motion of satellite

E = atan2(sin(nuo)*(1-ecc^2)^.5,ecc+cos(nuo));  %Eccentric anomaly
MA = E-ecc*sin(E); %Initial Mean anomaly

dT = 60; %time step in sec

time = [0:dT:Tf]; %time vector

bp(1) = 0;
bp(2) = length(time);
k = 2;

for j = 1:length(time)
    GMST = zeroTo360(GMSTo + w_earth*time(j),1);  %GMST in radians
    
    M = zeroTo360(MA + n*time(j),1);  %Mean anomaly in rad
    nu = nuFromM(M,ecc); %True anomaly in rad
    
    [ECI,Veci] = randv(a,ecc,inc,O,w,nu);
    ECEF = eci2ecef(ECI,GMST);
    
    [Lat_rad(j),Long_rad(j)] = RVtoLatLong(ECEF,GMST);
    
    if j > 1 && (Long_rad(j)-Long_rad(j-1)) < -pi
        bp(k) = j-1;
        k = k+1;
    end

end
bp(length(bp)+1) = length(time);
Lat = Lat_rad.*r2d;
Long = Long_rad.*r2d;

x = load('Coastline.dat');
plot(x(:,1),x(:,2),'g')
% polar(x(:,1)*pi/180,x(:,2),'g')
clear x

hold on
for j = 2:length(bp)
%     polar(Long(bp(j-1)+1:bp(j))*pi/180,Lat(bp(j-1)+1:bp(j)),'b')
	plot(Long(bp(j-1)+1:bp(j)),Lat(bp(j-1)+1:bp(j)),'b')
end

axis([-180,180,-90,90])

xlabel('Longitude (\circ)')
ylabel('Latitude (\circ)')