msattrack.m

航天工程工具箱

end
dt=dt*tcomp;     %something is seriously wrong with the time comression here.


function d=datenum2(YMDhms)
%this is actually not needed for MATLAB ver 5.3 or later.
d=datenum(YMDhms(1),YMDhms(2),YMDhms(3),YMDhms(4),YMDhms(5),YMDhms(6));


function v=visibility(pos,lon0,I,el,Ry,Rz,axunit)
%Create the visibility ring and return its cartesian coordinates
d2r=pi/180;
lon=dms2deg(pos.lon(I,:))*d2r;
lat=dms2deg(pos.lat(I,:))*d2r;
R=Rz(lon+lon0)*Ry(pi/2-lat);
a=linspace(0,2*pi);
h=pos.alt(I);
r=norm(pos.r(I,:));
Rr=r-h;
th=satvis(el*d2r,Rr,h);
o=ones(size(a));
v0=[Rr*cos(th)*o; Rr*sin(th)*cos(a); Rr*sin(th)*sin(a)];
v=(R*v0)';
%v=v0';
v=v/axunit;
%Rref(lat0,Re,Rp);
%phi(lat0,Re,Rp);


function y=fixlon(x)
%corrects the longitude so that it won't stretch out further than +-180
% eg: 185 -> -175,   -185 -> 175
y=rem(x,360);
if abs(y)>180                 %before: x
   y=sign(y)*(abs(y)-360);    %before: x
end


function pos=lightpos(l,zone,isdls)
%change light position due to earths orbit around the sun
% and tilt with respect to the ecliptic
AU=convert(1,'AU','m');        %distance Earth-sun
alpha=23.45;                   %tilt of Earth
aa=alpha*pi/180;
eps=0.0167;                    %eccentricity of orbit
time=clock;
t.Y=time(1);
t.M=time(2);
t.D=time(3);
t.h=time(4);
t.m=time(5);
t.s=time(6);
M=-3.18+0.98560*(now-datenum(t.Y,0,0));     %...,1,1) -> ...,0,0)
omega=-77.11;
L=M+omega;
[x,y,z]=sph2cart(L*pi/180,0,AU);
pos=[x y*cos(aa)-z*sin(aa) y*sin(aa)+z*cos(aa)];
set(l,'pos',pos)


function lon0=lontime(zone,isdls,epoch0,tcomp)
%initial longitude position relative to vernal equinox
if nargin<3, epoch0=clock;tcomp=1;end     %this routine doesn't work with tcomp>1 just yet.
time=compress(clock,epoch0,tcomp);
d=datenum(time(1),time(2),time(3))-datenum(time(1),0,0);        %...,1,1) -> ...,0,0)
h=hms2sec(time(4)-zone-daylight(zone,isdls),time(5),time(6))/3600;
lon0=99.698+0.98562*d+15.041*h;             %crude but seems to work.


function YMDhms2=compress(YMDhms,epoch0,tcomp)
%YMDhms2=tcomp*(YMDhms-epoch0)+epoch0;
%convert to daynum to seconds
te=datenum2(epoch0);          %daynum
te=day2sec(te);               %sec
tn=datenum2(YMDhms);          %daynum
tn=day2sec(tn);               %sec
t=(tn-te)*tcomp+te;
YMDhms2=datevec(sec2day(t));


function sec=day2sec(dayn)
Pe=24*3600;
sec=floor(dayn)*Pe+rem(dayn,1)*Pe;


function dayn=sec2day(sec)
Pe=24*3600;
d=floor(sec/Pe);
dayn=d+rem(sec/Pe-d,1);


function ax=maxis(oe,axunit)
ax=reshape(axis,2,3);
axmin=ax(1,:);
axmax=ax(2,:);
G=constant('G');
M=constant('Me');
mu=G*M;
for ii=1:size(oe,1)
   r=orb2cart(oe(ii,:),100,[],mu)/axunit;
   rmin=min(r);
   rmax=max(r);
   axmin=min(axmin,rmin);
   axmax=max(axmax,rmax);
end
ax=reshape([axmin-.1;axmax+.1],1,6);


