function [muo1, mu0]=occultquadVec(p,z,u1,u2,ToPlot,varargin)
% Aviv Comments
% ==================================================
% v.0.32 - vectorize rc,rj,rf
%        - Allow correcting behaviour near z==p, z=1-p
%        - Added optinal parameter TolFactor for tolerances control
%        - Removed all "CloseToOne" and added ToSmooth Parameter

%%%%%%%%%%%%%%%% - Mandel Agol comments
%
%   This routine computes the lightcurve for occultation
%   of a quadratically limb-darkened source without microlensing.
%   Please cite Mandel & Agol (2002) if you make use of this routine
%   in your research.  Please report errors or bugs to agol@tapir.caltech.edu
%       implicit none
%       integer i,nz
%       double precision z0(nz),u1,u2,p,muo1(nz),mu0(nz),
%      &       mu(nz),lambdad(nz),etad(nz),lambdae(nz),lam,
%      &       pi,x1,x2,x3,z,omega,kap0,kap1,q,Kk,Ek,Pk,n,ellec,ellk,rj
% 
%  Input:
% 
%  rs   radius of the source (set to unity)
%  z    impact parameter in units of rs
%  p    occulting star radius in units of rs
%  u1   linear    limb-darkening coefficient (gamma_1 in paper)
%  u2   quadratic limb-darkening coefficient (gamma_2 in paper)
% 
%  Optional input:
%
%  ToPlot:   Graphical output [default=false];
%  Tols:     A factor controling the three "ERRTOL" parameters for the three helper functions. If given, all three default ERRTOLs wil be scaled by it: [Default =1 ]
%  ToSmooth: if true, smooths the original Mandel-Agol model around z=p and around abs(z-p)=1 to produce better-behaving derivative. The value themselves are virtually unchanged. [default=true]
%
%  Output:
% 
%  muo1 fraction of flux at each z for a limb-darkened source
%  mu0  fraction of flux at each z for a uniform source
% 
%  Limb darkening has the form:
%   I(r)=[1-u1*(1-sqrt(1-(r/rs)^2))-u2*(1-sqrt(1-(r/rs)^2))^2]/(1-u1/3-u2/6)/pi
%  
%  To use this routine
% 
%  Now, compute pure occultation curve:

% Setting the extra inputs for the case in which only p,z,u1,u2 are given
if nargin<5
    ToPlot=false;
end;
TolFactor=1;
ToSmooth=true;
for i=1:2:numel(varargin)
    switch lower(varargin{i})
        case 'tolfactor'
            TolFactor=varargin{i+1};
        case 'tosmooth'
            ToSmooth=varargin{i+1};
        otherwise
            fprintf('WARNING: parameter %s is unknown, and so ignored.\n',varargin{i});
    end;
end;

%make z a column vector
if size(z,2)>size(z,1) % z is a column vector
    z=z.';
end;

%set zeros to lambdad, lambda3, 3tad, kap0, kap1
lambdad=zeros(size(z));
lambdae=lambdad;
etad=lambdad;
kap0=lambdad;
kap1=lambdad;

%if p is close to 0.5, make it 0.5
if abs(p-0.5)<1e-3
	p=0.5;
end;

omega=1-u1/3-u2/6;
x1=(p-z).^2;
x2=(p+z).^2;
x3=p.^2-z.^2;

% the source is unocculted:
% Table 3, I.
i=z>=1+p;
if any(i)
    lambdad(i)=0;
    etad(i)=0;
    lambdae(i)=0;
end;

% the source is partly occulted and the occulting object crosses the limb:
% Equation (26):
i=z>=abs(1-p) & z<=1+p;
if any(i)
    kap1(i)=acos(min((1-p^2+z(i).^2)/2./z(i),1));
    kap0(i)=acos(min((p*p+z(i).^2-1)/2/p./z(i),1));
    lambdae(i)=p^2*kap0(i)+kap1(i);
    lambdae(i)=(lambdae(i)-0.5*sqrt(max(4*z(i).^2-(1+z(i).^2-p^2).^2,0)))/pi;
