classdef orbit properties % Orbit properties and defaults a = 1; e = 0; i = 0; bigOmega = 0; omega = 0; T = 0; t = 0; unit_dist = 'AU'; unit_time = 'yr'; unit_angle = 'deg'; center = [0 0 0]; mu = 31556926^2 * (1/149598e6)^3 * ... % G in units of AU^3 / kg * yr^2 6.67300e-11 * 1.98911e30/(1 + 5.9736e24/1.98911e30)^2; n = 360; eccAnom = 0; meanAnom = 0; trueAnom = 0; lam0 = 0; curlypi = 0; end methods %%--------------------------------- %% CONSTRUCTOR %%================================= function orb = orbit(varargin) % ORBIT creates an object that holds orbital elements for % visualization and analysis. Simple, straightforward usage % involves creating an orbit object and then plotting it with % a call to the orbit/plot command. % % Orbit objects contain orbital element data, temporal data, % and unit data. There is also a small library of pre-defined % orbits to choose from. % % Calling syntax: % orb = orbit('elem1', val1, 'elem2', val2, ... , 'unit') % orb = orbit('object') % % Accepted orbital elements and parameters include: % Semimajor axis 'a' 'semimajor' 'axis' 'semi' % Eccentricity 'e' 'ecc' 'eccentricity' % Inclination 'i' 'inc' 'inclination' % Longitude of ascending node 'bigOmega' 'node' 'longitude' 'lon' % Argument of pericenter 'omega' 'arg' 'argument' 'peri' 'periapsis' 'pericenter' % Longitude of pericenter 'cp' 'curpi' 'curlypi' 'longperi' % Time of pericenter passage 'T' 'epoch' 'passage' % Mean longitude at epoch 'l0' 'lam0' 'lambda0' 'meanlong' 'longepoch' % Mean anomaly 'M' 'mean' 'meanAnom' 'meanAnomaly' % Eccentric anomaly 'eccAnom' 'eccentricAnom' 'eccAnomaly' 'eccentricAnomaly' % True anomaly 'f' 'true' 'trueAnom' 'trueAnomaly' % Gravitational parameter 'mu' 'gravity' % Position of central body 'center' 'focus' % (accepts either [x y z] or [r th]) % Current time 't' 'time' % % You can supply any number of those parameters to specify the % orbit; unspecified parameters will take default values. Be % aware that (1) if you supply a parameter that is % defined/derived by other parameters, be sure to supply those % first - e.g. supply bigOmega before supplying curlypi - and % (2) the orbit object stores mean anomaly values and % calculates true and eccentric anomaly from a third-order % series expansion of Kepler's Equation. % % Get orbital elements and parameters from an orbit object by % calling orb.element. Set orbital elements and parameters % either in the initial call to the constructor or by setting % orb.element = value. % % Accepted units include (just throw these into the input % stream, don't preface with 'unit' or anything): % Distance [{'AU'} 'km' 'm'] % Time [{'yr'} 's'] % Angle [{'deg'} 'rad'] % You can set orb.unit_dist, orb.unit_time, and orb.unit_angle % at any time from the command prompt or in a script. When you % do, all the orbital elements in the object will update % accordingly. % % The object library includes: % 'mercury' 'venus' 'earth' 'mars' 'jupiter' % 'saturn' 'uranus' 'neptune' 'pluto' 'moon' % (J2000 epoch, time values initialized to current time.) % % You can follow any object with other orbital elements to % modify it, e.g. orbit('earth', 't', 6) returns an orbit with % the properties of earth's but with time advanced six years % from the time of epoch. % % See also ORBIT/PLOT. % if nargin == 1 && isa(varargin{:}, 'orbit') % If given an orbit, return that orbit orb = varargin{:}; elseif nargin > 0 % Create orbit from given properties i = 1; while i <= nargin if ~isa(varargin{i}, 'char') error(['Bad argument (' num2str(i) '). Provide orbit with property/value pairs.']) end switch lower(varargin{i}) % Orbital Elements case {'a', 'axis', 'semi', 'semimaj', 'semimajor'} i = i + 1; orb.a = varargin{i}; case {'e', 'ecc', 'eccentricity'} i = i + 1; orb.e = varargin{i}; case {'i', 'inc', 'inclination'} i = i + 1; orb.i = varargin{i}; case {'bigOmega', 'node', 'longitude', 'lon'} i = i + 1; orb.bigOmega = varargin{i}; case {'omega', 'arg', 'argument', 'peri', 'periapsis', 'pericenter'} i = i + 1; orb.omega = varargin{i}; case {'T', 'epoch', 'passage'} i = i + 1; orb.T = varargin{i}; case {'t', 'time'} i = i + 1; orb.t = varargin{i}; case {'center', 'focus'} i = i + 1; orb.center = varargin{i}; case {'mu', 'gravity'} i = i + 1; orb.mu = varargin{i}; case {'m' 'mean' 'meananom' 'meananomaly'} i = i + 1; orb.meanAnom = varargin{i}; % Derived orbital elements case {'cp' 'curpi' 'curlypi' 'longperi'} i = i + 1; orb.omega = varargin{i} - orb.bigOmega; case {'l0' 'lam0' 'lambda0' 'meanlong' 'longepoch'} i = i + 1; orb.T = (varargin{i} - orb.curlypi)/orb.n; case {'eccanom' 'eccentricanom' 'eccanomaly' 'eccentricanomaly'} i = i + 1; orb.eccAnom = varargin{i}; case {'f' 'true' 'trueanom' 'trueanomaly'} i = i + 1; orb.trueAnom = varargin{i}; % Units case 'au' unit_dist = 'AU'; case {'m', 'meters'} unit_dist = 'm'; case {'km', 'kilometers'} unit_dist = 'km'; case {'s', 'sec', 'seconds'} unit_time = 's'; case {'y', 'yr', 'yrs', 'year', 'years'} unit_time = 'yr'; case {'d', 'deg', 'degree', 'degrees'} unit_angle = 'deg'; case {'rad', 'radian', 'radians'} unit_angle = 'rad'; % Orbit library case 'mercury' orb.a = 0.387098; orb.e = 0.205630; orb.i = 7.005; orb.bigOmega = 48.331; orb.omega = 29.124; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 3.3022e23/1.98911e30)^2; orb.lam0 = 252.25084; orb.t = 0; case 'venus' orb.a = 0.723332; orb.e = 0.0068; orb.i = 3.39471; orb.bigOmega = 76.67069; orb.omega = 54.85229; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 4.8685e24/1.98911e30)^2; orb.lam0 = 181.97973; orb.t = (now - datenum('01-Jan-2000'))/366; case 'earth' orb.a = 1.0000001124; orb.e = 0.016710219; orb.i = 1 + 34/60 + 43.3/3600; orb.bigOmega = 348.73936; orb.omega = 114.20783; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 5.9736e24/1.98911e30)^2; orb.lam0 = 100.46435; orb.t = (now - datenum('01-Jan-2000'))/366; case 'mars' orb.a = 1.523679; orb.e = 0.093315; orb.i = 1.850; orb.bigOmega = 49.562; orb.omega = 286.537; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 6.4185e23/1.98911e30)^2; orb.lam0 = 355.45332; orb.t = (now - datenum('01-Jan-2000'))/366; case 'jupiter' orb.a = 5.204267; orb.e = 0.048775; orb.i = 1.305; orb.bigOmega = 100.492; orb.omega = 275.066; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 1.8986e27/1.98911e30)^2; orb.lam0 = 34.40438; orb.t = (now - datenum('01-Jan-2000'))/366; case 'saturn' orb.a = 9.58201720; orb.e = 0.055723219; orb.i = 2.485240; orb.bigOmega = 113.642811; orb.omega = 336.013862; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 5.6846e26/1.98911e30)^2; orb.lam0 = 49.94432; orb.t = (now - datenum('01-Jan-2000'))/366; case 'uranus' orb.a = 19.22941195; orb.e = 0.044405586; orb.i = 0.772556; orb.bigOmega = 73.989821; orb.omega = 96.541318; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 8.681e25/1.98911e30)^2; orb.lam0 = 313.23218; orb.t = (now - datenum('01-Jan-2000'))/366; case 'neptune' orb.a = 30.10366151; orb.e = 0.011214269; orb.i = 1.767975; orb.bigOmega = 131.794310; orb.omega = 265.646853; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 1.0243e26/1.98911e30)^2; orb.lam0 = 304.88003; orb.t = (now - datenum('01-Jan-2000'))/366; case 'pluto' orb.a = 39.48168677; orb.e = 0.24880766; orb.i = 17.14175; orb.bigOmega = 110.30347; orb.omega = 113.76329; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 1.305e22/1.98911e30)^2; orb.lam0 = 238.92881; orb.t = (now - datenum('01-Jan-2000'))/366; case 'moon' orb.unit_dist = 'km'; orb.a = 384399; orb.e = 0.0549; orb.i = 5.145; orb.T = 0; orb.mu = 31556926^2 * (1/149598e6)^3 * ... 6.67300e-11 * 1.98911e30/(1 + 7.3477e22/1.98911e30)^2; orb.t = (now - datenum('01-Jan-2000'))/366; otherwise warning(['Property name ''' varargin{i} '''does not exist.']) end i = i + 1; end end end %%--------------------------------- %% DISPLAY %%================================= function display(orb) disp('ORBIT object.') disp(' Orbital elements:') disp(struct(... 'SemimajorAxis', [num2str(orb.a) ' ' orb.unit_dist], ... 'Eccentricity', num2str(orb.e), ... 'Inclination', [num2str(orb.i) ' ' orb.unit_angle], ... 'LongitudeOfAscendingNode', [num2str(orb.bigOmega) ' ' orb.unit_angle], ... 'ArgumentOfPericenter', [num2str(orb.omega) ' ' orb.unit_angle], ... 'Epoch', [num2str(orb.T) ' ' orb.unit_time] ... )); disp(' Current properties:') disp(struct(... 'Time', [num2str(orb.t) ' ' orb.unit_time], ... 'TrueAnomaly', [num2str(orb.trueAnom) ' ' orb.unit_angle], ... 'EccentricAnomaly', [num2str(orb.eccAnom) ' ' orb.unit_angle], ... 'MeanAnomaly', [num2str(orb.meanAnom) ' ' orb.unit_angle] ... )); end %%--------------------------------- %% CHANGE UNITS %%================================= function orb = set.unit_dist(orb, newDist) if ~(strcmpi(newDist,'m') || ... strcmpi(newDist,'km') || ... strcmpi(newDist,'AU')) error('Distance units should be m, km, or AU.') end if strcmpi(orb.unit_dist,'km') orb.a = orb.a * 1e3; orb.mu = orb.mu * (1e3)^3; end if strcmpi(newDist,'AU') && ... (strcmpi(orb.unit_dist,'m') || ... strcmpi(orb.unit_dist,'km')) orb.a = orb.a / 149598000000; orb.mu = orb.mu / (149598000000)^3; elseif (strcmpi(newDist,'m') || ... strcmpi(newDist,'km')) && ... strcmpi(orb.unit_dist,'AU') orb.a = orb.a * 149598000000; orb.mu = orb.mu * (149598000000)^3; end if strcmpi(newDist,'km') orb.a = orb.a / 1e3; orb.mu = orb.mu / (1e3)^3; end orb.unit_dist = newDist; end function orb = set.unit_time(orb, newTime) if ~(strcmpi(newTime,'s') || ... strcmpi(newTime,'yr')) error('Time units should be s or yr.') end if strcmpi(newTime,'yr') && ... strcmpi(orb.unit_time,'s') orb.T = orb.T / (3600*24*365); orb.t = orb.t / (3600*24*365); orb.mu = orb.mu * (3600*24*365)^2; %orb.n = orb.n * (3600*24*365); elseif strcmpi(newTime,'s') && ... strcmpi(orb.unit_time,'yr') orb.T = orb.T * (3600*24*365); orb.t = orb.t * (3600*24*365); orb.mu = orb.mu / (3600*24*365)^2; %orb.n = orb.n / (3600*24*365); end orb.unit_time = newTime; end function orb = set.unit_angle(orb, newAngle) if ~(strcmpi(newAngle,'deg') || ... strcmpi(newAngle,'rad')) error('Angle units should be deg or rad.') end if strcmpi(newAngle,'deg') && ... strcmpi(orb.unit_angle,'rad') orb.bigOmega = orb.bigOmega * 180/pi; orb.omega = orb.omega * 180/pi; orb.i = orb.i * 180/pi; orb.meanAnom = orb.meanAnom * 180/pi; elseif strcmpi(newAngle,'rad') && ... strcmpi(orb.unit_angle,'deg') orb.bigOmega = orb.bigOmega * pi/180; orb.omega = orb.omega * pi/180; orb.i = orb.i * pi/180; orb.meanAnom = orb.meanAnom * pi/180; end orb.unit_angle = newAngle; end function orb = set.unit(orb, newunit) switch lower(newunit) case {'m' 'km' 'au'} orb.unit_dist = newunit; case {'s' 'yr'} orb.unit_dist = newunit; case {'rad' 'deg'} orb.unit_angle = newunit; end end %%--------------------------------- %% SET/GET FUNCTIONS %%================================= function orb = set.t(orb, t) orb.t = t; if strcmpi(orb.unit_angle,'deg') orb.meanAnom = mod(orb.n*(t - orb.T), 360); else orb.meanAnom = mod(orb.n*(t - orb.T), 2*pi); end end % Set/get functions for parameters that are not actually stored in % the orbit object but are derived from stored parameters function n = get.n(orb) if strcmpi(orb.unit_angle, 'deg') n = 180/pi * sqrt(orb.mu / orb.a^3); else n = sqrt(orb.mu / orb.a^3); end end function orb = set.eccAnom(orb, anom) if strcmpi(orb.unit_angle,'deg') orb.meanAnom = mod(anom - orb.e*sind(anom), 360); else orb.meanAnom = mod(anom - orb.e*sin(anom), 2*pi); end end function anom = get.eccAnom(orb) if strcmpi(orb.unit_angle,'deg') anom = mod(... orb.meanAnom + orb.e*sind(orb.meanAnom) + ... orb.e^2*sind(2*orb.meanAnom)/2 + ... (3*sind(3*orb.meanAnom) - sind(orb.meanAnom))/8, ... 360); else anom = mod(... orb.meanAnom + orb.e*sin(orb.meanAnom) + ... orb.e^2*sin(2*orb.meanAnom)/2 + ... orb.e^3*(3*sin(3*orb.meanAnom) - sin(orb.meanAnom))/8, ... 2*pi); end end function orb = set.trueAnom(orb, anom) if strcmpi(orb.unit_angle,'deg') eccAnom = acosd((1 - (1-orb.e^2)/(1+orb.e*cosd(anom)))/orb.e); orb.meanAnom = mod(eccAnom - orb.e*sind(eccAnom), 360); else eccAnom = acos((1 - (1-orb.e^2)/(1+orb.e*cos(anom)))/orb.e); orb.meanAnom = mod(eccAnom - orb.e*sin(eccAnom), 2*pi); end end function anom = get.trueAnom(orb) if strcmpi(orb.unit_angle,'deg') anom = mod(... orb.meanAnom + 2*orb.e*sind(orb.meanAnom) + ... 5*orb.e^2*sind(2*orb.meanAnom)/4 + ... orb.e^3*(13*sind(3*orb.meanAnom) - 3*sind(orb.meanAnom))/12, ... 360); else anom = mod(... orb.meanAnom + 2*orb.e*sin(orb.meanAnom) + ... 5*orb.e^2*sin(2*orb.meanAnom)/4 + ... orb.e^3*(13*sin(3*orb.meanAnom) - 3*sin(orb.meanAnom))/12, ... 2*pi); end end function orb = set.curlypi(orb, curpi) if strcmpi(orb.unit_angle,'deg') orb.omega = mod(curpi - orb.bigOmega, 360); else orb.omega = mod(curpi - orb.bigOmega, 2*pi); end end function curpi = get.curlypi(orb) if strcmpi(orb.unit_angle,'deg') curpi = mod(orb.omega + orb.bigOmega, 360); else curpi = mod(orb.omega + orb.bigOmega, 2*pi); end end function orb = set.lam0(orb, lambda) if strcmpi(orb.unit_angle,'deg') orb.T = mod((lambda - orb.curlypi)/orb.n, 360); else orb.T = mod((lambda - orb.curlypi)/orb.n, 2*pi); end end function lambda = get.lam0(orb) if strcmpi(orb.unit_angle,'deg') lambda = mod(orb.T*orb.n + orb.curlypi, 360); else lambda = mod(orb.T*orb.n + orb.curlypi, 2*pi); end end %%--------------------------------- %% POSITION %%================================= function pos = position(orb) % ORBIT/POSITION takes an orbit object and returns the current % cartesian position vector of the object in its orbit. Syntax: % % pos = position(orbit) % f = orb.trueAnom; if orb.e == 1 % handle parabola specially r = orb.a .* secd((f - orb.omega)/2).^2; else r = orb.a * (1 - orb.e^2) ./ (1 + orb.e*cosd(f - orb.omega)); end if strcmpi(orb.unit_angle, 'rad') Q = expm( -cross(repmat([0 0 1]',1,3) * orb.bigOmega, eye(3)) ) * ... expm( -cross(repmat([1 0 0]',1,3) * orb.i, eye(3)) ); else Q = expm( -cross(repmat([0 0 1]',1,3) * pi/180 * orb.bigOmega, eye(3)) ) * ... expm( -cross(repmat([1 0 0]',1,3) * pi/180 * orb.i, eye(3)) ); end XYZ = Q * [r.*cosd(f); r.*sind(f); 0]; if numel(orb.center) == 3 x = orb.center(1); y = orb.center(2); z = orb.center(3); elseif numel(orb.center) == 2 if strcmpi(plotAngle, 'deg') x = orb.center(1) * cosd(orb.center(2)); y = orb.center(1) * sind(orb.center(2)); else x = orb.center(1) * cos(orb.center(2)); y = orb.center(1) * sin(orb.center(2)); end z = 0; end pos = [(XYZ(1,:)+x); (XYZ(2,:)+y); (XYZ(3,:)+z)]; end %%--------------------------------- %% PLOT %%================================= function plot(varargin) % ORBIT/PLOT takes an orbit or list of orbits and plots them % for visualization. Calling syntax: % % plot(orbit1, 'style1', orbit2, 'style2', ... ) % plot(orbit, 'property1', 'value1', 'property2', 'value2', ...) % % The string 'style' may be any conventional Matlab linestyle % string, such as 'r:' or 'bx--' just like PLOT or PLOT3. If % you use the property/value calling syntax, you must supply % only one orbit object per call to the plot command. The % property/value pairs are anything accepted by standard Matlab % line-plotting commands (see PLOT3 for details). For orbits, % there is an additional property: % % 'fill' [{'none'} 'plane' 'z'] % % which draws a transluscent shaded surface either filling in % the orbital plane or connecting the orbit to the (x,y) plane. % % The plot will be in consistent units, chosen from the first % orbit object in the calling sequence, and the plot axes will % be labelled by the unit. Note that the default Matlab axis % scalings are not to scale - you may want to use "axis equal" % to make small-inclination orbits actually look like % small-inclination orbits. % % See also ORBIT, PLOT3. idx = 1:nargin; idx = idx(cellfun('isclass', varargin, 'orbit')); plotDist = varargin{idx(1)}.unit_dist; plotTime = varargin{idx(1)}.unit_time; plotAngle = varargin{idx(1)}.unit_angle; plotargs = {}; plotplanets = {}; fill = 'n'; fillidx = 1:nargin; fillidx = fillidx(cellfun(@ischar, varargin) & ... cellfun(@(x)strcmpi(x, 'fill'), varargin)); if numel(fillidx) > 0 fillidx = fillidx(1); fill = varargin{fillidx + 1}; if fillidx + 1 < nargin varargin = {varargin{1:(fillidx-1)} varargin{(fillidx+2):end}}; else varargin = {varargin{1:(fillidx-1)}}; end end % Extract orbits from input args for i = 1:numel(idx) orb = varargin{idx(i)}; th = 0:360; if strcmpi(plotAngle, 'rad') th = th * pi/180; end orb.unit_dist = plotDist; orb.unit_time = plotTime; orb.unit_angle = plotAngle; if orb.e == 1 % handle parabola specially r = orb.a .* secd((th - orb.omega)/2).^2; else r = orb.a * (1 - orb.e^2) ./ (1 + orb.e*cosd(th - orb.omega)); end x = r.*cosd(th); y = r.*sind(th); z = zeros(size(r)); if strcmpi(plotAngle, 'rad') Q = expm( -cross(repmat([0 0 1]',1,3) * orb.bigOmega, eye(3)) ) * ... expm( -cross(repmat([1 0 0]',1,3) * orb.i, eye(3)) ); else Q = expm( -cross(repmat([0 0 1]',1,3) * pi/180 * orb.bigOmega, eye(3)) ) * ... expm( -cross(repmat([1 0 0]',1,3) * pi/180 * orb.i, eye(3)) ); end XYZ = Q * [reshape(x, [1 numel(x)]); reshape(y, [1 numel(y)]); reshape(z, [1 numel(z)]);]; if numel(orb.center) == 3 x = orb.center(1); y = orb.center(2); z = orb.center(3); elseif numel(orb.center) == 2 if strcmpi(plotAngle, 'deg') x = orb.center(1) * cosd(orb.center(2)); y = orb.center(1) * sind(orb.center(2)); else x = orb.center(1) * cos(orb.center(2)); y = orb.center(1) * sin(orb.center(2)); end z = 0; end plotargs = {plotargs{:} (XYZ(1,:)+x) (XYZ(2,:)+y) (XYZ(3,:)+z)}; if numel(varargin) > idx(i) && i == numel(idx) plotargs = {plotargs{:} varargin{(idx(i)+1):end} }; elseif numel(varargin) > idx(i) plotargs = {plotargs{:} varargin{(idx(i)+1):(idx(i+1) - 1)} }; end f = orb.trueAnom; if orb.e == 1 % handle parabola specially r = orb.a .* secd((f - orb.omega)/2).^2; else r = orb.a * (1 - orb.e^2) ./ (1 + orb.e*cosd(f - orb.omega)); end planet = Q * [r.*cosd(f); r.*sind(f); 0]; if numel(orb.center) == 3 x = orb.center(1); y = orb.center(2); z = orb.center(3); elseif numel(orb.center) == 2 if strcmpi(plotAngle, 'deg') x = orb.center(1) * cosd(orb.center(2)); y = orb.center(1) * sind(orb.center(2)); else x = orb.center(1) * cos(orb.center(2)); y = orb.center(1) * sin(orb.center(2)); end z = 0; end plotplanets = {plotplanets{:} (planet(1,:)+x) (planet(2,:)+y) (planet(3,:)+z), 'k.'}; end % Plot h = plot3(plotargs{:}); xlabel(plotDist); ylabel(plotDist); zlabel(plotDist); wasNotHeld = ~ishold; hold on % Shading if fill(1) == 'z' surf(x+[XYZ(1,:); XYZ(1,:)], y+[XYZ(2,:); XYZ(2,:)], ... [z+XYZ(3,:); zeros(size(XYZ(3,:)))], ... 'FaceAlpha', 0.2, 'EdgeColor', 'none', ... 'FaceColor', get(h, 'color')); line(x+[planet(1,:); planet(1,:)], y+[planet(2,:); planet(2,:)], ... [z+planet(3,:); zeros(size(planet(3,:)))], ... 'Color', get(h, 'color')); % Plot current position of planet plot3(plotplanets{1:3}, 'Color', get(h, 'color'), ... 'Marker', '.'); elseif fill(1) == 'p' surf(x+[XYZ(1,1:(end/2)); XYZ(1,end:(-1):(1+end/2))], ... y+[XYZ(2,1:(end/2)); XYZ(2,end:(-1):(1+end/2))], ... z+[XYZ(3,1:(end/2)); XYZ(3,end:(-1):(1+end/2))], ... 'FaceAlpha', 0.2, 'EdgeColor', 'none', ... 'FaceColor', get(h, 'color')); line(x+[planet(1,:); 0], y+[planet(2,:); 0], ... [z+planet(3,:); 0], ... 'Color', get(h, 'color')); % Plot current position of planet plot3(plotplanets{1:3}, 'Color', get(h, 'color'), ... 'Marker', '.'); else % Plot current position of planet plot3(plotplanets{:}); end if wasNotHeld hold off end end end end