function col=settings
%Get color settings from preset file.
fp=fopen('colors.cfg');
while 1
   line=fgetl(fp);
   if ~ischar(line), break;end
   if length(line) & ~strcmp(line(1),'%')
      name=deblank(line(1:find(line==':')-1));
      value=deblank(line(find(line==':')+1:end));
      value=value(find(~isspace(value)));
      if length(value) & value(1)=='#'
         val=hex2rgb(value(2:7));
      else
         switch lower(value)
         case 'yellow', val='y';
         case 'magenta', val='m';
         case 'cyan', val='c';
         case 'red', val='r';
         case 'green', val='g';
         case 'blue', val='b';
         case 'white', val='w';
         case 'black', val='k';
         case {'y','m','c','r','g','b','w','k'}, val=lower(value);
         otherwise, val='b';
         end
      end
      eval(['col.' name '=val;']);
   end
end


function [N,ismrk,isdyn,isorb,axunit,mquit]=menu1
hfig=dialog('Name','Basic Settings',...
            'HandleVisibility','on');
str1=' Number of points for innermost orbit.';
str2=' Satellite label(s).';
str3=' Dynamic mode.';
str4=' Show nominal orbit plane.';
str5=' Axis units.';
str5a={'R_e','Mm','kmile'};
dy=30;
h(1)=uicontrol('style','edit','min',1,'max',1,'string','100','pos',[20 380 60 20]);
uicontrol('style','text','string',str1,'pos',[110 380 250 20],'horiz','l')
%
h(2)=uicontrol('style','check','pos',[20 380-dy 40 20]);
uicontrol('style','text','string',str2,'pos',[110 380-dy 250 20],'horiz','l')
%
h(3)=uicontrol('style','check','pos',[20 380-2*dy 40 20]);
uicontrol('style','text','string',str3,'pos',[110 380-2*dy 250 20],'horiz','l')
%
h(4)=uicontrol('style','check','val',1,'pos',[20 380-3*dy 40 20],'horiz','l');
uicontrol('style','text','string',str4,'pos',[110 380-3*dy 250 20],'horiz','l')
%
h(5)=uicontrol('style','popup','string',str5a,'val',1,'pos',[20 380-4*dy 60 20],'horiz','l');
uicontrol('style','text','string',str5,'pos',[110 380-4*dy 250 20],'horiz','l')
%OK
h(6)=uicontrol('style','push','string','OK!','callback','uiresume');
%Cancel
global mquit
mquit=0;
h(7)=uicontrol('style','push','string','Cancel!','pos',get(h(6),'pos')+[80 0 0 0],'callback','global mquit,mquit=1;uiresume');
uiwait
N=str2num(get(h(1),'str'));
ismrk=get(h(2),'val');
isdyn=get(h(3),'val');
isorb=get(h(4),'val');
axunit=get(h(5),'val');
close


function [tcomp,dur,zone,isdls,issun,ismap,isgrd,isgwm,istrk,issat,istbl,israd, ...
      isdrift,istrail,isvis,elev,isssp,mquit]=menu2
hfig=dialog('Name','Appearence Settings',...
            'HandleVisibility','on');
