function FindConj()
% Uses fzero to estimates the spacing between Mercury-Earth Conjunctions.
% This happens when Mercury, the Sun, and Earth are co-linear with Mercury
% "in the middle."
% We want to estimate the time in between conjunctions.
% The function SineMercEarth(t) returns the sine of angle Earth-Sun-Mercury
% at time t so this function is zero at a conjunction. But it is also
% zero when the Sun is in between Mercury and Earth.
% To gain intuition, we plot this function

clc
t = linspace(0,600,500);
plot(t,SineMercEarth(t))
xlabel('t','fontsize',16)
grid on
title('Sine of the Earth-Sun-Mercury Angle','fontsize',18)
shg
% We note two things from the plot. (a) t = 0 is a conjunction and (b)
% the spacing between conjunctions is about 115 days.
% So we use fzero to go after the precise times...
GapEst = 115;
disp(sprintf('Next ten conjunctions, Spacing Estimate = %8.3f',GapEst))
disp(' ')
t = zeros(11,1);
disp('Conjunction Time Spacing')
disp('-----------------------------------')
for k=1:10;
t(k+1) = fzero(@SineMercEarth,k*GapEst);
disp(sprintf(' %2.0f %8.3f %8.3f',k,t(k+1),t(k+1)-t(k)))
end

function s = SineMercEarth(t)
% The sine of the Mercury-Sun-Earth angle at time t

% Mercury location:
xm = -11.9084 + 57.9117*cos(2*pi*t/87.97);
ym = 56.6741*sin(2*pi*t/87.97);

% Earth location:
xe = -2.4987 + 149.6041*cos(2*pi*t/365.25);
ye = 149.5832*sin(2*pi*t/365.25);

s = (xm.*ye - xe.*ym)./(sqrt(xm.^2 + ym.^2).*sqrt(xe.^2 + ye.^2));