%% TO solve this homework, I solved the planet problems in Center of Mass
%  coordinate. It means I put the Sun at origin and use reduced mass for
%  each planet is rotating around Sun. The mass of each planet is
%  negligible in comparision with the Mass of Sun.

%% Problem 1
type HW5_Solution_Planet

%% Problem 2
Earth=HW5_Solution_Planet();
Earth.a=1;             % AU
Earth.epsilon=0.0167
Earth.AphePeriPlanet()
fprintf('Aphelion of Earth is %d; Perihelion of Earth is %d\n',Earth.aphelion,Earth.perihelion)
Earth
%% Problem 3
% For this problem I choose aphelion as an initial position and the
% velocity is calculated by Velocity=sqrt(obj.GM*(2/Position-1/obj.a));
% You can choose either aphelion or perihelion position as the initial
% conditions.
Earth.T=10;      % Period of planet in Years
Earth.dt=.01;
Earth.calcOrbit()
% I use this function to show the result, it is arbitary for this problem.
% Earth.plotOrbit()

%% Problem 4

%clean up
clear all
close all
clear classes

Dt=1/365; %one day

% Mercury's orbit
Mercury=HW5_Solution_Planet();
Mercury.a=0.387;             % AU
Mercury.epsilon=0.5056;
Mercury.AphePeriPlanet()
Mercury.T=10;      % 10 Period of Earth in Years
Mercury.dt=Dt;
Mercury.calcOrbit()


% Venus's orbit
Venus=HW5_Solution_Planet();
Venus.a=0.723;             % AU
Venus.epsilon=0.0067
Venus.AphePeriPlanet()
Venus.T=-10;      % 10 Period of Earth in Years
Venus.dt=-Dt;
Venus.calcOrbit()

% Earth's orbit
Earth=HW5_Solution_Planet();
Earth.a=1;             % AU
Earth.epsilon=0.0167
Earth.AphePeriPlanet()
Earth.T=10;      % 10 Period of Earth in Years
Earth.dt=Dt;
Earth.calcOrbit()

% Mars's orbit
Mars=HW5_Solution_Planet()
Mars.a=1.52;
Mars.epsilon=0.0935;
Mars.AphePeriPlanet()
Mars.T=10;
Mars.dt=Dt;
Mars.calcOrbit()

% Jupiter's orbit
Jupiter=HW5_Solution_Planet()
Jupiter.a=5.204;
Jupiter.epsilon=0.049;
Jupiter.AphePeriPlanet()
Jupiter.T=10;
Jupiter.dt=Dt;
Jupiter.calcOrbit()

% Saturn's orbit
Saturn=HW5_Solution_Planet()
Saturn.a=9.582;
Saturn.epsilon=0.0557;
Saturn.AphePeriPlanet()
Saturn.T=10;
Saturn.dt=Dt;
Saturn.calcOrbit()

% Uranus's orbit
Uranus=HW5_Solution_Planet()
Uranus.a=19.2;
Uranus.epsilon=0.0472;
Uranus.AphePeriPlanet()
Uranus.T=10;
Uranus.dt=Dt;
Uranus.calcOrbit()

% Neptune's orbit
Neptune=HW5_Solution_Planet()
Neptune.a=30.01;
Neptune.epsilon=0.0087;
Neptune.AphePeriPlanet()
Neptune.T=10;
Neptune.dt=Dt;
Neptune.calcOrbit()

% Pluto's orbit
Pluto=HW5_Solution_Planet()
Pluto.a=39.26;
Pluto.epsilon=0.245;
Pluto.AphePeriPlanet()
Pluto.T=10;
Pluto.dt=Dt;
Pluto.calcOrbit()


% %Setup movie
filename='KeplerOrbits.mp4'
VidObj=VideoWriter(filename,'MPEG-4')

%Set FrameRate property
VidObj.FrameRate=30; %can't do real time so make it video rate

% %Open video file for writing
VidObj.open();
% Make a bigger figure
%fh=figure('units','normalized','outerposition',[0,0,1,1]);  %this may be too big for video writer on large displays
fh=figure('units','normalized','outerposition',[0,0,1/2,1/2]);  %Use this instead

xmax=50