end;

% the occulting object transits the source star (but doesn't
% completely cover it):
i=z<=1-p;
if any(i)
    lambdae(i)=p^2;
end;

% the edge of the occulting star lies at the origin- special 
% expressions in this case:

i=find(abs(z-p)<1e-4*(z+p));
if ~isempty(i)
% Table 3, Case V.:
    j=z(i)>=0.5;
    if any(j)
        %lam=0.5*pi;
        q=0.5/p;
        Kk=ellkVec(q);
        Ek=ellecVec(q);
% Equation 34: lambda_3
        lambdad(i(j))=1/3+16*p/9/pi*(2*p^2-1)*Ek-(32*p^4-20*p^2+3)/9/pi/p*Kk;
% Equation 34: eta_1
        etad(i(j))=1/2/pi*(kap1(i(j))+p^2*(p^2+2*z(i(j))).*kap0(i(j))-(1+5*p^2+z(i(j)).^2)./4.*sqrt((1-x1(i(j))).*(x2(i(j))-1)));
        if p==0.5
% Case VIII: p=1/2, z=1/2
            lambdad(i(j))=1/3-4/pi/9;
            etad(i(j))=3/32;
        end;
    end;
    j=~j;
    if any(j)
% Table 3, Case VI.:
        %lam=0.5*pi;
        q=2*p;
        Kk=ellkVec(q);
        Ek=ellecVec(q);
% Equation 34: lambda_4
        lambdad(i(j))=1/3+2/9/pi*(4*(2*p^2-1)*Ek+(1-4*p^2)*Kk);
% Equation 34: eta_2
        etad(i(j))=p^2/2*(p^2+2*z(i(j)).^2);
    end
end


% the occulting star partly occults the source and crosses the limb:
% Table 3, Case III:
i=find(z>=0.5+abs(p-0.5) & z<1+p | (p>0.5 & z>abs(1-p)*(1+10*eps(1)) & z<p));
if any(i)
    %lam=0.5*pi;
    q=sqrt((1-(p-z(i)).^2)/4./z(i)/p);
    Kk=ellkVec(q);
    Ek=ellecVec(q);
    n=1./x1(i)-1;
    Pk=Kk-n/3.*rj(0,1-q.^2,1,1+n,TolFactor);
% Equation 34, lambda_1:
    lambdad(i)=1/9/pi./sqrt(p*z(i)).*(((1-x2(i)).*(2*x2(i)+x1(i)-3)-3*x3(i).*(x2(i)-2)).*Kk+4*p*z(i).*(z(i).^2+7*p^2-4).*Ek-3*x3(i)./x1(i).*Pk);
    j=z(i)<p;
    if any(j)
        lambdad(i(j))=lambdad(i(j))+2/3;
    end;
% Equation 34, eta_1:
    etad(i)=1/2/pi*(kap1(i)+p^2*(p^2+2*z(i).^2).*kap0(i)-(1+5*p^2+z(i).^2)/4.*sqrt((1-x1(i)).*(x2(i)-1)));
end

% the occulting star transits the source:
% Table 3, Case IV.:
%i=find(p<=1 & z<=(1-p)*1.0001);
i=find(p<=1 & z<=(1-p)*1.0001 & abs(z-p)>10*eps);
if any(i)
    %lam=0.5*pi;
    q=sqrt((x2(i)-x1(i))./(1-x1(i)));
    Kk=ellkVec(q);
    Ek=ellecVec(q);
    n=x2(i)./x1(i)-1;
    Pk=Kk-n/3.*rj(0,1-q.^2,1,1+n,TolFactor);
    % Equation 34, lambda_2:
    lambdad(i)=2/9/pi./sqrt(1-x1(i)).*((1-5*z(i).^2+p^2+x3(i).^2).*Kk+(1-x1(i)).*(z(i).^2+7*p^2-4).*Ek-3*x3(i)./x1(i).*Pk);
    j=z(i)<p;
    if any(j)
        lambdad(i(j))=lambdad(i(j))+2/3;
    end
    j=abs(p+z(i)-1)<=1e-8;
    if any(j)
        lambdad(i(j))=2/3/pi*acos(1-2*p)-4/9/pi*sqrt(p*(1-p))*(3+2*p-8*p*p);
    end
    % Equation 34, eta_2:
    etad(i)=p^2/2.*(p^2+2*z(i).^2);