str1b='Duration of simulation.';
str1bu={'h','m','s'};
str2b='Time zone, +-h rel. to GMT/UTC.';
str3b='Daylight saving.';
str4='Show world map.';
str5='Show grid.';
str6='Show Greenwich meridian.';
str7='Show ground-track.';
str8='Show satellite.';
str9='Show state table.';
str10='Show radius vector.';
str11='Apply perturbations.';
str12a='Show trail.';
str13='Show sun position.';
str14='Time compression factor (larger=faster).';
str16a='Visibility : Elevation [deg]';
%str16b={'off','0 deg','5 deg','10 deg','15 deg','20 deg','25 deg','30 deg'};
str17='Show visibility';
str18='Show SSP';
dy=20;
dx=40;
global h
%Duration of simulation
h{1}(1)=uicontrol('style','edit','min',1,'max',1,'string','0','pos',[20 380 40 20]);
h{1}(2)=uicontrol('style','edit','min',1,'max',1,'string','3','pos',[20+dx 380 40 20]);
h{1}(3)=uicontrol('style','edit','min',1,'max',1,'string','0','pos',[20+2*dx 380 40 20]);
uicontrol('style','text','string',str1b,'pos',[200 380 250 20],'horiz','l')
uicontrol('style','text','string',str1bu{1},'pos',[20 380+dy 40 20],'horiz','c')
uicontrol('style','text','string',str1bu{2},'pos',[20+dx 380+dy 40 20],'horiz','c')
uicontrol('style','text','string',str1bu{3},'pos',[20+2*dx 380+dy 40 20],'horiz','c')
%time zone
h{1}(4)=uicontrol('style','edit','min',1,'max',1,'string','+1','pos',[20 380-dy 40 20]);
uicontrol('style','text','string',str2b,'pos',[200 380-dy 250 20],'horiz','l')
%daylight saving
h{1}(5)=uicontrol('style','check','val',1,'pos',[20 380-2*dy 40 20]);
uicontrol('style','text','string',str3b,'pos',[200 380-2*dy 250 20],'horiz','l')
%time compression
h{1}(7)=uicontrol('style','edit','min',1,'max',1,'string','1','pos',[20 380-3*dy 60 20]);
uicontrol('style','text','string',str14,'pos',[200 380-3*dy 250 20],'horiz','l')
%%%%%%%%%%%%%%%%%
ybase=280;
%sun position
uicontrol('style','text','string',str13,'pos',[60 ybase 250 20],'horiz','l')
h{1}(6)=uicontrol('style','check','val',1,'pos',[20 ybase 40 20]);
%world map
uicontrol('style','text','string',str4,'pos',[60 ybase-dy 250 20],'horiz','l')
h{2}(1)=uicontrol('style','check','val',1,'pos',[20 ybase-dy 40 20]);
%grid
uicontrol('style','text','string',str5,'pos',[60 ybase-2*dy 250 20],'horiz','l')
h{2}(2)=uicontrol('style','check','val',1,'pos',[20 ybase-2*dy 40 20]);
%greenwich meridian
uicontrol('style','text','string',str6,'pos',[60 ybase-3*dy 250 20],'horiz','l')
h{2}(3)=uicontrol('style','check','val',1,'pos',[20 ybase-3*dy 40 20]);
%ground track
uicontrol('style','text','string',str7,'pos',[60 ybase-4*dy 250 20],'horiz','l')
h{2}(4)=uicontrol('style','check','val',1,'pos',[20 ybase-4*dy 40 20]);
%satellite
uicontrol('style','text','string',str8,'pos',[60 ybase-5*dy 250 20],'horiz','l')
h{2}(5)=uicontrol('style','check','val',1,'pos',[20 ybase-5*dy 40 20]);
%state table
uicontrol('style','text','string',str9,'pos',[60 ybase-6*dy 250 20],'horiz','l')
h{2}(6)=uicontrol('style','check','val',1,'pos',[20 ybase-6*dy 40 20]);
%radius vector
uicontrol('style','text','string',str10,'pos',[60 ybase-7*dy 250 20],'horiz','l')
h{2}(7)=uicontrol('style','check','val',1,'pos',[20 ybase-7*dy 40 20]);
%trail
uicontrol('style','text','string',str12a,'pos',[60 ybase-8*dy 250 20],'horiz','l')
h{2}(8)=uicontrol('style','check','val',1,'pos',[20 ybase-8*dy 40 20]);
%%%%%%%%%%%%%%%%%%%%%%%
xbase=200;
%perturbations
uicontrol('style','text','string',str11,'pos',[xbase+100 ybase 250 20],'horiz','l')
h{3}(1)=uicontrol('style','check','val',1,'pos',[xbase+20 ybase 40 20]);
%show visibility
uicontrol('style','text','string',str17,'pos',[xbase+100 ybase-dy 250 20],'horiz','l')
h{3}(2)=uicontrol('style','check','val',1,'pos',[xbase+20 ybase-dy 40 20]);
%visibility : elevation
h{3}(3)=uicontrol('style','edit','min',1,'max',1,'string','0','pos',[xbase+20 ybase-2*dy 40 20]);
uicontrol('style','text','string',str16a,'pos',[xbase+100 ybase-2*dy 250 20],'horiz','l')
set(h{3}(2),'callback','global h, if get(h{3}(2),''val''),set(h{3}(3),''enable'',''on''),else, set(h{3}(3),''enable'',''off''),end');
%show SSP
uicontrol('style','text','string',str18,'pos',[xbase+100 ybase-3*dy 250 20],'horiz','l')
h{3}(4)=uicontrol('style','check','val',1,'pos',[xbase+20 ybase-3*dy 40 20]);
%%%%%%%%%%%%%%%%%%%%%%%
%OK
h{4}(1)=uicontrol('style','push','string','OK!','callback','uiresume');
%Cancel
global mquit
mquit=0;
h{4}(2)=uicontrol('style','push','string','Cancel!','pos',get(h{4}(1),'pos')+[80 0 0 0],'callback','global mquit,mquit=1;uiresume');
uiwait
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dur=hms2sec(str2num(get(h{1}(1),'str')),str2num(get(h{1}(2),'str')),str2num(get(h{1}(3),'str')));
trn=0;
tsc=0;
zone=str2num(get(h{1}(4),'str'));
isdls=get(h{1}(5),'val');
issun=get(h{1}(6),'val');
tcomp=str2num(get(h{1}(7),'str'));
%%%%%%%%%%%%%%%%%%%%%%%%%
ismap=get(h{2}(1),'val');
isgrd=get(h{2}(2),'val');
isgwm=get(h{2}(3),'val');
istrk=get(h{2}(4),'val');
issat=get(h{2}(5),'val');
istbl=get(h{2}(6),'val');
israd=get(h{2}(7),'val');
istrail=get(h{2}(8),'val');
isdrift=get(h{3}(1),'val');
isvis=get(h{3}(2),'val');
elev=str2num(get(h{3}(3),'str'));
isssp=get(h{3}(4),'val');
close


