function orbits(action)

%ORBITS  GUI demonstrating orbit calculations and map interaction
%
%  Synopsis:  orbits
%
%  ORBITS produces an interactive window allowing the analysis
%  of various satellite orbits and their resulting rotating
%  earth ground traces.  This example illustrates the utility
%  of a map in a numerical application.  The satellite queries
%  the map at each point and displays its current location.
%  This interaction happens within the simulation and represents
%  how a map can be used as data within an application.  This
%  script is executed by entering ORBITS.

%  Copyright 1996-1998 by Systems Planning and Analysis, Inc. and The MathWorks, Inc.
% $Revision: 1.13 $
%  Written by:  E. Byrns, E. Brown



if nargin == 0;   action = 'initialize';   end


switch action
case 'initialize'

%  Construct the interface window

     objhndls = orbpanl;

%  Load the world map and display

     load worldmtx
     countries = strvcat(nations(:).name);
	 colors = [0    0    0
	           0    0.7500    1.0000
		       0    0    0
			   ones(size(countries,1)-3,1)*[0    0    0]];

     axes(objhndls(1))
	 axesm behrmann; colormap(colors);
	 axis([ -2.7207    2.7207   -1.1547    1.1547]) % numbers hardcoded for 
	                                                % behrmann; axis tight 
													% causes orbit line to be 
													% lost upon stop
	 set(objhndls(1),'ButtonDownFcn','')
	 meshm(map,maplegend,[30 50],'ButtonDownFcn','');
	 
%  Display initial orbit

     apogee      = get(objhndls(5),'Value');
     perigee     = get(objhndls(8),'Value');
     inclination = get(objhndls(11),'Value');

     [newcoord,lat,long] = orbcalcs('nonrotating',apogee,perigee,...
                                     inclination,(0:360)*pi/180);

     hold on;  axes(objhndls(1));
	 objhndls(17) = plot3m(lat,long,1.01*ones(size(lat)),'r.',...
	                       'EraseMode','xor','MarkerSize',6,...
						   'ButtonDownFcn','');

     axes(objhndls(2));
     objhndls(18) = plot3(newcoord(1,:), newcoord(2,:), ...
               newcoord(3,:), 'r-',...
			   'EraseMode','normal');
     hold on;

     plot3(0,0,0,'b.', 'MarkerSize',60);  hold off
     axis([-6 6 -6 6 -6 6]);  axis('off'); drawnow
     
%  Save the object handles and map data

     set(objhndls(3),'UserData',map);
     set(objhndls(4),'UserData',maplegend);
     set(objhndls(5),'UserData',countries);
     set(gcf,'UserData',objhndls,'Visible','on');


case  'run'    % Run button pushed

    objhndls = get(gcf,'UserData');
	qryhndls = get(objhndls(3),'UserData');

	set(objhndls(5),  'Enable', 'off')
	set(objhndls(8),  'Enable', 'off')
	set(objhndls(11), 'Enable', 'off')
	set(objhndls(14), 'Enable', 'off')
	set(objhndls(15), 'Enable', 'on')
	set(objhndls(16), 'Enable', 'off')

%  Get necessary data

    apogee      = get(objhndls(5),'Value');
    perigee     = get(objhndls(8),'Value');
    inclination = get(objhndls(11),'Value');

     map       = get(objhndls(3),'UserData');
     maplegend = get(objhndls(4),'UserData');
	 countries = get(objhndls(5),'UserData');

%  Loop until the user stops.  Display position on both
%  the map and 3D graph.

	loop = 0;   tottime = 0;   theta = 0;
	set(objhndls(2),'UserData',1);  hold on;

	while get(objhndls(2),'UserData')

%  Compute the new position

		loop = loop + 1;
        [newcoord,lat,long,theta,tottime] = orbcalcs('simulate',...
                    apogee,perigee,inclination,theta,tottime);