end

% Now, using equation (33):
muo1=1-((1-u1-2*u2)*lambdae+(u1+2*u2)*lambdad+u2.*etad)/omega;
% Equation 25:
mu0=1-lambdae;

% the  source is completely occulted:
% Table 3, II.
if p>=1
  i=z<=p-1;
  if any(i)
      %lambdad(i)=1;
      %etad(i)=1;
      %lambdae(i)=1;
      muo1(i)=0;
      mu0(i)=0;
  end;
end

if ToSmooth>0
    deg=2;
    CorrectionRange=2e-4;
    FitRange=5*CorrectionRange;
%     i=find(z>=1-p-CorrectionRange & z<=1-p+CorrectionRange); % inner than II/III contact
%     if ~isempty(i)
%         if p<1
%             tmpz=1-p+[linspace(-FitRange,-CorrectionRange,deg*10).'; linspace(CorrectionRange,FitRange,deg*10).'];
%         else
%             tmpz=p-1+[linspace(-FitRange,-CorrectionRange,deg*10).'; linspace(CorrectionRange,FitRange,deg*10).'];
%         end
%         tmpz(tmpz<0)=[];
%         tmp_muo1=occultquadVec(p,tmpz,u1,u2,0,'ToSmooth',false);
%         if ~isempty(tmpz)
%             pol=IterativePolyfit(tmpz,tmp_muo1,deg,3);
%             muo1(i)=polyval(pol,z(i));
%         end
%     end;

    % inner than II/III contact
    i=z>=1-p-CorrectionRange & z<=1-p; 
    if ~isempty(i)
        %if p<1
        %    tmpz=1-p+linspace(-FitRange,-CorrectionRange,deg*10).';
        %else
        %    tmpz=p-1+linspace(-FitRange,-CorrectionRange,deg*10).';
        %end
        tmpz=abs(1-p)+linspace(-FitRange,-CorrectionRange,deg*10).';
        tmpz(tmpz<0)=[];
        tmp_muo1=occultquadVec(p,tmpz,u1,u2,0,'ToSmooth',false);
        if ~isempty(tmpz)
            pol=polyfit(tmpz,tmp_muo1,deg);
            muo1(i)=polyval(pol,z(i));
        end
    end;
    
    % outer than II/III contact
    j=z>=1-p & z<=1-p+CorrectionRange; 
    if ~isempty(j)
        %if p<1
        %    tmpz=1-p+linspace(CorrectionRange,FitRange,deg*10).';
        %else
        %    tmpz=p-1+linspace(CorrectionRange,FitRange,deg*10).';
        %end
        tmpz=abs(1-p)+linspace(CorrectionRange,FitRange,deg*10).';
        tmpz(tmpz<0)=[];
        tmp_muo1=occultquadVec(p,tmpz,u1,u2,0,'ToSmooth',false);
        if ~isempty(tmpz)
            pol=polyfit(tmpz,tmp_muo1,deg);
            muo1(j)=polyval(pol,z(j));
        end
    end;

    % close to z=p
    if p>0.5
        deg=3;
    end;
    CorrectionRange=1e-2;
    FitRange=5*CorrectionRange;
    if p<0.5  
        i=find(z-p>-CorrectionRange & z-p<min(CorrectionRange,1-2*p)); % near z~p
    else
        i=find(z-p>max(-CorrectionRange,1-2*p) & z-p<CorrectionRange); % near z~p
    end
    if ~isempty(i)
        if p<0.5
            tmpz=p+linspace(-FitRange,min(FitRange,1-2*p),deg*10).';
        else
            tmpz=p+linspace(max(-FitRange,1-2*p),FitRange,deg*10).';
        end
        tmpz(abs(tmpz-p)<CorrectionRange | tmpz<0)=[];
        tmp_muo1=occultquadVec(p,tmpz,u1,u2,0,'ToSmooth',false);
        if ~isempty(tmpz)
            pol=polyfit(tmpz,tmp_muo1,deg);
            muo1(i)=polyval(pol,z(i));
        end;
    end;

