Code covered by the BSD License

# orbits - plot orbits around Earth in an interactive manner.

by

### Michael (view profile)

23 Jun 2011 (Updated )

This GUI lets you plot orbits around Earth interactively.

orbits(action);
```function [] = orbits(action);
% This function plots an Earth orbit based on
% the following orbtial elements:
% perigee = altitude at perigee (miles)
% apogee = altitude at apogee (miles)
% i = inclination of orbital plane to equatorial plane (degrees)
% Omega = argument of ascending node (degrees)
% omega = argument of periapsis (degrees)
%
% Written with Matlab 5.3 by:
% Michael Hanchak  (mhanchak AT yahoo DOT com)
% Dayton, Ohio
% August 9, 2001
%
% Reference Text:
% Fundamentals of Astrodynamics
%
% Users of this function assume all risk in using the calculations
% contained within. This function is not to be distributed without
% the author's full consent.
%
% For best results, requires "coast.mat" and "topo.mat"

% standard switchyard programming logic
if nargin == 0;
action = 'build';
end

colors = [.6 .7 .8];
colors2 = 'k';
switch action

% build GUI
case 'build'

strings = {'w','W','i','Apogee','Perigee'};
tags = {'o','O','i','alt_a','alt_p'};
values = {'45','45','30','345','115'};
strings2 = {'Plot Orbit','Clear Orbits','Center Earth',...
'Zoom All','Flyby','Help','Toggle Earth','Quit'};
callbacks = {'orbits(''plot'')','orbits(''clear'')',...
'camva(15);view(120,30);camlookat(findobj(''tag'',''earth''));',...
'camva(15);camlookat','orbits(''flyby'')','help orbits',...
'orbits(''earth'')','close(gcf)'};

if isempty(findobj('tag','orbits'))
www = figure('tag','orbits');
else
www = findobj('tag','orbits');
figure(www);
clf
end
set(www,'position',[25   75   600   500],'color',colors2);

for kk = 1:length(strings),
ppp = uicontrol('Units','pixels','Position',[2 (20*(kk-1)+16) 50 20],...
'String',strings{kk},'style','text','backgroundcolor',colors2,...
'foregroundcolor','w','fontsize',9);
uicontrol('Units','pixels','Position',[55 (20*(kk-1)+16) 50 20],...
'tag',tags{kk},'style','edit','string',values{kk},'backgroundcolor',[1 1 1]);
if kk <=2
set(ppp,'fontname','symbol','fontsize',12);
end
end
for kk =1:length(strings2),
uicontrol('Units','pixels','Position',[10 (25*(kk-1)+150) 70 22],...
'string',strings2{kk},'callback',callbacks{kk},'backgroundcolor',colors);
end
uicontrol('Units','pixels','Position',[10 370 100 80],...
'string',['Left click and drag Earth for dynamic viewing.',...
' Right click for zooming. Double-Click to center. '],...
'style','text','backgroundcolor',colors2,'fontsize',8,...
'HorizontalAlignment','left','foregroundcolor','w');
axes('position',[.2 .05 .75 .9],'units','normalized')
axis equal
axis off
axis vis3d
hold on

% Plot reference frame axes
h3 = plot3([0 0],[0 0],[0 rad+500],'color',[0 0 .8]);
set([h1 h2 h3],'linewidth',4);

view(120,30);
camva(15);
rot3d;
data.simple = 1;
data.handles = [];
data.earthhandle = [];
set(gcf,'userdata',data);

orbits('earth') % call routine to draw earth
camlookat(findobj('tag','earth'));

case 'earth'
% determine level of graphics

data = get(gcf,'userdata');
simple = data.simple;
delete(data.earthhandle);
data.earthhandle = [];
% Plot the earth (texture map or simple)
if simple == 0 & ~isempty(which('topo.mat'))

[X,Y,Z] = sphere(50);
topo = [topo(:,181:360) topo(:,1:180)];
mat.dull.AmbientStrength = 0.4;
mat.dull.DiffuseStrength = .6;
mat.dull.SpecularColorReflectance = .5;
mat.dull.backfacelighting = 'reverselit';
mat.dull.SpecularExponent = 20;
mat.dull.SpecularStrength = .8;
mat.dull, ...
'FaceColor','texturemap',...
'EdgeColor','none',...
'FaceLighting','phong',...
'Cdata',topo,'tag','earth');
colormap(topomap1)
%light('position',rad*[-10 -10 -10], 'color', [.6 .2 .2]);
%set(gcf,'renderer','opengl');
data.simple = 1;
else

[X,Y,Z] = sphere(24);
%set(data.earthhandle,'tag','earth','facecolor',[.6 .7 .9],'edgecolor',[1 1 1]);
set(data.earthhandle(1),'tag','earth','facecolor',[0 0 1],'edgecolor',[.3 .3 1]);
if ~isempty(which('coast.mat'))

ncst = ncst *pi/180;
all =zeros(length(ncst),3);
for j = 1:length(ncst)
theta = ncst(j,1);
phi = ncst(j,2);
all(j,:) = [cos(theta)*cos(phi),...
sin(theta)*cos(phi),...
-sin(phi)];
end

set(data.earthhandle(2),'color',[0 .9 0]);

end
data.simple = 0;
end

set(gcf,'userdata',data);

% plot orbits
case 'plot'

data = get(gcf,'userdata');
handles = data.handles;

% get orbital elements from GUI
alt_p = str2num(get(findobj('tag','alt_p'),'string'));
alt_a = str2num(get(findobj('tag','alt_a'),'string'));
inc = str2num(get(findobj('tag','i'),'string'));
Omega = str2num(get(findobj('tag','O'),'string'));
omega = str2num(get(findobj('tag','o'),'string'));

% check for correctness of input data
if alt_p > alt_a
error1 = errordlg('Perigee must be smaller than apogee');
waitfor(error1);
else

% Orbital elements
a = (alt_p + alt_a + 2*rad)/2;
c = a - alt_p - rad;
e = c/a;
p = a*(1 - e^2);

th = linspace(0,2*pi,200);
r = p./(1 + e*cos(th));
xx = r.*cos(th);
yy = r.*sin(th);

% Coordinate Transformations
ZZ = [cos(Omega) -sin(Omega) 0;
sin(Omega) cos(Omega) 0;
0 0 1];
XX = [1 0 0;
0 cos(inc) -sin(inc);
0 sin(inc) cos(inc)];
ZZ2 = [cos(omega) -sin(omega) 0;
sin(omega) cos(omega) 0;
0 0 1];

% actual plot
vec = ZZ*XX*ZZ2*[xx;yy;zeros(1,length(xx))];
h1 = plot3(vec(1,:),vec(2,:),vec(3,:));
set([h1],'linewidth',1,'color',[1 1 1]);

% line of ascending node
h2 = plot([0 vec1(1)],[0 vec1(2)],'r-');
% line of periapsis
h3 = plot3([0 vec2(1)],[0 vec2(2)],[0 vec2(3)],'color',...
[1 1 0]);
% line of inclination
h4 = plot3([0 vec3(1)],[0 vec3(2)],[0 vec3(3)],...
'color',[0 0 .8]);
set([h2 h3 h4],'linewidth',2);

data.handles = [data.handles h1 h2 h3 h4];
set(gcf,'userdata',data);
%camlookat;
end

case 'clear'
data = get(gcf,'userdata');
handles = data.handles;
for rr = 1:length(handles)
if ishandle(handles(rr))
delete(handles(rr));
end
end
data.handles = [];
set(gcf,'userdata',data);

case 'flyby'
camlookat(findobj('tag','earth'));
camva(15);
camup([0 0 1]);
for x = -300000:1000:60000
campos([x,30000,30000])
drawnow
end

% here are the recursive callbacks for the rotation of the axes
case 'rot'
rot3d('rot');
case 'down'
rot3d('down');
case 'up'
rot3d('up');
case 'zoom'
rot3d('zoom');

end

% this function below allows for dynamics "click and drag"
%  rotation of the plot.
function rot3d(huh)

if nargin<1
set(gcf,'WindowButtonDownFcn','orbits(''down'')');
set(gcf,'WindowButtonUpFcn','orbits(''up'')');
set(gcf,'WindowButtonMotionFcn','');
else

switch huh
case 'down'
if strcmp(get(gcf,'SelectionType'),'normal')
set(gcf,'WindowButtonMotionFcn','orbits(''rot'')');
elseif strcmp(get(gcf,'SelectionType'),'alt')
set(gcf,'WindowButtonMotionFcn','orbits(''zoom'')');
elseif strcmp(get(gcf,'SelectionType'),'open') %center point
temp1 = get(gca,'currentpoint');
temp1 = (temp1(1,:) + temp1(2,:))/2; % average points
%newvec = temp1 - campos
%oldvec = camtarget-campos
%dir = cross(oldvec,newvec)
%ang = atan2(norm(dir),dot(oldvec,newvec))+2
%campan(-ang,0,'data',dir);
camtarget([temp1]);
end
rdata.oldpt = get(0,'PointerLocation');
set(gca,'userdata',rdata);
case 'up'
set(gcf,'WindowButtonMotionFcn','');
case 'rot'
rdata = get(gca,'userdata');
oldpt = rdata.oldpt;
newpt = get(0,'PointerLocation');
dx = (newpt(1) - oldpt(1))*.5;
dy = (newpt(2) - oldpt(2))*.5;

%direction = [0 0 1];
%coordsys  = 'camera';
%pos  = get(gca,'cameraposition' );
%targ = get(gca,'cameratarget'   );
%dar  = get(gca,'dataaspectratio');
%up   = get(gca,'cameraupvector' );
%[newPos newUp] = camrotate(pos,targ,dar,up,-dx,-dy,coordsys,direction);
%set(gca,'cameraposition', newPos, 'cameraupvector', newUp);

camorbit(gca,-dx,-dy,'camera');

rdata.oldpt = newpt;
set(gca,'userdata',rdata);
case 'zoom'
rdata = get(gca,'userdata');
oldpt = rdata.oldpt;
newpt = get(0,'PointerLocation');
dy = (newpt(2) - oldpt(2))/abs(oldpt(2));
camzoom(gca,1+dy)
rdata.oldpt = newpt;
set(gca,'userdata',rdata)
end

end
```