function [tMin,fMin,tBreak1,tBreak2,tMinAll]=essential_distMinAnglePair_base(Q211,Q212)
flagCheckFirstDer=true;
flagUseNewton=true;     %Note: requires flagCheckFirstDer=true
tolMZero=1e-15;
tMinAll=[];

[tBreak1,~,~,c1,m1,p1]=essential_distMinAnglePair_discontinuityDistance(Q211);
[tBreak2,~,~,c2,m2,p2]=essential_distMinAnglePair_discontinuityDistance(Q212);

%check for the degenerate case where the cost is constant
if abs(m1)<tolMZero && abs(m2)<tolMZero
    tMin=0;
    fMin=2*pi^2;
    tMinAll=0;
else
    %ft=@(t)  acos((m1*sin(t+p1)+c1-1)/2)^2+acos((m2*sin(t+p2)+c2-1)/2)^2;

    if abs(modAngle(tBreak1-tBreak2))<1e-8
        tMin=tBreak1+pi;
        fMin=0;
%         theta1=@(t) acos((m1*sin(t+p1)+c1-1)/2);
%         theta2=@(t) acos((m2*sin(t+p2)+c2-1)/2);
% 
%         ft=@(t) 0.5*(theta1(t)^2+theta2(t)^2);
%         [tMin,fMin]=fminbnd(ft,tBreak1,tBreak1+2*pi);
    else
        tSearch1=tBreak1;
        tSearch2=tBreak2;
        if tSearch1>tSearch2
            tSearch1=tSearch1-2*pi;
        end

        if flagCheckFirstDer
            %compute derivatives of each term at discontinuity points
            df1Break1=essential_distMinAnglePair_computeDfBreak(tBreak1,Q211);
            df2Break2=essential_distMinAnglePair_computeDfBreak(tBreak2,Q212);
%             disp('[df1Break1 df2Break2]')
%             disp([df1Break1 df2Break2])
            %compute derivative of each term at other's discontinuity
            %(unroll two calls to dfi)
            theta1Break2=acos(clip((m1*sin(tBreak2+p1)+c1-1)/2));
            df1Break2=-theta1Break2*(m1*cos(tBreak2+p1))/(2*sin(theta1Break2));
            theta2Break1=acos(clip((m2*sin(tBreak1+p2)+c2-1)/2));
            df2Break1=-theta2Break1*(m2*cos(tBreak1+p2))/(2*sin(theta2Break1));

            %compute left and right derivatives of sum of the two terms
            dfBreak1n=+df1Break1+df2Break1;
            dfBreak1p=-df1Break1+df2Break1;
            dfBreak2n=+df2Break2+df1Break2;
            dfBreak2p=-df2Break2+df1Break2;

            flagSearch1=false;
        %     plot([tBreak1 tBreak2],[dfBreak1p dfBreak2p],'cx','MarkerSize',10)
        %     plot([tBreak1 tBreak2],[dfBreak1n dfBreak2n],'mx','MarkerSize',10)
            if sign(dfBreak1p)~=sign(dfBreak2n)
                if flagUseNewton
                    %parabolic prediction of min
                    tMin0=tSearch1-dfBreak1p*(tSearch2-tSearch1)/(dfBreak2n-dfBreak1p);
                    %tMin0=(tSearch1+tSearch2)/2;
                    [tMin,fMin]=essential_distMinAnglePair_dfNewton(m1,p1,c1,m2,p2,c2,tMin0,tSearch1,tSearch2);
                    %fMin=essential_distMinAnglePair_ft(m1,p1,c1,m2,p2,c2,tMin);
                else
                    [tMin,fMin]=fminbnd(essential_distMinAnglePair_ft,tSearch1,tSearch2);
                end
                tMinAll=[tMinAll tMin];
                flagSearch1=true;
            end
            tSearch1=tSearch1+2*pi;
            if sign(dfBreak2p)~=sign(dfBreak1n)
                if flagUseNewton
                    %parabolic prediction of min
                    tMin0=tSearch2-dfBreak2p*(tSearch1-tSearch2)/(dfBreak1n-dfBreak2p);
                    %tMin0=(tSearch1+tSearch2)/2;
                    [tMin2,fMin2]=essential_distMinAnglePair_dfNewton(m1,p1,c1,m2,p2,c2,tMin0,tSearch2,tSearch1);
                    %fMin2=essential_distMinAnglePair_ft(m1,p1,c1,m2,p2,c2,tMin2);
                else
                    [tMin2,fMin2]=fminbnd(essential_distMinAnglePair_ft,tSearch2,tSearch1);
                end
                if ~flagSearch1 || (flagSearch1 && fMin2<fMin)
                    tMin=tMin2;
                    fMin=fMin2;
                end
                tMinAll=[tMinAll tMin2];
            end
        else
            [tMin1,fMin1]=fminbnd(essential_distMinAnglePair_ft,tSearch1,tSearch2);
            tSearch1=tSearch1+2*pi;
            [tMin2,fMin2]=fminbnd(essential_distMinAnglePair_ft,tSearch2,tSearch1);
            if fMin1<fMin2
                tMin=tMin1;
                fMin=fMin1;
            else
                tMin=tMin2;
                fMin=fMin2;
            end
        end
    end
end

function v=clip(v)
v=min(1,max(-1,v));


% function f=fi(m,p,c,t)
% f=acos((m*sin(t+p)+c-1)/2);
% 
% function d=dfi2(m,p,theta,t)
% dtheta= -(m*cos(t+p))/(2*sin(theta));
% d=theta*dtheta;
% 
% function dd=ddfi2(m,p,theta,t)
% eztuSq=(m*cos(t+p)/(2*sin(theta)))^2;
% dd=eztuSq+theta/2*cot(theta/2)*(1-eztuSq);
% 
% function d=dfi(m,p,c,t)
% theta=acos((m*sin(t+p)+c-1)/2);
% dtheta= -(m*cos(t+p))/(2*sin(theta));
% d=theta*dtheta;
% 
% function dd=ddfi(m,p,c,t)
% theta=acos((m*sin(t+p)+c-1)/2);
% eztuSq=(m*cos(t+p)/(2*sin(theta)))^2;
% dd=eztuSq+theta/2*cot(theta/2)*(1-eztuSq);



