%% Problem 1
clc
clear all

type HW4_Solution_ODESolver.m

%% Problem 2

% This is exponential decay:  N(t)=A*exp(-t/tau).  
% t_half = ln(2)*tau
t_half=30;          %days
tau=t_half/log(2);   %days

%our differential equation is N'(t)=-1/tau*N(t)
N_p=@(N)-1/tau*N;
X0=100;

%Create object and setup problem
ODE=HW4_Solution_ODESolver();
ODE.f=N_p;
ODE.X0=X0;  %percent
ODE.dt=.1;  %days
ODE.T=365;  %days
%Run Euler method
ODE.Euler();
fprintf('Amount remaining after one year: %g percent\n',ODE.X(end))

%% Problem 3
%Its a good idea to visualize result
ODE.plotX1('Time (days)','Percent Remaining');

%% Problem 4
% A mass of 1 kg is attached to a spring that has spring constant 10 N/m.  At
% time t=0, the mass is let go from 1 m above it's equilibrium position.
% Use Euler's Method and a dt=.001 s to find and plot the position over time
% for 100 s.  Plot the analytical solution on the same plot using a
% different color. 
m=1;        %kg
k=10;       %N/m

%I'm taking equilibrium position at y=0;

%our differential equation is Y'(t)=-k/m*Y(t)
y_p=@(y)[y(2);-k/m*y(1)];
y0=1;   %m
v0=0;   %m/s
X0=[y0; v0];

%Create object and setup problem
ODESpring=HW4_Solution_ODESolver();
ODESpring.f=y_p;
ODESpring.X0=X0;  %initial position and velocity
ODESpring.dt=.001;  %s
ODESpring.T=100;  %s
%Run Euler method
ODESpring.Euler();
%plot:
fh=ODESpring.plotX1('Time (s)','Y (m)')

%Add analytical solution to plot
figure(fh);     %in case this figure is not current
hold on;        %so we don't erase old plot
y_A=@(t)y0*cos(sqrt(k/m)*t)     %Analytical Solution
plot(ODESpring.Tvec,y_A(ODESpring.Tvec),'r--')
legend('ODE Solver','Analytical')
title('Euler')


%% Problem 5

%Re-Run using RK4 method
ODESpring.RK4();
%plot:
fh=ODESpring.plotX1('Time (s)','Y (m)')

%Add analytical solution to plot
figure(fh);     %in case this figure is not current
hold on;        %so we don't erase old plot
y_A=@(t)y0*cos(sqrt(k/m)*t)     %Analytical Solution
plot(ODESpring.Tvec,y_A(ODESpring.Tvec),'r--')
legend('ODE Solver','Analytical')
title('RK4')

%% Problem 6

%I'm using the object above with these changes:

ODESpring.dt=.05;
%this is the period
Period=2*pi/(sqrt(k/m))
ODESpring.T=2*Period

%Re-run K4
ODESpring.RK4();

%Position:
Yposition=ODESpring.X(1,:);
%Velocity:
Yvelocity=ODESpring.X(2,:);
%Force:
Yforce=-k*ODESpring.X(1,:);

%make figure
fH=figure;

%Setup movie: 
%If you are having trouble with 'mp4' try another compression method:
% help VideoWriter
filename='HW4_MassOnSpringMovie.mp4'
VidObj=VideoWriter(filename,'MPEG-4')

%Set FrameRate property so it runs at real time
VidObj.FrameRate=1/ODESpring.dt;

%Open video file for writing
VidObj.open();

for ii=1:length(ODESpring.Tvec)
    figure(fH) %make 'fH' the current figure
    %Plot the position of the mass
    plot(0,Yposition(ii),'ko','linewidth',5)
    axis([-1 1 -10 10])
    hold on %so that other plots will add to existing plot

    %Draw velocity vector
    quiver(0,Yposition(ii),0,Yvelocity(ii),0,'g','linewidth',4)
    
    %Draw Force vector
    quiver(0,Yposition(ii),0,Yforce(ii),0,'r','linewidth',3)
    
    %Draw y0
    plot([-3 3],[0 0],'k--')
    
    hold off %so new plot will won't add to this one
    
    xlabel('X(m)')
    ylabel('Y(m)')
    title('Mass on a spring')
    legend('Mass','Velocity (m/s)','Force (N)')
    axis equal
    %now add frame to movie
    currFrame = getframe();
    VidObj.writeVideo(currFrame);
    
end

%close video file
VidObj.close() 











