function P1()
% When is Venus Brightest?
close all
% Plot the brightness function and ascertain a good starting
% interval, keeping in mind that -B(t) must be convex on [L,R].
n= 400;
L = 0;
R = 200;
t = linspace(L,R,n);
for k=1:n
    bval(k) = B(t(k));
end
plot(t,bval);
shg
title('Looks like [20 60] would be a good intial interval...','fontsize',14)
set(gcf,'position',[100 100 1200 800])
pause(3)
L = 20;
R = 60;
tStar = Golden(@B,L,R);
V = Venus(tStar);
E = Earth(tStar);
thetaStar = acos((V-E)'*V/(norm(V)*norm(E-V)))*180/pi;

%%% Set up a simulation...
close all
figure
set(gcf,'position',[100 100 1200 800])
% Display the two orbits...
thetaVals = linspace(0,2*pi,100); c = cos(thetaVals); s = sin(thetaVals);
subplot(1,2,1)
plot(0,0,'.k','markersize',10)
hold on
plot(66.01*c,65.55*s,'m',94.5*c,91.5*s,'b');
axis equal square off
% Display the B(t)...
subplot(1,2,2)
tF = 400;
tvals = linspace(0,tF,tF+1);
Bvals = zeros(1,tF+1);
for k=1:tF+1
    Bvals(k) = B(tvals(k));
end
plot(tvals,Bvals)
axis square
hold on
    
% Simulate...
for k=1:tF
    V = Venus(tvals(k));
    E = Earth(tvals(k));
    subplot(1,2,1)
    h1 = plot(V(1),V(2),'om',E(1),E(2),'ob',[0 V(1) E(1)],[0 V(2) E(2)],'k');
    theta = acos((V-E)'*V/(norm(V)*norm(E-V)))*180/pi;
    title(sprintf('theta(t) = %6.2f',theta),'fontsize',16)
    subplot(1,2,2)
    h2 = plot(tvals(k),Bvals(k),'.r','Markersize',15);
    title(sprintf('t = %3d   B(t) = %10.3e',k,-Bvals(k)),'fontsize',16)
    pause(.05)
    delete(h1)
    delete(h2)
end
% Display the optimal position...
subplot(1,2,1)
V = Venus(tStar);
E = Earth(tStar);
h1 = plot(V(1),V(2),'om',E(1),E(2),'ob',[0 V(1) E(1)],[0 V(2) E(2)],'k');
title(sprintf('theta(Star) = %6.2f',thetaStar),'fontsize',16)
% Display the Optimal time and B-value
subplot(1,2,2)
title(sprintf('tStar = %10.6f  B(tStar) = %10.3e',tStar,-B(tStar)),'fontsize',16)
plot(tStar,B(tStar),'.r','Markersize',15);
hold off

  function z = B(t)
% negative of the brightness function
V = Venus(t);
E = Earth(t);
E2V = V - E;
d = norm(E2V);
cos_mu = V'*E2V/(d*norm(V));
F = (1+cos_mu)/2;
z = -F/d^2;

  function V = Venus(t);
% Location of Venus
theta = 2*pi*t/224.7;
x = 66.01*cos(theta);
y = 65.55*sin(theta);
V = [x;y];

function E= Earth(t);
% Location of Earth
theta = 2*pi*t/365.25;
x = 94.5*cos(theta);
y = 91.5*sin(theta);
E = [x;y];

% -1 for starting interval unnecessarily large (length ~100 or more)
% up to -2 for inadequate justification of starting interval choice