%  Update the current lat/long display

        set(objhndls(21),'String',num2str(fix(lat)));
        set(objhndls(22),'String',num2str(fix(long)));

%  Determine the map cell position.  Update display if
%  not on a border.

        cellval = ltln2val(map,maplegend,lat,long);
        if cellval ~= 1
		        set(objhndls(13),'String',countries(cellval,:))
		end

%  Update the display on the map and delete the old display
%  on the 3D graph.

        if loop > 1
              set(hndl1,'MarkerSize',10);
			  set(hndl2,'Xdata',newcoord(1),'Ydata',newcoord(2),...
			            'Zdata',newcoord(3) )
        else
              axes(objhndls(2));
              hndl2 = line(newcoord(1),newcoord(2),newcoord(3),...
				     'MarkerFaceColor','red','MarkerEdgeColor','red',...
                     'Marker','.',...
			         'EraseMode','xor', 'MarkerSize',20);
		end


%  Display new point on the map

        axes(objhndls(1));
        hndl1 = plot3m(lat,long,1.01*ones(size(lat)),'r.',...
               'MarkerSize',20, 'EraseMode','xor','ButtonDownFcn','');

		drawnow;

	end

%  At this point, the user has stopped the simulation.
%  Set the slider and buttons to be enabled again.

	set(objhndls(5),  'Enable', 'on')
	set(objhndls(8),  'Enable', 'on')
	set(objhndls(11), 'Enable', 'on')
	set(objhndls(14), 'Enable', 'on')
	set(objhndls(15), 'Enable', 'off')
	set(objhndls(16), 'Enable', 'on')

%  Delete the satellite display

%  Clear old orbit traces.

	delete([hndl1 hndl2]);
    oldorbits = get(objhndls(1),'Children');
	if length(oldorbits) > 6
	      delete(oldorbits(1:length(oldorbits)-6));  clear oldorbits
	end


case 'sliders'     % User has activated a slider

%  Retrieve data and round off the slider value

    objhndls = get(gcf,'UserData');

    apogee   = get(objhndls(5),'Value');
	apogee   = fix(apogee/100)*100;

    perigee    = get(objhndls(8),'Value');
	perigee    = fix(perigee/100)*100;

%  Set the slider value displays and sync the slider limits

	set(objhndls(4),'String',num2str(apogee));
	set(objhndls(5),'Value',max(apogee,110), 'Min', max(perigee,110));

	set(objhndls(7),'String',num2str(perigee));
	set(objhndls(8),'Value',perigee, 'Max', max(min(apogee,40000),110));

%  Get and display the inclination slider information

    inclination = get(objhndls(11),'Value');
	inclination = fix(inclination);
	set(objhndls(10),'String',num2str(inclination));
	set(objhndls(11),'Value',inclination);

%  Display new orbit

    [newcoord,lat,long] = orbcalcs('nonrotating',apogee,perigee,...
                                    inclination,(0:360)*pi/180);

    datastruct = gcm(objhndls(1));
    [x,y,z] = mfwdtran(datastruct,lat,long,0,'none');

    set(objhndls(17),'Xdata',x,'Ydata',y,'Zdata',1.01*ones(size(x)));

    set(objhndls(18),'Xdata',newcoord(1,:), ...
	    'Ydata',newcoord(2,:), 'Zdata', newcoord(3,:));


case 'stop'   % Stop button selected

%  Set the while loop test variable to false

	 objhndls = get(gcf,'UserData');
	 set(objhndls(2),'UserData',0);

end


%*************************************************************************
%*************************************************************************
%*************************************************************************


function [newcoord,lat,long,theta,tottime] = orbcalcs(action,...
              apogee,perigee,inclination,theta,tottime)

%ORBCALCS  Calculates orbital position (used with ORBITS)
%
%  ORBCALCS determines the orbital position of a satellite
%  in both 3 dimensional space and relative to the Earth
%  latitude and longitude.  These calculations are performed
%  for both a non-rotating earth and a rotating earth.  The
%  function uses conic section properties to determine
%  the satellite position at a specified time in its orbit.