function popup(ver,coastl,mu)

if exist('msattrack.key')
   load msattrack.key -mat
end
if ~exist('msattrack.key') | ~strcmp(key,'registered')
%   disp(' ')
%   disp('Thank you for using this software.')
%   disp('This is freeware but any kind of')
%   disp('donation is very welcome:')
%   disp('   -------------------')
%   disp('   -------------------')
%   disp(' ')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   d2r=pi/180;
%   mlon=coastl(:,1)*d2r;
%   mlat=coastl(:,2)*d2r;
%   popup=dialog('Name','Registering information',...
%      'HandleVisibility','on');
%   set(popup,'color',[0 0 1])
%   ppos=get(popup,'pos');
%   ptitle=uicontrol('style','text','pos',[0 ppos(4)-50 ppos(3) 50],'str',['MSatTrack ver. ' ver ]);
%   set(ptitle,'background',[0 0 1],'foreground',[1 1 1],'fontsize',20)
%   set(ptitle,'fontweight','bold','fontangle','normal')
%   paxes=axes('units','pix','pos',[300 100 200 200]);
%   [xx,yy,zz]=sphere(50);
%   fubar=surf(xx,yy,zz);
%   hold on
%   shading interp
%   set(fubar,'facecolor',[0 1 0],'edgecolor','none')
%   lighting phong
%   light('pos',[1 -.5 0])
%   [xx,yy,zz]=sph2cart(mlon,mlat,ones(size(mlon)));
%   fubar=plot3(xx,yy,zz);
%   set(fubar,'color',[0 0 0])
%   %r=orb2cart([2,.2,30,180,0,0],100,[],mu);r=[r;r(1,:)];
%   %plot3(r(:,1),r(:,2),r(:,3),'w')
%   %r=orb2cart([2,.2,30,90,0,0],100,[],mu);r=[r;r(1,:)];
%   %plot3(r(:,1),r(:,2),r(:,3),'w')
%   axis equal
%   view(3)
%   set(paxes,'color',[0 0 1],'xcol',[0 0 1],'ycol',[0 0 1],'zcol',[0 0 1])
%   str={'Thank you for showing your interrest in this software.', ...
%         'This program is shareware. A personal licence costs $30.', ...
%         'You can register this software online over the Internet', ...
%         'via the PayPal system.', ...
%         'As a registered user you will automatically get', ...
%         'future versions free of charge.'};
%   uicontrol('style','text','pos',[20 250 ppos(3)/2 100],'horiz','left','str',str,'background',[0 0 1],'foreground',[1 1 1]);
%   pl8r=uicontrol('style','push','string','Later','callback','uiresume');
%   uicontrol('style','push','string','Register!','pos',get(pl8r,'pos')+[80 0 0 0],'callback','web(''http://www.paypal.com'');uiresume');
%   uiwait
%   close
end