Assignment 5 Scripts
P1
% Script P1Hints
close all
% >>>>>>>>>>>>> Show the Tangent Lines >>>>>>>>>>>>>>>
% Draw the Sun
R = 1; % The radius of the Sun
theta = linspace(0,2*pi);
plot(R*cos(theta),R*sin(theta))
axis([-2 10 -4 4])
axis equal
hold on
% Draw the inner planet and the tangent lines to the Sun
x = 4; y = 2; % Inner planet is at (x,y);
plot(x,y,'o')
rho = sqrt(x^2+y^2);
mu = R/rho;
% The left and right tangent points ...
Left = mu*[x*mu+y*sqrt(1-mu^2);y*mu-x*sqrt(1-mu^2)];
Right = mu*[x*mu-y*sqrt(1-mu^2);y*mu+x*sqrt(1-mu^2)];
% Connect...
plot([x Left(1)],[y Left(2)],'r')
plot([x Right(1)],[y Right(2)],'g')
hold off
title('The left (red) and right (green) tangent lines.')
% >>>>>>>>>>>>>>> Detecting the start and the finish of a transit >>>>>>>>>>
Event = cell(5,1);
Event{1} = 'Before Transit'; Event{2} = 'Start of Transit'; Event{3} = 'In Transit';
Event{4} = 'End of Transit'; Event{5} = 'After Transit';
k = 0;
for t = [0 1.38 2 2.6 4]
figure
R = 1; % The radius of the Sun
theta = linspace(0,2*pi);
plot(R*cos(theta),R*sin(theta))
axis([-2 10 -4 4])
axis equal
hold on
x = 4*cos(2*(t-2)*pi/15); y = 4*sin(2*(t-2)*pi/15); % Inner planet
u = 8*cos(2*(t-2)*pi/30); v = 8*sin(2*(t-2)*pi/30); % Outer planet
plot(x,y,'o',u,v,'*',[x u],[y v],'k')
% The left and right tangent points ...
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)];
Right = mu*[x*mu-y*sqrt(1-mu^2);y*mu+x*sqrt(1-mu^2)];
text(x,y-.4,'P'); text(u,v+.3,'Q'); text(Left(1),Left(2)-.3,'L');text(Right(1),Right(2)+.3,'R')
plot([x Left(1)],[y Left(2)],'r')
plot([x Right(1)],[y Right(2)],'g')
sRed = Sine(Left, [x;y], [u;v]);
sGreen = Sine(Right,[x;y], [u;v]);
k = k+1; text(6,2,Event{k});
hold off
title(sprintf('sin(angle QPL)= %4.2f sin( angle QPR) = %4.2f',sRed,sGreen))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% P1: Transits
% Orbit data for Mercury, Venus, Earth, Mars, Jupiter, and Saturn
Names = cell(6,1);
Name{1} = ' Mercury '; Name{2} = ' Venus '; Name{3} = ' Earth ';
Name{4} = ' Mars '; Name{5} = ' Jupiter '; Name{6} = ' Saturn ';
A = [69.8200; 108.9441; 152.1028; 249.2205; 816.5574; 1513.3940];
P = [46.0032; 107.4819; 147.1053; 206.6782; 740.6071; 1353.6339];
Y = [ 87.9668; 224.7018; 365.2467; 686.9878; 4332.6759; 10759.4359] ;
phi = [ 77.5158; 131.6175; 103.0033; 336.1308; 14.3933; 93.1322];
tau = [-40.2570; -86.8795; 61.9183; 607.7758; -1688.6236; -134.5013];
Planet = cell(6,1);
for k=1:6
Planet{k} = struct('Y',Y(k),'tau',tau(k),'A',A(k),'P',P(k),'phi',phi(k));
end
clc
rand('seed',0)
for i=1:5
for j=i+1:6
% Look for transits of Planet(i) as an observer on Planet{j}.
% Pick 5 random t values during the next several revolutions of Planet{i}.
t0 = 0;
disp(['Transits of' Name{i} 'as seen from' Name{j} ' ... '])
disp(' ')
disp(' Current Time till next Transit Transit Duration' )
disp(' Date Days Hours Minutes Days Hours Minutes')
disp('------------------------------------------------------------------------')
for k=1:5
t0 = floor(t0 + (1+rand(1,1))*Planet{i}.Y);
[tStart,tFinish] = NextTransit(Planet{i},Planet{j},t0);
% How long till the start of the transit?
T = tStart-t0;
D = floor(T); H = floor((T - D)*24); M = round( (T - D - H/24)*24*60);
% How long was the transit?
T1 = tFinish - tStart;
D1 = floor(T1); H1 = floor((T1 - D1)*24); M1 = round( (T1 - D1 - H1/24)*24*60);
disp(sprintf( '%10.0f %10d %5d %5d %10d %5d %5d',t0,D,H,M,D1,H1,M1))
end
disp(' ');disp('>>>>> ');disp(' ')
end
end