msattrack.m

航天工程工具箱

function [cfg1,pos1,P1,gwm1]=msattrack(sat,name,epoch)
%MSATTRACK  Satellite tracker.
%   MSATTRACK(SAT), where SAT is an array of structs:
%
%      SAT(1:N).oe
%      SAT(1:N).name
%      SAT(1:N).epoch
%
%   as will be explained below. N is the number of satellites.
%   Alternatively:
%
%   MSATTRACK(OE[,NAME[,EPOCH]]), where OE is the orbital elements (Nx6)
%   for each satellite with name NAME {N} (strcell), and EPOCH (Nx6)
%   is in YMDhms format, eg:
%
%      [2003 4 13 18 6 55.8708]
%
%   The orbital elements have the format:
%
%      a [m]   : semi-major axis
%      e []    : eccentricity
%      i [deg] : inclination
%      W [deg] : longitude of the ascending node
%      w [deg] : argument of perigee
%      M [deg] : mean anomaly at EPOCH
%
%   If in dynamic mode you can get following outputs:
%
%   [CFG,POS,P,GWM]=MSATTRACK(...)
%   where:
%
%      CFG : configuration settings
%      POS : positions as shown during simulation
%            altitude        [m]
%            longitude       [dms]
%            latitude        [dms]
%            time point      [YMDhms] (UTC)
%            radius vector   [m]
%            velocity vector [m/s]
%      P   : period-time(s) for satellite(s)
%      GWM : Greenwich meridian pos. relative to vernal equinox
%
%   To end the simulation at any time simply left click with the mouse
%   on the simulation figure.
%
%   See also TLE2ORB, ORB2CART, CART2ORB, SATVIS, SATSORT, MA2TP, TP2MA.

% Copyright (c) 2002-12-07, Rasmus Anthin.
% Revision 2002-12-08 - 2002-12-10, 2002-12-18,
%          2002-12-19 - 2002-12-23, 2002-12-30, 2003-03-25,
%          2003-04-06 - 2003-04-10, 2003-04-13,
%          2003-04-15 - 2003-04-19,
%          2003-06-03 - 2003-06-09,
%          2003-06-18 - 2003-06-19, 2003-06-21, 2003-07-21.





% Model and Program Limitations:
%  3. Earth/geoid is ellipsoidal/oblative: a=b=Re, c=Rp.
%  5. no luni-solar perturbations.
%  6. no orbital decay due to drag.
%  9. no east-west drift.

% Bugs:
%  -When using perturbations, satellite will go too fast
%   (=go berserk)
%   if including mean anomaly drift (dM/dt).
%  -SSP longitude for ISS is shifted about 2 deg compared
%   to the NASA satellite tracker. (earth mass??)
%  -Time compression doesn't have a well defined time origin.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   CONSTANTS   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r2d=180/pi;
d2r=1/r2d;
z=[0 0];
Re=constant('Ree');       %equatorial earth radius  (used for normalizing)
Rp=constant('Rep');       %polar earth radius
J2=1082.626523e-6;        %Earth zonal harmonic coefficient #2.
J3=-.254e-5;
J4=-.161e-5;
mu=constant('Me')*constant('G');      %specific gravity constant of Earth
Pse=convert(1,'yr_sid','s');          %Earth orbital period     (sideral)
Pe=hms2sec([23 56 4.09054]);          %Earth rotational period  (sideral)
we=2*pi/Pe;                           %Earth rotational frequency (sideral)
lat0=linspace(-pi/2,pi/2)';    %GWM latitudes
ver='5.2';
load map


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   INLINE FUNCTIONS   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rref=inline('sqrt((E*cos(f)).^2+(P*sin(f)).^2)','f','E','P');          %Earth radius at ssp
phi=inline('atan(Re/Rp*tan(eta))','eta','Re','Rp');                    %ellip lat as func of spheric lat
dWdt=inline('-3*n*J2*Re^2*cos(i)/(2*a^2*(1-e^2)^2)','a','e','i','n','J2','Re');      %regression of line of nodes
dwdt=inline('3*n*J2*Re^2*(4-5*sin(i)^2)/(4*a^2*(1-e^2)^2)','a','e','i','n','J2','Re');   %perigee drift
dMdt=inline('n*(1+3/2*J2*(Re/a)^2*[1-3/2*sin(i)^2]/sqrt(1-e^2)^3)','a','e','i','n','J2','Re');   %drift of mean anomaly
Rx=inline('[1 0 0;0 cos(x) -sin(x);0 sin(x) cos(x)]');
Ry=inline('[sin(x) 0 cos(x);0 1 0;cos(x) 0 -sin(x)]');
Rz=inline('[cos(x) -sin(x) 0;sin(x) cos(x) 0;0 0 1]');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   PROCESS INPUTS   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% oe, name, epoch, issync, N.sat
if ~nargin | isempty(sat)
   sat=[(mu*(Pe/2/pi)^2)^(1/3) zeros(1,5)];         %demo: GEO-orbit
   name={'GEO'};
