  function [StartTime,Duration] = Transits(i,j,T)
% StartTime and Duration that specify the start times (days)  and durations
% (in minutes) of all the transits of planet i as seen from planet j.
% We assumes 1 <= i < j <= 4 and that 1=Mercury, 2=Venus, 3=Earth, and
% 4=Mars.

% To get approximate transit start times, we use the information about the
% planet's initial angle and its period. Here is the information for
% the four planets obtained from Location.m...
 phi = [ 77.45645; 131.53298; 102.94719; 336.04084];
 Y   = [  87.968; 224.695; 365.242; 686.930];
% For planet i...
 phi_i = phi(i)*pi/180;
 Y_i = Y(i);
 % For planet j ...
 phi_j = phi(j)*pi/180;
 Y_j = Y(j);
 
 % Assume circular orbits, a very reasonable assumption. Let's reason about 
 % conjunction times given this assumption. Planet i will
 % be on the Sun-to-Planet j line segment if
 %
 %  phi_i + 2*pi*(t/Y_i) = phi_j + 2*pi*(t/Y_j) + 2*pi*k for integers k.
 % 
 % Solving this for t we see that t = t0 + k*del where
 
 t0 = ((phi_j-phi_i)/(2*pi))/(1/Y_i - 1/Y_j);
 del = 1/(1/Y_i - 1/Y_j);
 
 % Notice that del is positive since Y_j > Y_i.
 
 % Here is a vector of approximate conjunction times around
 % the interval [0,T]...
 
 ApproxConj = t0:del:T+del;
 
 % Now let us visit each of these approximate conjunction times and examine the
 % transit that is in that vicinity...
 
 StartTime = [];
 Duration = [];
 % We'll use Approximate Conjunction time plus/minus tau for bracketing
 % interval.
 tau = 20;
 for k=1:length(ApproxConj)
     % Compiute the start time...
     tStart = fzero(@(t) Left(t,i,j),[ApproxConj(k)-tau,ApproxConj(k)+tau]);
     if 0<=tStart && tStart<=T
         % The start time is in [0,T] so compute the end time...
         tEnd = fzero(@(t) Right(t,i,j),[ApproxConj(k)-tau,ApproxConj(k)+tau]);
         StartTime = [StartTime;tStart];
         Duration  = [Duration;tEnd-tStart];
     end
 end

   function sigma = Left(t,i,j)
 % sigma is the sine of angle Left-Outer-Inner where at time t Outer is
 % the location of planet j, Inner is the location of planet i (i<j), and 
 % Left is the Left tangent point on the Sun.
 r = .6955;
 Inner = Location(t,i);
 Outer = Location(t,j);
 [Left,Right] = TP(Outer,r);
 % Unit vector in direction  Left-Outer...
 u = [Left.x-Outer.x ; Left.y-Outer.y]; u = u/norm(u);
 % Unit vector in direction  Inner-Outer...
 v = [Inner.x-Outer.x ; Inner.y-Outer.y]; v = v/norm(v);
 % The sine of the angle between u and v...
 sigma = u(1)*v(2) - u(2)*v(1);
 
   function sigma = Right(t,i,j)
 % sigma is the sine of angle Right-Outer-Inner where at time t Outer is
 % the location of planet j, Inner is the location of planet i (i<j), and 
 % Right is the Right tangent point on the Sun.
 r = .6955;
 Inner = Location(t,i);
 Outer = Location(t,j);
 [Left,Right] = TP(Outer,r);
 % Unit vector in direction  Left-Outer...
 u = [Right.x-Outer.x ; Right.y-Outer.y]; u = u/norm(u);
 % Unit vector in direction  Inner-Outer...
 v = [Inner.x-Outer.x ; Inner.y-Outer.y]; v = v/norm(v);
 % The sine of the angle between u and v...
 sigma = u(1)*v(2) - u(2)*v(1);
 % method not working -4
 % method not complete:
 % left-right confusion, negative duration -2
 % only works for 1st transit -3
 % initial time beyond range, wrong radius -1
 % slow but working: full credit
 % P3solve thought that output Duration of transits is measured by days
 % while the specification is minutes. 
