%% Chapter 4: Kinematics
%% Kepler's Laws of Planetary Motion
% In the first part of this M-file, we produce a picture illustrating two of Kepler's laws. First, planets follow elliptical orbits. We start with the equation of a particular ellipse.

syms x y
ellipse = (x^2/13^2 + y^2/12^2 == 1);

%
% Now we'll move one focus to the origin and change to polar coordinates.

newellipse = subs(ellipse, x, x+5);
syms r theta
newellipse = simplify(subs(newellipse,[x,y],[r*cos(theta),r*sin(theta)]))
rsols = solve(newellipse, r)

%
% Since we want r to be positive, we'll use the first solution.

polarellipse = simplify(rsols(1))

%
% Now we compute area of a sector in polar coordinates.

polararea = int(polarellipse^2/2)

%
% We plan to plot two shaded areas: one between angles 5*pi/6 and 11*pi/12, and the second from 0 to an angle alpha that encloses the same area as the first sector. Here is how we compute alpha.

sector = subs(polararea, 11*pi/12) - subs(polararea, 5*pi/6)
f = @(alpha) double(subs(polararea, theta, alpha) - sector);
endvalue =  fzero(f, 0.1)

%
% Now we plot the polar region bounded by newellipse from 11*pi/12 to 5*pi/6 and also the polar region from 0 to endvalue. We need a vectorized function for plotting and vectors of values of theta.

ellfun = @(theta) 144./(5*cos(theta) + 13);
t0 = linspace(0, 2*pi, 100);
t1 = linspace(0, endvalue, 20);
t2 = linspace(5*pi/6, 11*pi/12, 20);
polarplot(t0, ellfun(t0), 'k'), hold on
for j = 1:20
  polarplot([t1(j),t1(j)], [0, ellfun(t1(j))], 'r')
  polarplot([t2(j),t2(j)], [0, ellfun(t2(j))], 'r')
end
hold off

%
%% Studying Equations of Motion with MATLAB
% We solve the equation of motion for a charged particle with charge q, mass m, and initial velocity v(0)=(v1, v2, v3) in a uniform magnetic field B = B*k*. First we compute the right-hand side of the equation.

syms q c m t v1 v2 v3 B
syms x(t) y(t) z(t)
Dx=diff(x); Dy=diff(y); Dz=diff(z);
rhs = (-q*B/(c*m))*cross([0,0,1],[Dx, Dy, Dz])
sol = dsolve([diff(x,2),diff(y,2),diff(z,2)] == rhs, ...
    [x(0),y(0),z(0)]==[0,0,0], [Dx(0),Dy(0),Dz(0)]==[v1,v2,v3]);

%
% Since there is no force in the z-direction, we expect constant-speed motion in the z-direction. Indeed:

sol.z

%
% Now let's simplify further to see what this looks like. We're take v1=0, v2=1, and q=c=m=B=1.

assume(t, 'real'); assume(v3, 'real')
simplesol = ...
    simplify(subs([sol.x,sol.y,sol.z],[c,m,q,B,v1,v2],[1,1,1,1,0,1]))

%
% We check that this is real and write it as follows in terms of trig functions.

simplify(imag(simplesol))
simplify(real(simplesol))

%
% Now we redo the same example using ode45 and plot the result. Here w is an array with six columns, corresponding to the 6 coordinates in phase space (x, y, z, and the corresponding components of the velocity). Note that the initial condition becomes w0 = [0; 0; 0; 0; 1; 1], and w' = [w(4); w(5); w(6); w(5); -w(4); 0].

[t,w] = ode45(@(t,w) [w(4); w(5); w(6); w(5); -w(4); 0],...
    [0 6*pi],[0; 0; 0; 0; 1; 1]);
plot3(w(:,1),w(:,2),w(:,3))

% Note the characteristic shape of a helix.