end
if isstruct(sat)
   [oe,name,epoch]=splitstruct(sat);                %extract structure element
else
   oe=sat;
end
if nargin>1 & ~iscell(name) & (ischar(name) | isempty(name))
   name={name};
elseif nargin<2 & ~exist('name')
   for i=1:size(oe,1)
      name{i}=['sat' num2str(i-1)];
   end
end
N.sat=size(oe,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   GET BASIC SETTINGS   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% N.pts, islbl, isdyn, isorb
[N.pts,islbl,isdyn,isorb,axunit,mquit]=menu1;                   %basic preferences
if mquit, return,end
units=[Re 1e6 convert(1e3,'mile','m')];
ustr={'R_\oplus','Mm','kmile'};
ustr=ustr{axunit};
axunit=units(axunit);
usc=Re/axunit;
cfg.basic.info=char({'BASIC CONFIG', ...
                     '------------', ...
                     'N      : number of points along orbit in one revolution for each satellite [#]', ...
                     'islbl  : show labels on each satellite [0-1]', ...
                     'isdyn  : dynamic mode [0-1]', ...
                     'isorb  : show nominal (unperturbed) orbital planes for each satellite [0-1]', ...
                     'axunit : unit for the axis in the graph, 1 [R_e, Mm, kmiles]'});
cfg.basic.N=N.pts;
cfg.basic.islbl=islbl;
cfg.basic.isdyn=isdyn;
cfg.basic.isorb=isorb;
cfg.basic.axunit=axunit;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   DISPLAY PERIOD TIMES   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:N.sat
   [r,v,P]=orb2cart(oe(ii,:),0,mu);
   Phms=sec2hms(P);
   Pstr=['P = [' ...
         num2str(Phms(1)) ' h, ' ...
         num2str(Phms(2)) ' min, ' ...
         num2str(Phms(3)) ' s]' ...
         '  : ' name{ii}];
   disp([blanks(3) Pstr ])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   GET DYNAMICS SETTINGS   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tcomp, dur, zone, isdls, issun, ismap, isgrd, istrk,
% issat, istbl, israd, isdrift, istrail, isvis, elev,
% isssp, mquit
if isdyn                                           %dynamics preferences
   [tcomp,dur,zone,isdls,issun,ismap,isgrd,isgwm,istrk,issat, ...
         istbl,israd,isdrift,istrail,isvis,elev,isssp,mquit]=menu2;
   if ~isvis, elev=NaN;end
else
   zone=0;
   isdls=0;
end
if mquit, return, end
cfg.app.info=char({'APPEARENCE CONFIG', ...
                   '-----------------', ...
                   'tcomp   : time compression factor [#]', ...
                   'dur     : duration of simulation [hms]', ...
                   'zone    : time zone [h]', ...
                   'isdls   : daylight saving [0-1]', ...
                   'issun   : show sun direction vector [0-1]', ...
                   'ismap   : show coastline contour map [0-1]', ...
                   'isgrd   : show longitude-latitude grid [0-1]', ...
                   'isgwm   : mark out the greenwhich meridian by using another color [0-1]', ...
                   'istrk   : show ground track [0-1]', ...
                   'issat   : show satellite (as small spheres) [0-1]', ...
                   'istbl   : show current state vector as table during runtime [0-1]', ...
                   'israd   : show the radius vector for each satellite [0-1]', ...
                   'isdrift : apply orbital perturbations due to the earth''s bulge [0-1]', ...
                   'istrail : show trail of points behind each satellite [0-1]', ...
                   'isvis   : show visibility regions [0-1]', ...
                   'elev    : minimum elevation for visibility (see: isvis) [deg]', ...
                   'isssp   : show subsatellite point [0-1]'});
if isdyn
   cfg.app.tcomp=tcomp;
   cfg.app.dur=sec2hms(dur);
   cfg.app.zone=zone;
   cfg.app.isdls=isdls;
   cfg.app.issun=issun;
   cfg.app.ismap=ismap;
   cfg.app.isgrd=isgrd;
   cfg.app.isgwm=isgwm;
   cfg.app.istrk=istrk;
   cfg.app.issat=issat;
   cfg.app.istbl=istbl;
   cfg.app.israd=israd;
   cfg.app.isdrift=isdrift;
   cfg.app.istrail=istrail;
   cfg.app.isvis=isvis;
   cfg.app.elev=elev;
   cfg.app.isssp=isssp;
else
   cfg.app.tcomp=NaN;
   cfg.app.dur=NaN;
   cfg.app.zone=zone;
   cfg.app.isdls=isdls;
   cfg.app.issun=NaN;
   cfg.app.ismap=NaN;
   cfg.app.isgrd=NaN;
   cfg.app.isgwm=NaN;
   cfg.app.istrk=NaN;
   cfg.app.issat=NaN;
   cfg.app.istbl=NaN;
   cfg.app.israd=NaN;
   cfg.app.isdrift=NaN;
   cfg.app.istrail=NaN;
   cfg.app.isvis=NaN;
   cfg.app.elev=NaN;
   cfg.app.isssp=NaN;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   GET PARAMETER SETTINGS   %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%----------- GET COLORS ----------
col=settings;                      %load color settings

%----------- EPOCH ---------------
if ~exist('epoch')
   epoch=utc(zone,isdls,'YMDhms');
   epoch=repmat(epoch,size(oe,1),1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   PREPARE 3D ENVIRONMENT   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hfig=figure;
hold on
set(gcf,'name',['MSatTrack ver. ' ver],'num','off')     %set title
[xe,ye,ze]=sphere(100);xe=xe*Re;ye=ye*Re;ze=ze*Rp;   %create earth structure
earth=surf(xe/axunit,ye/axunit,ze/axunit);                       %plot earth
shading interp                                       %smoothen the surface
set(earth,'facecolor',col.earth,'edgecolor','none')  %set surface color
lighting phong                                       %create light reflection
lh=light('pos',[1 0 0]);    %light from sun in spring time
set(gca,'color',col.bg)     %set background color
he=circle(usc,z);             %set equator color
set(he,'color',col.equator)
orbpts=100;                 %used to be =N.pts
if isorb
   for i=1:N.sat
      r=orb2cart(oe(i,:),orbpts,[],mu);
      r=r([1:orbpts 1],:);
      ho=plot3(r(:,1)/axunit,r(:,2)/axunit,r(:,3)/axunit);    %satellite nominal orbits
      set(ho,'color',col.orbit)
   end
end
if islbl
   for i=1:N.sat
      p=orb2cart(oe(i,:),180,mu);
      th(i)=text(p(1)/axunit,p(2)/axunit,p(3)/axunit,['' name{i}]);
      set(th(i),'color',col.labels)
   end
end
dim=[' [' ustr ']'];
xlabel(['x' dim]),ylabel(['y' dim]),zlabel(['z' dim])
plot3(z,z,[-1.3 1.3]*usc,'k')              %line of poles
plotvec([-1.8 0]*usc,[1.8 0]*usc,'k')      %line of vernal equinox
text(2*usc,0,0,'\gamma')                   %vernal equinox label
text(0,0,1.5*usc,'N','fontw','b')          %north pole label
text(0,0,-1.5*usc,'S','fontw','b')         %south pole label
view(3)
ax=maxis(oe,axunit);
axis equal
axis(ax)
rotate3d on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   SATELLITE DYNAMICS   %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isdyn
   if islbl, delete(th);end
   [xs ys zs]=sphere(10);
   sc=.05*usc;
   %sc=.01*max(diff(reshape(ax,2,3)));           %compensate for axis dependency
   xs=xs*sc;ys=ys*sc;zs=zs*sc;
   
   sat=synchronize(oe,mu,axunit,N,1,clock,epoch,zone,isdls,1);
   % sat.r   current position
   % sat.v   current velocity
   % sat.P   period time
   % sat.t   time elapsed since epoch
   % sat.n   mean angular frequency
   
   %PREPARE TABLE
   if istbl
      disp(' ')
      disp(' SATELLITE            ALT[km]    VEL[km/s]  LON[dms]    LAT[dms]')
      disp(repmat('=',1,64))
   end
  
   mlon=coastl(:,1)*d2r;
   mlat=coastl(:,2)*d2r;
   
   tstr='Satellite dynamics ';
   if isdrift
      tstr=[tstr 'with '];
   else
      tstr=[tstr 'w/o '];
   end
   tstr=[tstr 'perturbations'];
   figure(hfig)
   title(tstr)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %  CONSIDER EARLIER PERTURBATIONS  %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if isdrift   %synchronized mode:  account for earlier regressions/perturbations
      sat=synchronize(oe,mu,axunit,N,0,clock,epoch,zone,isdls,1);
      [sat,oe]=perturb(sat,oe,N,[],J2,Re,dWdt,dwdt,dMdt);
   end

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %  GENERATE LONGITUDE LINES %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if isgrd    %prepare grid. only plot longitude lines
      grdy=(-90:30:90)*d2r;
      grdy=grdy(~~grdy);      %in order to leave out the equatorial line
      grdx=linspace(-pi,pi,100);
      o=ones(size(grdx));
      for i=1:length(grdy)
         [grd.x grd.y grd.z]=sph2cart(grdx,phi(grdy(i),Re,Rp)*o,Rref(grdy(i),Re,Rp)/axunit*o);