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
