A5 Solutions
P1
function [tStart,tFinish] = NextTransit(Planet1,Planet2,t0)
% Planet1 and Planet2 are structures that specify the orbits of two planets, e.g.,
% Planet1.Y, Planet1.tau, Planet1.A, Planet1.P, Planet1.phi
% Planet2.Y, Planet2.tau, Planet2.A, Planet2.P, Planet2.phi
% Y and tau are in days, A and P are in millions of kilometers, and phi is in degrees.
% Assume that the orbit of Planet1 is inside the orbit of Planet2 and that
% Planet1.Y < Planet2.Y
% t0 is the current time (in days).
%
% The next transit of Planet 1 begins at tStart and ends at tFinish.
% The absolute error in tStart and tFinish should be less than one minute.
% Some critical quantities for estimating the time of the next transit.
tau1 = Planet1.tau; Y1 = Planet1.Y; phi1 = Planet1.phi;
tau2 = Planet2.tau; Y2 = Planet2.Y; phi2 = Planet2.phi;
% To get a reasonable estimate, assume the orbits are circular.
kappa = 1/((1/Y1)-(1/Y2)); % Periodicity of transits if we assume circular orbits
% At time t0 the polar angles (in radians) of the two planets are
%
% Planet1: theta1 = 2*pi*phi1/360 + 2*pi*(t0-tau1)/Y1
% Planet2: theta2 = 2*pi*phi2/360 + 2*pi*(t0-tau2)/Y2;
%
% After time delta_T the polar angles of the two planets are
%
% Planet1: psi1 = theta1 + 2*pi*delta_T/Y1
% Planet2: psi2 = theta2 + 2*pi*detla_T/Y2
%
% What is the smallest positive delta_T so that (psi1-psi2) is a whole number multiple
% of 2*pi? Simplifying the equation psi1 = psi2 +2*pi*k (k initeger) we get
%
% z + delta_T*kappa = k
% where
z = ( phi1/360 + (t0-tau1)/Y1 ) - (phi2/360 + (t0-tau2)/Y2 );
% Thus, by adding delta_T*kappa to k we are to obtain an integer. The
% smallest value that does this is 1-f where
f = z-floor(z);
% Thus,
delta_T = (1-f)*kappa;
% Find the start time of the next transit
tNext = t0 + delta_T;
del = kappa/10;
Interval = [tNext-del,tNext+del];
% Look for next transit in [tNext-del,tNext+del]
tol = .000000001; % Doesn't have to be this small
tStart = fzero(@TransitStart,Interval,optimset('TolX',tol),Planet1,Planet2);
if tStart<t0
% This can happen if a transit is going on at t0 or it t0 is just after
% a transit has occurred. In this case, look for the next transit in the
% vicinity of tNext+kappa
Interval = Interval + kappa;
tStart = fzero(@TransitStart,Interval,optimset('TolX',tol),Planet1,Planet2);
end
% Transits are never more than 2 days long
tFinish = fzero(@TransitFinish,[tStart,tStart+2],optimset('TolX',tol),Planet1,Planet2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = TransitStart(t,Planet1,Planet2)
% Planet1 specifies the orbit parameters Y, tau, A, P, and phi of the inner planet.
% Planet2 specifies the orbit parameters Y, tau, A, P, and phi of the outer planet.
% s is the sine of the angle defined by the left tangent point,
% Planet1, and Planet2 at time t.
% This function is zero at the start of the transit.
Inner = Locate(t,Planet1); % Inner planet location
Outer = Locate(t,Planet2); % Outer planet location
% Find the left tangent point
R = .69597; % Sun radius
x = Inner(1);
y = Inner(2);
rho = sqrt(x^2+y^2);
mu = R/rho;
Left = mu*[x*mu+y*sqrt(1-mu^2);y*mu-x*sqrt(1-mu^2)];
% Return the sine of the angle defined by points Left, Inner, and Outer
s = Sine(Left,Inner,Outer);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = TransitFinish(t,Planet1,Planet2)
% Planet1 specifies the orbit parameters Y, tau, A, P, and phi of the inner planet.
% Planet2 specifies the orbit parameters Y, tau, A, P, and phi of the outer planet.
% s is the sine of the angle defined by the right tangent point,
% Planet1, and Planet2 at time t.
% This function is zero at the end of the transit.
Inner = Locate(t,Planet1); % Inner planet location
Outer = Locate(t,Planet2); % Outer planet location
% Find the left tangent point
R = .69597; % Sun radius
x = Inner(1);
y = Inner(2);
rho = sqrt(x^2+y^2);
mu = R/rho;
Right = mu*[x*mu-y*sqrt(1-mu^2);y*mu+x*sqrt(1-mu^2)];
% Return the sine of the angle defined by points Right, Inner, and Outer
s = Sine(Right,Inner,Outer);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = Sine(A,B,C)
% A, B, and C are column 2-vectors
% Let theta be the angle formed as C-B is rotated counterclockwise
% to align with to A-B.
% s = sin(theta)
u = C-B; u = u/norm(u);
v = A-B; v = v/norm(v);
s = -v(1)*u(2)+v(2)*u(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function P = Locate(t,Planet)
% Planet specifies an orbit with parameters Y, tau, A, P, and phi
% P is a column 2-vector that specifies its location at time t
theta = 2*pi*(t - Planet.tau)/Planet.Y;
x = (Planet.P - Planet.A)/2 + ((Planet.P + Planet.A)/2)*cos(theta);
y = sqrt(Planet.A*Planet.P)*sin(theta);
c = cos(pi*Planet.phi/180);
s = sin(pi*Planet.phi/180);
P = [c -s ; s c]*[x ; y];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Transits of Mercury as seen from Venus ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
107 44 10 8 0 6 38
199 104 4 40 0 6 52
346 99 21 16 0 5 42
493 92 23 58 0 6 48
663 75 2 14 0 6 45
>>>>>
Transits of Mercury as seen from Earth ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
121 43 4 53 0 6 31
254 32 11 9 0 6 39
415 96 8 10 0 6 39
506 5 8 10 0 6 39
598 35 14 30 0 6 31
>>>>>
Transits of Mercury as seen from Mars ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
134 5 1 50 0 7 23
281 55 8 50 0 7 7
369 69 20 40 0 6 48
490 53 18 40 0 6 46
583 63 20 23 0 7 4
>>>>>
Transits of Mercury as seen from Jupiter ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
124 30 22 2 0 7 50
272 62 8 43 0 7 49
411 13 3 40 0 7 48
580 23 20 59 0 7 47
742 41 18 56 0 7 45
>>>>>
Transits of Mercury as seen from Saturn ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
134 6 13 34 0 8 2
230 87 18 25 0 8 1
375 31 8 56 0 8 1
499 84 14 10 0 8 1
648 24 4 55 0 8 1
>>>>>
Transits of Venus as seen from Earth ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
429 367 22 21 0 7 51
824 555 15 19 0 7 59
1107 272 15 19 0 7 59
1342 37 15 19 0 7 59
1732 233 11 54 0 7 56
>>>>>
Transits of Venus as seen from Mars ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
298 129 13 15 0 8 37
664 108 1 12 0 8 35
1058 37 13 28 0 8 41
1505 258 17 9 0 8 44
1811 296 10 17 0 8 27
>>>>>
Transits of Venus as seen from Jupiter ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
280 117 8 4 0 10 5
725 145 2 22 0 10 4
1112 231 9 14 0 10 2
1505 75 8 12 0 10 1
1876 178 22 22 0 9 59
>>>>>
Transits of Venus as seen from Saturn ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
241 108 1 9 0 10 24
607 201 11 29 0 10 24
1030 8 3 55 0 10 24
1315 182 10 57 0 10 24
1637 90 1 24 0 10 25
>>>>>
Transits of Earth as seen from Mars ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
645 94 15 17 0 9 18
1184 333 4 38 0 9 52
1636 653 2 20 0 10 7
2101 188 2 20 0 10 7
2597 462 20 59 0 10 2
>>>>>
Transits of Earth as seen from Jupiter ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
426 96 6 37 0 11 32
968 350 0 2 0 11 27
1661 56 7 53 0 11 24
2358 160 8 45 0 11 22
2745 173 22 46 0 11 23
>>>>>
Transits of Earth as seen from Saturn ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
695 120 18 39 0 12 0
1244 328 14 58 0 12 1
1797 153 21 31 0 12 1
2278 51 1 21 0 12 2
3003 82 0 47 0 12 3
>>>>>
Transits of Mars as seen from Jupiter ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
1026 108 0 58 0 13 39
1895 47 10 0 0 13 44
2644 108 4 30 0 13 36
3982 415 12 2 0 13 23
4719 498 9 13 0 13 35
>>>>>
Transits of Mars as seen from Saturn ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
1030 605 22 27 0 14 14
1980 395 11 4 0 14 18
2857 255 17 1 0 14 24
4171 408 0 44 0 14 35
5221 87 18 13 0 14 38
>>>>>
Transits of Jupiter as seen from Saturn ...
Current Time till next Transit Transit Duration
Date Days Hours Minutes Days Hours Minutes
------------------------------------------------------------------------
6344 7072 8 1 0 22 7
14753 5784 9 55 0 22 27
19302 1235 9 55 0 22 27
26934 942 20 41 0 23 22
34603 580 23 3 0 22 13
P2
% P2
% Some initializations
close all
i = sqrt(-1);
theta = linspace(0,2*pi,1000);
% Build intuition about g_r ...
n = 50; r = (1:n)/100; g_r = zeros(1,n);
for k=1:n
z = 1 + r(k)*(cos(theta) + i*sin(theta));
g = (((3*z -1).*z - 1).*z - 1)./(4*z.^3);
g_r(k) = max(abs(g));
end
plot(r,g_r)
ylabel('g_r')
xlabel('r')
title('g_r = max | (3z^3 - z^2 -z -1 )/(4z^3) | where z = 1 + r*exp(i*theta)')
grid
clc
disp(' rho r(rho)')
disp('--------------------')
for rho=(1:9)/10;
% Find a zero of f(r) = g_r - rho
L = r(1) ; fL = g_r(1) - rho; % f(1/n) < 0
R = r(n) ; fR = g_r(n) - rho; % f(.5) > 0
% Apply bisection...
while R-L > .0001
rm = (L+R)/2;
% Evaluate f at the midpoint of [L,R]
z = 1 + rm*(cos(theta) + i*sin(theta));
fm = max(abs( (((3*z -1).*z - 1).*z - 1)./(4*z.^3) )) - rho;
if fL*fm<=0
R = rm; fR = fm;
else
L = rm; fL = fm;
end
end
r_of_rho = (L+R)/2;
disp(sprintf(' %5.1f %8.4f',rho,r_of_rho))
end
%%%%%%%%%%%%%%%%%
rho r(rho)
--------------------
0.1 0.0600
0.2 0.1095
0.3 0.1511
0.4 0.1867
0.5 0.2176
0.6 0.2449
0.7 0.2690
0.8 0.2907
0.9 0.3103
P3
% P3
xstar = [0;1];
clc
% Newton's Method
disp('Newton''s Method')
disp(' k x(1) x(2) Error')
disp('-----------------------------------------------------------------')
k = 0;
x = [-.5 ; 1.4];
while (norm(x-xstar,inf) > eps)
F = [ (x(1)+3)*(x(2)^3 - 7)+18 ; sin(x(2)*exp(x(1))-1) ];
s = -(Jac(x)\F);
x = x + s;
k = k+1;
disp(sprintf('%5.0f %20.16f %20.16f %8.3e',k,x(1),x(2),norm(x-xstar,inf)))
end
% Broyden's Method
disp(' ')
disp('Broyden''s Method')
disp(' k x(1) x(2) Error')
disp('-----------------------------------------------------------------')
k = 0;
x = [-.5 ; 1.4];
J = Jac(x);
F = [ (x(1)+3)*(x(2)^3 - 7)+18 ; sin(x(2)*exp(x(1))-1) ];
while (norm(x-xstar,inf) > eps)
s = -(J\F);
x = x + s;
FNew = [ (x(1)+3)*(x(2)^3 - 7)+18 ; sin(x(2)*exp(x(1))-1) ];
J = J + (FNew - F - J*s)*(s/(s'*s))';
F = FNew;
k = k+1;
disp(sprintf('%5.0f %20.16f %20.16f %8.3e',k,x(1),x(2),norm(x-xstar,inf)))
end
%%%%%%%%%%%%%%%%%%%%%
Newton's Method
k x(1) x(2) Error
-----------------------------------------------------------------
1 -0.0553151357177096 1.0280665838357432 5.532e-002
2 -0.0001403508964315 1.0001574043269821 1.574e-004
3 -0.0000000177907770 1.0000000055513854 1.779e-008
4 0.0000000000000001 1.0000000000000000 1.181e-016
Broyden's Method
k x(1) x(2) Error
-----------------------------------------------------------------
1 -0.0553151357177096 1.0280665838357432 5.532e-002
2 0.0005099530701465 1.0001236434918623 5.100e-004
3 -0.0002338478636407 1.0000765608693296 2.338e-004
4 -0.0000408262210374 1.0000135978907720 4.083e-005
5 -0.0000001327508960 1.0000000453383613 1.328e-007
6 -0.0000000005391206 1.0000000001806633 5.391e-010
7 0.0000000000016679 0.9999999999994428 1.668e-012
8 -0.0000000000000007 1.0000000000000004 7.435e-016
9 0.0000000000000000 0.9999999999999999 1.110e-016