classdef HW5_Solution_Planet < handle
    %Planet Calculates orbit and angular momentum of planet
    %To solve this problem we consider Center of Mass frame.
    %The first approximation is ignoring the mass of planet as M(Sun)>>> m(Planet)
    
    properties
        a          %major axis (AU)
        epsilon    %eccentricity
        dt         %timestep for RK4
        T          %Period (Years)
        aphelion   %(AU)
        perihelion %(AU)
        Tvec;	   %time vector
        GM=4*pi^2; %gravitational constant*Mass(Sun) in AU
        X0;        %initial conditions vector
        X;         %position and velocity matrix
        Accel;	   %accelaration
    end
    
    methods
    % find aphelion and perihelion
    function AphePeriPlanet(obj,a,epsilon)
        obj.aphelion=obj.a*(1+obj.epsilon);
        obj.perihelion=obj.a*(1-obj.epsilon); 
        obj.T=sqrt(obj.a^3) % GM=4*pi*2 in Astronomical Unit
    end
        

        
        function calcOrbit(obj)
            P=HW5_Solution_ODESolver();
	    % A is -GM for Kepler's Orbit
            P.f=@(X)HW5_Solution_CentralPotential(X,-obj.GM);
            P.dt=obj.dt;
            P.T=obj.T;
            v0=sqrt(obj.GM*(2/obj.aphelion-1/obj.a));
            P.X0=[obj.aphelion;0;0;0;v0;0];
            obj.X0=P.X0
            P.RK4();
            obj.X=P.X;
            obj.Tvec=P.Tvec;
	    % Calculate the accelaration
            r=sqrt( sum((P.X(1:3,:)).^2) );
            obj.Accel=[-obj.GM*P.X(1,:)./r.^3;...
                -obj.GM*P.X(2,:)./r.^3;...
                -obj.GM*P.X(3,:)./r.^3]
        end
        
        function plotOrbit(obj)
           figure

	   % plot Sun at origin
           plot3(0,0,0,'yo','linewidth',5)
           hold on

	   % Plot the planet 
           plot3(obj.X(1,:),obj.X(2,:),obj.X(3,:))
          
	   % Put accelaration vector at each point
%             quiver3(obj.X(1,:),obj.X(2,:),obj.X(3,:),...
%                 obj.Accel(1,:),obj.Accel(2,:),obj.Accel(3,:))
           hold off
           axis equal
           xlabel('X (AU)')
           ylabel('Y (AU)')
           zlabel('Z (AU)')
        end
        
    end
    
end

