%HW5_Solutions
%% Problem #1

% declaring symbolic variables
syms x(t) C W;
Dx = diff(x); %first derivative of x
D2x = diff(x,2); %second derivative of x

%General solution
SoluGen = dsolve(D2x + 2*C*W*Dx + (W^2)*x == 0 )
simplify(SoluGen,'IgnoreAnalyticConstraints', true,'Steps',50)
pretty(SoluGen)

%suppose that the mass is streached out by x(0) = 5 and it is held steady
%so Dx(0)/Dt = 0

Solu = dsolve(D2x + 2*C*W*Dx + (W^2)*x == 0, x(0) == 5, Dx(0) == 0)
simplify(Solu,'IgnoreAnalyticConstraints', true,'Steps',50)
pretty(Solu)


% This system is underdamped if b^2 < 4*a*c
% for our problem the underdamped case is given by (2*eps*omeg)^2 < 4*omeg^2
% which gives C < 1. I take C = 0.1
% W = sqrt(k/m) = 1

Xsol = subs(Solu,{W,C},{1,0.1})

% using ezplot to plot the symbolic function x(t)
figure
ezplot(Xsol,[0,20])
title('X(t)')
xlabel('t')
ylabel('x')

%% Problem #2
%The class TK4 utilized 

%Planet motion
Orb = RK4();
Orb.Dt = 0.01;
Orb.T = 10; %Total Time
Orb.X = zeros(4,length(0:Orb.Dt:Orb.T));
Orb.X0 = [10;0;0;5]; %initial x-position and velocities
M = 1; %mass

%Do this is two steps...
Fx = @(x,y) -500*x/(x^2+y^2)^(3/2); %x-force
Fy = @(x,y) -500*y/(x^2+y^2)^(3/2); %y-force
Orb.f = @(X)[X(3);X(4);Fx(X(1),X(2))/M;Fy(X(1),X(2))/M ];


[Orb] = Orb.calc(); %calling Calc-function inside the class
Orb.plotFig1() %calling PlotFig-function inside the class

%Spring
Spr = RK4();
Spr.Dt = 0.1;
Spr.X0 = [5;0]; %initial x-position and velocity
Spr.T = 20; %Total Time
Spr.X = zeros(2,length((0:Spr.Dt:Spr.T)));
C = 0.1;
W = 1;
k=1;
m=1;
Spr.f = @(X)[X(2); -2*C*W*X(2)-W^2*X(1)];
[Spr] = Spr.calc(); %calling Calc-function inside the class
Spr.plotFig2 %calling PlotFig-function inside the class




%% Problem #3

X=Spr.X(1,:)
V=Spr.X(2,:)
T=Spr.Tvec;
NPoints=length(T)
DT=T(2)-T(1);

C = 0.1;
W = 1;
k=1;
m=1;

F=m*(-2*C*W*V-W^2*X);

%plotting the figure
figure
axis([-1 1 min(X)*2  2*max(X)]);
xlabel('X')
ylabel('T')
for ii = 1:NPoints
plot(0,X(ii),'rO','LineWidth',2)
hold on
quiver(0,X(ii),0,V(ii),0) %adding the velocity vector to the plot
quiver(0,X(ii),0,F(ii),0) %adding the Force vector to the plot
xlabel('X')
ylabel('Y')
axis([-1 1 min(X)*2  2*max(X)]);
legend('Mass Location','Velocity','Force')
pause(DT) %after plotting each point there is a 0.05 s pause.
hold off
end
    
    
    