%  The conic section approach still requires that the dependent
%  variable be time and not angle so that the simulation
%  display has the appropriate velocity characteristics.

%  Copyright 1996-1998 by Systems Planning and Analysis, Inc. and The MathWorks, Inc.
% $Revision: 1.13 $
%  Written by:  E. Byrns, E. Brown


mu       = 3.986E5;
radearth = 6371;
omega    = 7.2921158553e-5;


e         = abs(apogee-perigee)/(apogee + perigee + 2*radearth);
semimajor = (apogee+perigee)/2 + radearth;
totarea   = pi*semimajor*sqrt(semimajor^2 * (1-e^2));
period    = 2*pi*sqrt(semimajor^3/mu);

rotation = [cos(inclination*pi/180)  0 -sin(inclination*pi/180)
                  0           1        0
            sin(inclination*pi/180)  0  cos(inclination*pi/180)];


if strcmp(action,'nonrotating')

     radius   = semimajor*(1-e^2) ./ (1+e*cos(theta));
     oldcoord = [radius.*cos(theta); radius.*sin(theta); zeros(size(theta))];
     newcoord = (rotation * oldcoord)/radearth;

     long = atan2(newcoord(2,:),newcoord(1,:));
     long = long*180/pi;

     lat  = atan(newcoord(3,:)./sqrt(newcoord(1,:).^2 + newcoord(2,:).^2));
     lat  = lat*180/pi;

elseif strcmp(action,'simulate')

     deltime = period/50;
     delarea = deltime * totarea / period;
     tottime = tottime + deltime;

     radius   = semimajor*(1-e^2) ./ (1+e*cos(theta));
     delta    = 2*delarea/radius^2;
	 theta    = theta + delta;

     radius   = semimajor*(1-e^2) ./ (1+e*cos(theta));

     oldcoord = [radius.*cos(theta); radius.*sin(theta); zeros(size(theta))];
     newcoord = (rotation * oldcoord)/radearth;

     long = atan2(newcoord(2,:),newcoord(1,:));
     long = (long - omega*tottime);
	 long = atan2(sin(long),cos(long))*180/pi;

     lat  = atan(newcoord(3,:)./sqrt(newcoord(1,:).^2 + newcoord(2,:).^2));
     lat  = lat*180/pi;
end


%*************************************************************************
%*************************************************************************
%*************************************************************************


function savehndl = orbpanl

%ORBPANL  Creates the interface objects for the ORBITS GUI
%
%  The data for the operation of this demo is stored in
%  various user data positions in the interface objects.

%  Copyright 1996-1998 by Systems Planning and Analysis, Inc. and The MathWorks, Inc.
% $Revision: 1.13 $
%  Written by:  E. Byrns, E. Brown


limits = [100  5000 40000];    %  Apogee, perigee slider data
degrees = [-180 180   30];     %  Inclination slider data

%  Compute the Pixel and Font Scaling Factors so
%  GUI figure windows and fonts look OK across all platforms

PixelFactor = guifactm('pixels');
FontScaling =  guifactm('fonts');


figure('Visible','off', 'Units','Points',...
       'Resize','on',   'Name','Satellite Orbits', 'Numbertitle','off',...
       'Renderer','zbuffer');
colordef(gcf,'white')
figclr = get(gcf,'Color');

%  Map Axes

savehndl(1) = axes('Units','Normalized', 'Position',[0.05  0.45  0.90  0.50]);

%  Space 3D Axes

savehndl(2) = axes('Units','Normalized', 'Position',[0.45  0.05  0.50  0.40],...
            'View',[-37.5 30], 'Visible','off', 'DrawMode','Fast');

%  Apogee Slider Object