%display only every Nth frame
N=3;
for ii=1:N:length(Earth.Tvec) % Choose the the smallest Planet.Period to show the movie
    % In this case Earth.Period < Mars.Period
    figure(fh)
    hold off
    
    % As I solved it in Center of Mass coordinate, it means I put the Sun
    % at origin and use reduced mass for each planet is rotating around
    plot(0,0,'yo','linewidth',5)
    hold on
    
    %Plot the position of the Planets - I'll do in 2D 
    plot(Mercury.X(1,ii),Mercury.X(2,ii),'Marker','o','Color',[1 0 0],'linewidth',2)
    plot(Venus.X(1,ii),Venus.X(2,ii),'Marker','o','Color',[0 0 0.25],'linewidth',2)
    plot(Earth.X(1,ii),Earth.X(2,ii),'Marker','o','Color',[0.5 1 0.25],'linewidth',2)
    plot(Mars.X(1,ii),Mars.X(2,ii),'Marker','o','Color',[1 0.75 0.25],'linewidth',2)
    plot(Jupiter.X(1,ii),Jupiter.X(2,ii),'Marker','o','Color',[0.5 0.5 0.5],'linewidth',2)
    plot(Saturn.X(1,ii),Saturn.X(2,ii),'Marker','o','Color',[0.5 0.5 0.75],'linewidth',2)
    plot(Uranus.X(1,ii),Uranus.X(2,ii),'Marker','o','Color',[0.5 0.5 1],'linewidth',2)
    plot(Neptune.X(1,ii),Neptune.X(2,ii),'Marker','o','Color',[0.5 1 0.5],'linewidth',2)
    plot(Pluto.X(1,ii),Pluto.X(2,ii),'Marker','o','Color',[0.75 0.5 0.5],'linewidth',2)
    
    %Draw velocity vectors. This part is not necessary for homework solution.
    quiver(Mercury.X(1,ii),Mercury.X(2,ii),Mercury.X(4,ii),Mercury.X(5,ii),0,'b','linewidth',.5)
    quiver(Venus.X(1,ii),Venus.X(2,ii),Venus.X(4,ii),Venus.X(5,ii),0,'b','linewidth',.5)
    quiver(Earth.X(1,ii),Earth.X(2,ii),Earth.X(4,ii),Earth.X(5,ii),0,'b','linewidth',.5)
    quiver(Mars.X(1,ii),Mars.X(2,ii),Mars.X(4,ii),Mars.X(5,ii),0,'b','linewidth',.5)
    quiver(Jupiter.X(1,ii),Jupiter.X(2,ii),Jupiter.X(4,ii),Jupiter.X(5,ii),0,'b','linewidth',.5)
    quiver(Saturn.X(1,ii),Saturn.X(2,ii),Saturn.X(4,ii),Saturn.X(5,ii),0,'b','linewidth',.5)
    quiver(Uranus.X(1,ii),Uranus.X(2,ii),Uranus.X(4,ii),Uranus.X(5,ii),0,'b','linewidth',.5)
    quiver(Neptune.X(1,ii),Neptune.X(2,ii),Neptune.X(4,ii),Neptune.X(5,ii),0,'b','linewidth',.5)
    quiver(Pluto.X(1,ii),Pluto.X(2,ii),Pluto.X(4,ii),Pluto.X(5,ii),0,'b','linewidth',.5)
    
    % Show time
    s=sprintf('%4.2g years',Earth.Tvec(1,ii));
    text(10,-8,0,s);
    axis equal
    axis([-12 xmax -10 12])
    
    hold off    %so new plot will won't add to this one
    
    xlabel('X (AU)')
    ylabel('Y (AU)')
    
    title('Solar System')
    legend('Sun','Mercury','Venus','Earth','Mars','Jupiter','Saturn','Uranus','Neptune','Pluto')
    
    %now add frame to movie
    currFrame = getframe();
    VidObj.writeVideo(currFrame);
    
end

%close video file
VidObj.close()




%% Problem 5
% To simplify this problem, I assume the particles collide the X-Y plane
% only in the x-axis. Look at X0.

clear all
close all

Z1=2;   %Charge on alpha
Z2=79;  %Charge on Gold

m_a=6.64*10^-27;    %mass of alpha in kg
e=1.6*10^-19;       %C
k=8.9*10^9;         %kg m^3 s^-2 C^-2
A=k*Z1*Z2*e^2/m_a;  %m^3/s^2
V_0=2*10^7;         %m/s

%Setup ODE solver
P=HW5_Solution_ODESolver();

%Use the central potential for scattering
P.f=@(X)HW5_Solution_CentralPotential(X,A)

%guess some intial conditions and parameters to give good scattering plots
Z_0 = -1e-11;
P.T=-2*Z_0/V_0;
P.dt=P.T/1000;

%This is the vector of impact parameters to be tested. 
%Let's use log spaced values to sample several orders of magnitude
X0_vec=logspace(-15,-11);

%Make figure and show trajectories
figure
plot(0,0,'yo','linewidth',2)
hold on
for ii=1:length(X0_vec)
    X0=[X0_vec(ii) 0 Z_0 0 0 V_0];
    P.X0=X0;
    P.RK4()
    plot(P.X(1,:),P.X(3,:))
    xlabel('X (m)')
    ylabel('Y (m)')
    zlabel('Z (m)')
    
    %this is the scattering angle calculated from the final velocity
    theta(ii)=atan2(P.X(4,end),P.X(6,end)); 
    
    b(ii)=X0_vec(ii);  %record the impact parameter. 
end
axis equal

%compare to theory. This is Scattering angle versus impact parameter
RS=@(b)2*atan(k*Z1*Z2*e^2./(m_a*V_0^2*b))

%plot the simulation and theory together. 
figure
semilogx(b,theta,'ro')
hold on
semilogx(b,RS(b),'k');
legend('Simulated','Theory')
xlabel('Impact Parameter (m)')
ylabel('Scattering Angle (radians)')



