%Support function for essential_distMinAnglePair implementing Newton's search
function [tMin,fMin]=essential_distMinAnglePair_dfNewton(m1,p1,c1,m2,p2,c2,tMin,tLow,tHigh)
tolDist=1e-8;
for i=1:100
%     d=dfi(m1,p1,c1,tMin)+dfi(m2,p2,c2,tMin);
%     dd=ddfi(m1,p1,c1,tMin)+ddfi(m2,p2,c2,tMin);
    %The code below unrolls the following calls
    %     f1=fi(m1,p1,c1,tMin);
    %     f2=fi(m2,p2,c2,tMin);
    %     d=dfi2(m1,p1,f1,tMin)+dfi2(m2,p2,f2,tMin);
    %     dd=ddfi2(m1,p1,f1,tMin)+ddfi2(m2,p2,f2,tMin);
    mc1=m1*cos(tMin+p1);
    mc2=m2*cos(tMin+p2);
    f1=acos(clip((m1*sin(tMin+p1)+c1-1)/2));
    f2=acos(clip((m2*sin(tMin+p2)+c2-1)/2));
    sf1=2*sin(f1);
    sf2=2*sin(f2);
    d1=-f1*mc1/sf1;
    d2=-f2*mc2/sf2;
    d=d1+d2;
    eztuSq1=(mc1/sf1)^2;
    dd1=eztuSq1+f1/2*cot(f1/2)*(1-eztuSq1);
    eztuSq2=(mc2/sf2)^2;
    dd2=eztuSq2+f2/2*cot(f2/2)*(1-eztuSq2);
    dd=dd1+dd2;
        
            
    tOld=tMin;
    tMin=max(tLow+tolDist,min(tHigh-tolDist,tOld-d/dd));
    if abs(tMin-tOld)<tolDist
        break
    end
end
fMin=f1^2+f2^2;

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

