orbit.m

%% Calculate the orbit of planet around star over x amount of days
% Assume orbiting around the Sun. Input corresponding planet mass and
% distance from the Sun.
function [T, Y] = orbit(days)
    pos_star = [0; 0];
    mass_star = 2 * 1e30; % kg
    m = 1.8986 * 1e27; % kg
    init_velocity = [0; 13.06 * 1e3]; % m / s
    dist_from_sun = 7.779 * 1e11; % m
    init_pos = pos_star + [dist_from_sun; 0];
    [T, Y] = simulate();

    function [T, Y] = simulate()
        [T, Y] = ode45(@change, [0, days * 24 * 3600], [init_pos; init_velocity]);
        %comet(Y(:, 1), Y(:, 2))
    end

    %% To model orbit, necessary to know position of orbiting body and its velocity.
    % These change due to velocity and acceleration, respectively.
    % Assume star being orbited does not move
    function res = change(t, params)
        % Unpack variables
        Pos = params(1:2);
        Velocity = params(3:4);
        dPosdt = Velocity; % change in position is velocity
        dVelocitydt = acceleration(t, params); % change in v is acceleration
        res = [dPosdt; dVelocitydt];
    end

    function res = acceleration(t, params)
        % Unpack variables
        Pos = params(1:2);
        Velocity = params(3:4);
        % returns vector of gravity force from star on planet
        Fg = gravity_force_func(pos_star, mass_star, Pos, m);
        res = [Fg(1) / m; Fg(2) / m];
    end

end