end


if ToPlot
    plot([flipud(-z); z],[flipud(mu0); mu0],'.',[flipud(-z); z],[flipud(muo1); muo1],'.r');
    a=axis;
    axis([a(1) a(2) a(3) 1+(1-a(3))*0.1]);
    title({'Uniform source (blue) and quadratic model for the limb darkening (red)';sprintf('p=%g, u1=%g, u2=%g',p,u1,u2)});
end;


function rc=rc(x,y,ERRTOL)
xt=zeros(size(x)); yt=xt;
w=xt; s=xt;
alamb=xt; ave=xt;

ERRTOL=.04*ERRTOL;
THIRD=1/3;
C1=.3;
C2=1/7;
C3=.375;
C4=9/22;
i=y>0;
if any(i)
    xt(i)=x(i);
    yt(i)=y(i);
    w(i)=1;
end;
i=~i;
if any(i)
    xt(i)=x(i)-y(i);
    yt(i)=-y(i);
    w(i)=sqrt(x(i))./sqrt(xt(i));
end
i=true(size(x));
while any(i)
    alamb(i)=2*sqrt(xt(i)).*sqrt(yt(i))+yt(i);
    xt(i)=.25*(xt(i)+alamb(i));
    yt(i)=.25*(yt(i)+alamb(i));
    ave(i)=THIRD*(xt(i)+yt(i)+yt(i));
    s(i)=(yt(i)-ave(i))./ave(i);
    i(i)=abs(s(i))>ERRTOL;
end;
rc=w.*(1+s.^2.*(C1+s.*(C2+s.*(C3+s.*C4))))./sqrt(ave);
%  (C) Copr. 1986-92 Numerical Recipes Software

function rj=rj(x,y,z,p,ERRTOL)
% built for x,z scalars and y,p vectors
fac=ones(size(y));
x=x*fac; z=z*fac; rj=0*fac;
xt=rj; yt=rj; zt=rj; pt=rj; a=rj; b=rj; rho=rj; tau=rj; rcx=rj;
sqrtx=rj; sqrty=rj; sqrtz=rj; alamb=rj; alpha=rj; beta=rj; Sum=rj; ave=rj;
ave4=ave*ones(1,4);

ERRTOL=.05*ERRTOL;
C1=3/14;
C2=1/3;
C3=3/22;
C4=3/26;
C5=.75*C3;
C6=1.5*C4;
C7=.5*C2;
C8=C3+C3;
i=p>0;
if any(i)
    xt(i)=x(i);
    yt(i)=y(i);
    zt(i)=z(i);
    pt(i)=p(i);
end;
i=~i;
if any(i)
    xt(i)=min([x(i),y(i),z(i)],[],2);
    zt(i)=max([x(i),y(i),z(i)],[],2);
    yt(i)=x(i)+y(i)+z(i)-xt(i)-zt(i);
    a(i)=1./(yt(i)-p(i));
    b(i)=a(i).*(zt(i)-yt(i)).*(yt(i)-xt(i));
    pt(i)=yt(i)+b(i);
    rho(i)=xt(i).*zt(i)./yt(i);
    tau(i)=p(i).*pt(i)./yt(i);
    rcx(i)=rc(rho(i),tau(i),ERRTOL);
