msattrack.m

航天工程工具箱

         hl=plot3(grd.x,grd.y,grd.z);
         set(hl,'color',col.grid)
      end
   end
   
   t=0;
   I=0;
   satpt=0;
   tbeg=clock;
   lon0=fixlon(lontime(zone,isdls))*d2r;
   butt='normal';
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %       START OF ITERATION        %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   while dur>(now-datenum2(tbeg))*24*3600 & strcmp(butt,'normal')
      butt=get(hfig,'selectiontype');
      lon0p=lon0;
      %lon0=fixlon(lontime(zone,isdls,epoch,tcomp))*d2r;
      [sat,satpt,dt,smin,lon0]=synchronize(oe,mu,axunit,N,satpt,tbeg,epoch,zone,isdls,tcomp);  %I update lon0 via this function for improving speed
      if dt          %update, only if position is changed
         %SHOW SATELLITE POSITION
         I=I+1;       %increase "anomaly" index
         for ii=1:N.sat
            [TH(ii) PHI(ii) R(ii)]=cart2sph(sat(ii).r(1),sat(ii).r(2),sat(ii).r(3));    %spherical coordinates for satellit
            pos(ii).alt(I,:)=R(ii)-Rref(PHI(ii),Re,Rp);                            %satellite altitude
            pos(ii).lon(I,:)=round(deg2dms(fixlon((TH(ii)-lon0)*r2d)));    %satellite longitude
            pos(ii).lat(I,:)=round(deg2dms(phi(PHI(ii),Re,Rp)*r2d));             %satellite latitude
            pos(ii).t(I,:)=compress(utc(zone,isdls,'YMDhms'),epoch(smin,:),tcomp);
            %pos(ii).t(I,:)=utc(zone,isdls,'YMDhms');
            pos(ii).r(I,:)=sat(ii).r;
            pos(ii).v(I,:)=sat(ii).v;
         end
         
         
         if I>1                    %don't delete before first point is plotted
            if issat, delete(sh),end
            if israd, delete(rh),end
            if islbl, delete(th),end
            if istrk, delete(trh),end
            if isgwm, delete(gwh),end
            if ismap, delete(mh),end
            if isgrd, delete(grh),end
            if issun, delete(sunh),end
            if isvis, delete(vish),end
            if isssp, delete(ssph),end
         end
         
         %THE SATELLITE ITSELF
         if issat
            for ii=1:N.sat
               sh(ii)=surf(xs+sat(ii).nr(1),ys+sat(ii).nr(2),zs+sat(ii).nr(3));   %plot satellite
               set(sh(ii),'facec',col.satellite,'edgec','none','facel','flat')
            end
         end
         %SATELLITE TRAIL
         if istrail
            for ii=1:N.sat
               tlh=plot3(sat(ii).nr(1),sat(ii).nr(2),sat(ii).nr(3));         %plot satellite trail
               set(tlh,'color',col.trail)
            end
         end
         %RADIUS VECTOR
         if israd
            for ii=1:N.sat
               rh(ii)=plot3([0 sat(ii).nr(1)],[0 sat(ii).nr(2)],[0 sat(ii).nr(3)]);  %plot radius vector
               set(rh(ii),'color',col.radius)
            end
         end
         %SATELLITE LABEL
         if islbl
            for ii=1:N.sat
               th(ii)=text(sat(ii).nr(1),sat(ii).nr(2),sat(ii).nr(3),['  ' name{ii}]);           %plot satellite name
               set(th(ii),'color',col.labels)
            end
         end
         %SATELLITE POSITION TABLE
         if istbl
            disp([datestr(datenum2(pos(smin).t(I,:)),'dd-mmm-yyyy HH:MM:SS') '  :'])
            for ii=1:N.sat
               tname=fixstrlen(name{ii},20);
               talt=fixstrlen(num2str(pos(ii).alt(I)/1e3,6),10);
               tvel=fixstrlen(num2str(norm(pos(ii).v(I,:))/1e3,6),10);
               disp([' ' tname(1:20) ' ' talt ' ' tvel ' ' ...
                     fixstrlen([num2str(pos(ii).lon(I,1)) ',' ...
                                num2str(pos(ii).lon(I,2)) ',' ...
                                num2str(pos(ii).lon(I,3))],12) ...
                     fixstrlen([num2str(pos(ii).lat(I,1)) ',' ...
                                num2str(pos(ii).lat(I,2)) ',' ...
                                num2str(pos(ii).lat(I,3))],10)])
            end
         end
         %SUN POSITION
         svec=normalize(lightpos(lh,zone,isdls))*1.2*usc;
         if issun       %issync is a necessary condition for issun.
            sunh=plot3([0 svec(1)],[0 svec(2)],[0 svec(3)],'y');
            set(sunh,'linew',5)
         end
         %MAP CONTOUR
         if ismap
            [map.x map.y map.z]=sph2cart(mlon+lon0,phi(mlat,Re,Rp),Rref(mlat,Re,Rp)/axunit);      %evaluate cartesian coordinates for coastline
            mh=plot3(map.x,map.y,map.z);       %plot the coastline
            set(mh,'color',col.map)
         end
         %GREENWICH MERIDIAN
         gw(:,3)=phi(lat0,Re,Rp);                  %Greenwich meridian latitudes
         gw(:,2)=lon0;                             %GW meridian longitude
         gw(:,1)=Rref(lat0,Re,Rp)/axunit;              %GW meridian radius (norm)
         gwm.phi(I,:)=round(deg2dms(lon0*r2d));    %store info for output
         gwm.t(I,:)=pos(1).t(I,:);                 % "     "    "    "
         [gw(:,1) gw(:,2) gw(:,3)]=sph2cart(gw(:,2),gw(:,3),gw(:,1));    %Cartesian coords for GW meridian
         if isgwm
            gwh=plot3(gw(:,1),gw(:,2),gw(:,3));                      %plot GW meridian
            set(gwh,'color',col.gwm)
         end
         %LON-LAT GRID
         if isgrd    %grid mode
            grdx=(-180:30:180)*d2r+lon0;
            if isgwm
               grdx=grdx(~~(grdx-lon0));             %in order to leave out the GW meridian.
            end
            grdy=linspace(-pi/2,pi/2,100);
            o=ones(size(grdy));
            for ii=1:length(grdx)
               [grd.x grd.y grd.z]=sph2cart(grdx(ii)*o,phi(grdy,Re,Rp),Rref(grdy,Re,Rp)/axunit);
               grh(ii)=plot3(grd.x,grd.y,grd.z);
               set(grh(ii),'color',col.grid)
            end
         end
         %GROUND TRACK
         for ii=1:N.sat
            [tr(ii).c(I,1) tr(ii).c(I,2) tr(ii).c(I,3)]= ...
               sph2cart(TH(ii),phi(PHI(ii),Re,Rp),Rref(PHI(ii),Re,Rp)/axunit);    %evaluate ssp in Cart coords
            if istrk
               trh(ii)=plot3(tr(ii).c(:,1),tr(ii).c(:,2),tr(ii).c(:,3));                %plot satellite track
               set(trh(ii),'color',col.track)
            end
            [tr(ii).s(1:I,1) tr(ii).s(1:I,2) tr(ii).s(1:I,3)]= ... 
               cart2sph(tr(ii).c(:,1),tr(ii).c(:,2),tr(ii).c(:,3));       %convert all ssp:s
            %[tr(ii).c(1:I,1) tr(ii).c(1:I,2) tr(ii).c(1:I,3)]= ...
               %sph2cart(tr(ii).s(:,1)+we*dt,tr(ii).s(:,2),tr(ii).s(:,3)); %update ssp longitude
            [tr(ii).c(1:I,1) tr(ii).c(1:I,2) tr(ii).c(1:I,3)]= ...
               sph2cart(tr(ii).s(:,1)+lon0-lon0p,tr(ii).s(:,2),tr(ii).s(:,3));
         end
         %SHOW SSP
         if isssp
            for ii=1:N.sat
               [ssp(ii,1) ssp(ii,2) ssp(ii,3)]= ...
                  sph2cart(TH(ii),phi(PHI(ii),Re,Rp),Rref(PHI(ii),Re,Rp)/axunit);
               ssph(ii)=plot3(ssp(ii,1),ssp(ii,2),ssp(ii,3),'.');
               set(ssph(ii),'color',col.ssp)
            end
         end
         %VISIBILITY REGIONS
         if isvis
            for ii=1:N.sat
               vis=visibility(pos(ii),lon0,I,elev,Ry,Rz,axunit);
               vish(ii)=plot3(vis(:,1),vis(:,2),vis(:,3));
               set(vish(ii),'color',col.vis)
            end
         end
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
         %GENERATE ALTITUDE RECORD
         for ii=1:N.sat
            alt(ii).x(I)=pos(ii).alt(I);
            alt(ii).t(I)=(datenum2(pos(ii).t(I,:))-datenum2(pos(ii).t(1,:)))*24*3600;
         end
         
         %axis auto
         axis equal
         axis(ax)         
      end   %end of dt~=0
      
      if getver>=6.5         %only observed this bug in matlab 6.5
         rotate3d on
      end
      drawnow
      
      if isdrift
         [sat,oe]=perturb(sat,oe,N,dt,J2,Re,dWdt,dwdt,dMdt);
      end
   end   %end of while loop
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%   COLLECT DATA FOR OUTPUT   %%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if nargout
      %[cfg1,pos1,P1,gwm1]
      cfg1=cfg;   % see above
      pos1=pos;   % see above
      for ii=1:N.sat
         P1(ii)=sat(ii).P;
      end
      gwm1=gwm;   %see above
   end
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%   PREPARE 2D MAPS   %%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if I>1
      % PLOT LONG-LAT MAP...
      clear t x
      for ii=1:N.sat
         lon=dms2deg(pos(ii).lon);
         lat=dms2deg(pos(ii).lat);
         lonc{ii}=[]; latc{ii}=[];
         j=[0 find(abs(diff(lon))>90)' length(lon)];
         for J=1:length(j)-1
            lonc{ii}=[lonc{ii}; NaN; lon(j(J)+1:j(J+1))];
            latc{ii}=[latc{ii}; NaN; lat(j(J)+1:j(J+1))];
         end
         clear lon lat
         len(ii)=length(lonc{ii});
         t(:,ii)=alt(ii).t(:);
         x(:,ii)=alt(ii).x(:);
      end
      for ii=1:N.sat
         lon(:,ii)=[lonc{ii}; repmat(nan,max(len)-length(lonc{ii}),1)];
         lat(:,ii)=[latc{ii}; repmat(nan,max(len)-length(lonc{ii}),1)];
      end
      clear lonc latc
      
      %plot the 2D track curves.
      figure, hold on
      plot(lon,lat)
      plot3(mlon*r2d,mlat*r2d,-ones(size(mlat)),'g')
      for ii=1:N.sat
        j=find(~isnan(lon(:,ii)));
        i=find(diff(j)==1);         %to assure that you can see which ring corresponds to a certain track
        i=j(i(1));
        plot(lon(i,ii),lat(i,ii),'ko')
      end
      legend(char(name),0)
      xlabel('\lambda [deg]'), ylabel('\phi [deg]')
      title(['Ground Track Map'])
      axis equal
      axis tight
      axis([-180 180 -90 90])
      set(gca,'xtick',-180:60:180)
      set(gca,'ytick',-90:30:90)
      grid
      hold off
      
      %ALTITUDE(TIME) RECORD
      figure
      plot(t,x/1e3)
      legend(char(name),0)
      xlabel('t [s]'), ylabel('h [km]')
      title(['Altitude Record'])
      grid
      hold off
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%    THE END    %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
popup(ver,coastl,mu)
clear global mquit h


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   FUNCTIONS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[y]                   = NORMALIZE   (x)
%[oe,name,epoch]       = SPLITSTRUCT (sat)
%[sat,oe]              = PERTURB     (sat,oe,N,dt,J2,Re,dWdt,dwdt,dMdt)
%[sat,i,dt,iii,lon0]   = SYNCHRONIZE (oe,mu,axunit,N,i0,YMDhms_now,epoch,zone,isdls,tcomp)
%[d]                   = DATENUM2    (YMDhms)
%[v]                   = VISIBILITY  (pos,lon0,I,el,Ry,Rz,axunit)
%[y]                   = FIXLON      (x)
%[pos]                 = LIGHTPOS    (l,zone,isdls)
%[lon0]                = LONTIME     (zone,isdls,epoch0,tcomp)
%[YMDhms2]             = COMPRESS    (YMDhms,epoch0,tcomp)
%[sec]                 = DAY2SEC     (dayn)
%[dayn]                = SEC2DAY     (sec)
%[ax]                  = MAXIS       (oe)
%[col]                 = SETTINGS
%[N,ismrk,isdyn,isorb] = MENU1
%[tcomp,dur,zone,isdls,
% issun,ismap,isgrd,
% isgwm,istrk,issat,
% istbl,israd,isdrift,
% istrail,isvis,elev,
% isssp,mquit]         = MENU2
%                        POPUP       (ver,coastl,mu)

function y=normalize(x)
y=x/norm(x);


function [oe,name,epoch]=splitstruct(sat)
%break down structure into its contents and gather the data in matrices
for i=1:length(sat);
   oe(i,:)=sat(i).oe;
   name{i}=sat(i).name;
   epoch(i,:)=sat(i).epoch;
end


function [sat,oe]=perturb(sat,oe,N,dt,J2,Re,dWdt,dwdt,dMdt)
d2r=pi/180;
for ii=1:N.sat
   n=sat(ii).n; t=sat(ii).t;
   if isempty(dt), dt=t;end
   a=oe(ii,1); e=oe(ii,2); i=oe(ii,3)*d2r;
   oe(ii,4)=oe(ii,4)+dWdt(a,e,i,n,J2,Re)*dt;
   oe(ii,5)=oe(ii,5)+dwdt(a,e,i,n,J2,Re)*dt;
   %oe(ii,6)=oe(ii,6)+dMdt(a,e,i,n,J2,Re)*dt;
end


function [sat,i,dt,iii,lon0]=synchronize(oe,mu,axunit,N,i0,YMDhms_now,epoch,zone,isdls,tcomp)
Pe=24*3600;
d2r=pi/180;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Get period times for satellites
[r,v,P]=orb2cart(oe(1,:),0,mu);
Pmin=P;
for ii=1:N.sat
   [r,v,P]=orb2cart(oe(ii,:),0,mu);
   sat(ii).P=P;
   sat(ii).n=360/sat(ii).P;
   Pmin=min(Pmin,P);
   if Pmin==sat(ii).P
      iii=ii;
   end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[r,v,P,dt]=orb2cart(oe(iii,:),N.pts,1,mu);
lon0=fixlon(lontime(zone,isdls,YMDhms_now,tcomp))*d2r;
t_n=day2sec(datenum2(compress(utc(zone,isdls,'YMDhms'),YMDhms_now,tcomp)));
%t_n=utc(zone,isdls,'s');
for ii=1:N.sat
   %tested and seems to work
   t_e=day2sec(datenum2(epoch(ii,:)));
   t=t_n-t_e;
   tt=mod(t,sat(ii).P);
   [r,v]=orb2cart(oe(ii,:),oe(ii,6)+sat(ii).n*tt,mu);
   if ii==iii                %for the fastest orbit...
      i=mod(round(tt/Pmin*(N.pts-1))+round(oe(iii,6)/360*N.pts),N.pts-1)+1;         %tested and works
      %[r,v]=orb2cart(oe(iii,:),N.pts,i,mu);           %calculate current position and velocity
      if i0==i, dt=0;end                              %let system stand still if index remain constant
   %else
      %[r,v]=orb2cart(oe(ii,:),sat(ii).n*tt,mu);
   end
   sat(ii).r=r;
   sat(ii).nr=r/axunit;
   sat(ii).v=v;
   sat(ii).t=t;