savehndl(3) = uicontrol(gcf, 'Style','Text', 'String','Apogee (km)',...
          'Units','Normalized', 'Position',[0.05  0.35  0.13  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','black');

savehndl(4) = uicontrol(gcf, 'Style','Text', 'String',num2str(limits(2)),...
          'Units','Normalized',  'Position',[0.19  0.35  0.12  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','black');

savehndl(5) = uicontrol(gcf, 'Style','slider',...
			 'Units','Normalized', 'Position',[0.05  0.31  0.18  0.03], ...
			 'Min',limits(2), 'Max',limits(3)+10, 'Value',limits(2), ...
			 'CallBack','orbits(''sliders'');');

%  Perigee Slider Object

savehndl(6) = uicontrol(gcf, 'Style','Text', 'String','Perigee (km)',...
          'Units','Normalized', 'Position',[0.05  0.25  0.15  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','black');

savehndl(7) = uicontrol(gcf, 'Style','Text', 'String',num2str(limits(2)),...
          'Units','Normalized', 'Position',[0.19  0.25  0.12  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','black');

savehndl(8) = uicontrol(gcf, 'Style','slider',...
			 'Units','Normalized', 'Position',[0.05  0.21  0.18  0.03], ...
			 'Min',limits(1), 'Max',limits(2), 'Value',limits(2), ...
			 'CallBack','orbits(''sliders'');');

%  Inclination Slider Object

savehndl(9) = uicontrol(gcf, 'Style','Text', 'String','Inclination (deg)',...
          'Units','Normalized',  'Position',[0.05  0.15  0.18  0.04],...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','black');

savehndl(10) = uicontrol(gcf, 'Style','Text', 'String',num2str(degrees(3)),...
          'Units','Normalized', 'Position',[0.21  0.15  0.08  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','black');

savehndl(11) = uicontrol(gcf, 'Style','slider',...
			 'Units','Normalized', 'Position',[0.05  0.11  0.18  0.03], ...
			 'Min',degrees(1), 'Max',degrees(2), 'Value',degrees(3), ...
			 'CallBack','orbits(''sliders'');');

%  Text box for Location Title

savehndl(12) = uicontrol(gcf, 'Style','Text', 'String','Current Location:', ...
		  'Units','Normalized', 'Position',[0.05  0.05  0.17  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','red' );

%  Text box for Current Location

savehndl(13) = uicontrol(gcf, 'Style','Text', 'String','', ...
		  'Units','Normalized', 'Position',[0.22  0.05  0.29  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','red' );

%  Text box for Lat/Long Title

savehndl(20) = uicontrol(gcf, 'Style','Text', 'String','Lat/Long:', ...
		  'Units','Normalized', 'Position',[0.05  0.01  0.09  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','red' );

%  Text box for Latitude

savehndl(21) = uicontrol(gcf, 'Style','Text', 'String','', ...
		  'Units','Normalized', 'Position',[0.15  0.01  0.08  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','red' );

%  Text box for Longitude

savehndl(22) = uicontrol(gcf, 'Style','Text', 'String','', ...
		  'Units','Normalized', 'Position',[0.24  0.01  0.09  0.04], ...
		  'HorizontalAlignment','left',...
		  'BackgroundColor',figclr, 'ForegroundColor','red' );

% Run button

savehndl(14) = uicontrol(gcf, 'Style','push', 'String','Run',...
			'Units','Normalized', 'Position',[0.29  0.31  0.10  0.06], ...
			'Callback','orbits(''run'');', 'Interruptible','on');

% Stop button

savehndl(15) = uicontrol(gcf, 'Style','push', 'String','Stop',...
			'Units','Normalized', 'Position',[0.29  0.24  0.10  0.06], ...
			'Callback','orbits(''stop'');', 'Enable','off');

% Exit button

savehndl(16) = uicontrol(gcf, 'Style','push', 'String','Exit',...
			'Units','Normalized', 'Position',[0.29  0.17  0.10  0.06], ...
			'Callback','close');


%*************************************************************************
%*************************************************************************
%*************************************************************************