end
%continue                      % ************** 1 ***************
i=true(size(y));
while any(i)
    sqrtx(i)=sqrt(xt(i));
    sqrty(i)=sqrt(yt(i));
    sqrtz(i)=sqrt(zt(i));
    alamb(i)=sqrtx(i).*(sqrty(i)+sqrtz(i))+sqrty(i).*sqrtz(i);
    alpha(i)=(pt(i).*(sqrtx(i)+sqrty(i)+sqrtz(i))+sqrtx(i).*sqrty(i).*sqrtz(i)).^2;
    beta(i)=pt(i).*(pt(i)+alamb(i)).^2;
    Sum(i)=Sum(i)+fac(i).*rc(alpha(i),beta(i),ERRTOL);
    fac(i)=.25*fac(i);
    xt(i)=.25*(xt(i)+alamb(i));
    yt(i)=.25*(yt(i)+alamb(i));
    zt(i)=.25*(zt(i)+alamb(i));
    pt(i)=.25*(pt(i)+alamb(i));
    ave(i)=.2*(xt(i)+yt(i)+zt(i)+pt(i)+pt(i));
    ave4(i,:)=ave(i)*ones(1,4);
    %delx=(ave-xt)/ave;
    %dely=(ave-yt)/ave;
    %delz=(ave-zt)/ave;
    %delp=(ave-pt)/ave;
    del4=(ave4-[xt yt zt pt])./ave4;
    i=max(abs(del4),[],2)>ERRTOL;
end;
ea=del4(:,1).*(del4(:,2)+del4(:,3))+del4(:,2).*del4(:,3);
eb=del4(:,1).*del4(:,2).*del4(:,3);
ec=del4(:,4).^2;
ed=ea-3*ec;
ee=eb+2*del4(:,4).*(ea-ec);
rj=3*Sum+fac.*(1+ed.*(-C1+C5*ed-C6*ee)+eb.*(C7+del4(:,4).*(-C8+del4(:,4).*C4))+del4(:,4).*ea.*(C2-del4(:,4)*C3)-C2*del4(:,4).*ec)./(ave.*sqrt(ave));
i=p<=0;
if any(i)
    rj(i)=a.*(b.*rj(i)+3*(rcx(i)-rf(xt(i),yt(i),zt(i),ERRTOL)));
end;
%  (C) Copr. 1986-92 Numerical Recipes Software

function rf=rf(x,y,z,ERRTOL)
xyz=[x y z];
sqrtxyz=zeros(size(xyz));
alamb=sqrtxyz;
ave=sqrtxyz;
ave3=sqrtxyz;
del3=sqrtxyz;

ERRTOL=.08*ERRTOL;
THIRD=1/3;
C1=1/24;
C2=.1;
C3=3/44;
C4=1/14;
xyzt=xyz;
%continue                  % **************** 1 *****************
i=true(size(x));
while any(i)
    sqrtxyz(i,:)=sqrt(xyz(i,:));
   %alamb(i,:)=sqrtxyz(i,1).*(sqrtxyz(i,2)+sqrtxy(i,3))+sqrtxyz(i,2).*sqrtxyz(i,3);
    alamb(i,:)=sqrtxyz(i,1).*(sqrtxyz(i,2)+sqrtxyz(i,3))+sqrtxyz(i,2).*sqrtxyz(i,3);
    xyzt(i,:)=.25*(xyzt(i,:)+alamb(i,:)*ones(1,3));
    ave(i)=THIRD*(sum(xyzt(i,:),2));
    ave3(i,:)=ave(i)*ones(1,3);
    del3(i,:)=(ave3(i,:)-xyzt(i,:))./ave3(i,:);
    i(i)=max(abs(del3(i,:)),[],2)>ERRTOL;
end;
e2=del3(:,1).*del3(:,2)-del3(:,3)^2;
e3=del3(:,1).*del3(:,2).*del3(:,3);
rf=(1+(C1*e2-C2-C3*e3).*e2+C4*e3)./sqrt(ave);
%  (C) Copr. 1986-92 Numerical Recipes Software 0NL&WR2.