A3 Scripts and Functions
A3A
% Script A3A
% Illustrates maxSpline
close all
bigscreen
n = 17;
randn('seed',0);
x = sort(linspace(-2,2,n)' + .01*randn(n,1));
y = randn(n,1);
S = spline(x,y);
z = linspace(x(1),x(n),200);
sVals = ppval(S,z);
[zStar,maxS] = maxSpline(S);
plot(z,sVals,x,y,'o',zStar,ppval(S,zStar),'*r')
title(sprintf('xStar = %18.12f, S(xStar) = %18.12f',zStar,maxS))
ShowCrescent
% ShowCresent
% Takes a look at the Venus Cresent during a complete revolution
close all
bigscreen
n = 73;
theta = linspace(0,2*pi,n);
dE = 92.9570;
dV = 67.2322;
R = 3759;
for k=1:n
Crescent(theta(k),dV,dE,R)
pause(.3)
end
Crescent
function Crescent(theta,dV,dE,R)
% Shows what an inner planet V looks like from an outer planet E where theta is the
% E-Sun-V angle in radians. (Think: V = Venus, E = Earth)
% dV = radius of inner planet orbit (assumed circular)
% dE = radius of outer planet orbit (assumed circular)
% R = radius of the inner planet
% theta is the polar angle that defines the location of V
% in a fixed (non moving) coordinate system that puts the Sun
% at (0,0), and E at (dE,0):
subplot(1,2,1)
% In this window we "look down" on the E-Sun-V system.
% We use the fixed coordinate system
% Display the Sun and the connecting line segments
r = dE/12; % V's Radius in the display. (Not its radius as a planet)
E = [dE;0]; % E's position
V = dV*[cos(theta);sin(theta)]; % V's position
plot(0,0,'*r')
hold on
plot([0 E(1)],[0 E(2)],'--b')
plot([0 V(1)],[0 V(2)],'--m')
plot([E(1),V(1)],[E(2),V(2)],'--k')
% Display the E as blue dot and V (for now) as a black disk
kappa = linspace(0,2*pi); cVals = cos(kappa); sVals = sin(kappa);
fill(V(1)+r*cVals,V(2)+r*sVals,'k')
fill(E(1)+r*cVals/2,E(2)+r*sVals/2,'b')
text(90,-12,'E')
text(-3,-12,'S')
% Color the appropriate semicircle of Venus...
mu = theta+pi/2;
tau = linspace( mu, mu+pi);
u = V(1) + r*cos(tau);
v = V(2) + r*sin(tau);
fill([u u(1)],[v v(1)],'y')
axis equal
axis(1.1*[-dE dE -dE dE])
% Display theta = V-Sun-E angle and phi = V-E-Sun angle
plot(20*cos(linspace(0,theta)),20*sin(linspace(0,theta)),'m')
phi = atan(sin(theta)*dV/(dE - cos(theta)*dV));
plot(dE+20*cos(linspace(pi-phi,pi)),20*sin(linspace(pi-phi,pi)),'k')
title(sprintf('theta = %3.1f (deg) phi = %3.1f (deg)',theta*180/pi,phi*180/pi))
% Display the coordinates of S and E in the "moving" coordinate system
% In this system, V = (0,0) and E is at (0,-dEV) where dEV is the
% V-to-E distance. Since dEV varies with theta, it is a "moving" coordinate
% system.
dEV = (dE - dV*cos(theta))/cos(phi);
S_move = [-dE*sin(phi) ; -dEV+dE*cos(phi) ];
E_move = [ 0 ; -dEV];
xlabel(sprintf('S_{move} = [%5.2f,%5.2f] E_{move} = [ 0 , %5.2f]',S_move(1),S_move(2),E_move(2)))
hold off
subplot(1,2,2)
% In this window we display the Crescent.
% In the language of the moving coordinate system, the observer is at E,
% i.e., (0,-dEV,0).
% The observer is looking at V, a sphere of radius R
% centered at (0,0,0).
% The sphere is illuminated by S located at (S_move(1),S_move(2),0)
% The x-values of the right and left edge of the V-disk as a function of z ...
z = linspace(-R,R);
RightEdge = sqrt(R^2 - z.^2);
LeftEdge = -RightEdge;
% Display a full black disk...
fill([LeftEdge RightEdge(100:-1:1)],[z z(100:-1:1)],'k')
% Some trigonometry associated with the location of
% V is necessary to compute the terminator.
Q = E-V; Q = Q/sqrt(Q(1)^2 + Q(2)^2);
P = [sin(theta);-cos(theta)];
s = Q(1)*P(2) - Q(2)*P(1);
% Color in the crescent
hold on
if theta <=pi
% The crescent is in between the terminator and the left edge.
Term = s*RightEdge;
fill([ Term LeftEdge(100:-1:1)], [z z(100:-1:1)],'y')
else
% The crescent is in between the terminator and the right edge.
Term = s*LeftEdge;
fill([ Term RightEdge(100:-1:1)], [z z(100:-1:1)],'y')
end
axis equal off
axis([-R R -R R])
title('View of V from E')
%title(sprintf('Crescent Energy = %10.6f',CrescentEnergy(theta,dV,dE,R)))
hold off
shg
A3B
% A3B
% Answers the question, when does an inner planet appear brightest
% to an outer planet.
% The four planets nearest the Sun...
Planet(1) = struct('name','Mercury','d',35.9773*10^6,'R',1516);
Planet(2) = struct('name','Venus','d',67.2322*10^6,'R',3759);
Planet(3) = struct('name','Earth','d',92.9570*10^6,'R',3963);
Planet(4) = struct('name','Mars','d',141.6102*10^6 ,'R',2111);
close all
n = 37;
theta = linspace(0,pi,n)';
tvals = linspace(0,pi,500)';
for i=1:3
for j = i+1:4
% Inner planet = Planet(i)
% Outer planet = Planet(j)
% Compute the elongation angle and the energy for each theta
dV = Planet(i).d;
dE = Planet(j).d;
R = Planet(i).R;
Energy = zeros(n,1);
phi = zeros(n,1);
for k=1:n
phi(k) = atan(sin(theta(k))*dV/(dE - cos(theta(k))*dV));
Energy(k) = CrescentEnergy(theta(k),dV,dE,R);
end
figure
subplot(2,1,1)
% Fit a spline through the elogation values.
% Display the spline and its maximum value.
S_phi = spline(theta,phi);
[thetaStar,phiMax] = maxSpline(S_phi);
plot(tvals,ppval(S_phi,tvals),thetaStar,phiMax,'*r');
title(sprintf('thetaStar = %7.3f (deg) phiMax = %7.3f (deg)',thetaStar*180/pi,phiMax*180/pi))
ylabel('phi')
xlabel(['theta ' '(' Planet(i).name ' and ' Planet(j).name ')'])
subplot(2,1,2)
% Fit a spline through the energy values.
% Display the spline and its maximum value.
S_Energy = spline(theta,Energy);
[thetaStar,EnergyMax] = maxSpline(S_Energy);
plot(tvals,ppval(S_Energy,tvals),thetaStar,EnergyMax,'*r');
title(sprintf('thetaStar = %7.3f (deg) EnergyMax = %7.3e',thetaStar*180/pi,EnergyMax))
ylabel('Energy')
xlabel(['theta ' '(' Planet(i).name ' and ' Planet(j).name ')'])
shg
end
end