%function bhcoat calculates efficiencies for extinction, 
%total scattering and radar-backscattering for a given size parameters of core and coat and relative
%refractive indices. This code is a translation of a Fortran code of Bohren & Huffman.
%Emmanuel Boss Sep 20 2001

function[Qscat,Qext,Qback]=bhcoat(x,y,refrel1,refrel2);

del=1.0e-8; %inner sphere convergence criterion

x1=refrel1*x;
x2=refrel2*x;
y2=refrel2*y;
ystop=y+4*y^.3333+2;
refrel=refrel2/refrel1;
nstop=ystop;  %series terminated after nstop terms
D0x1=cos(x1)/sin(x1);
D0x2=cos(x2)/sin(x2);
D0y2=cos(y2)/sin(y2);
PSI0y=cos(y);
PSI1y=sin(y);
CHI0y=-sin(y);
CHI1y=cos(y);
XI0y=PSI0y-CHI0y*sqrt(-1);
XI1y=PSI1y-CHI1y*sqrt(-1);
CHI0y2=-sin(y2);
CHI1y2=cos(y2);
CHI0x2=-sin(x2);
CHI1x2=cos(x2);
Qscat=0.0;
Qext=0.0;
XBACK=0+0*sqrt(-1);
n=1;
Iflag=0;
while(n-1-nstop<0)
   RN=n;
   PSIy=(2*RN-1)*PSI1y/y-PSI0y;
   CHIy=(2*RN-1)*CHI1y/y-CHI0y;
   XIy=PSIy-CHIy*sqrt(-1);
   D1y2=1./(RN/y2-D0y2)-RN/y2;
   if (Iflag==1)
       DNBAR=D1y2-BRACK*CHIPy2;
       DNBAR=DNBAR/(1.-BRACK*CHIy2);
       GNBAR=D1y2-CRACK*CHIPy2;
       GNBAR=GNBAR/(1.-CRACK*CHIy2);
       AN=(DNBAR/refrel2+RN/y)*PSIy-PSI1y;
       AN=AN/((DNBAR/refrel2+RN/y)*XIy-XI1y);
       BN=(refrel2*GNBAR+RN/y)*PSIy-PSI1y;
       BN=BN/((refrel2*GNBAR+RN/y)*XIy-XI1y);
       Qscat=Qscat+(2.*RN+1)*(abs(AN)^2+abs(BN)^2);
       XBACK=XBACK+(2*RN+1)*(-1)^n*(AN-BN);
       Qext=Qext+(2*RN+1)*(real(AN)+real(BN));
       PSI0y=PSI1y;
       PSI1y=PSIy;
       CHI0y=CHI1y;
       CHI1y=CHIy;
       XI1y=PSI1y-sqrt(-1)*CHI1y;
       CHI0x2=CHI1x2;
       CHI1x2=CHIx2;
       CHI0y2=CHI1y2;
       CHI1y2=CHIy2;
       D0x1=D1x1;
       D0x2=D1x2;
       D0y2=D1y2;
       n=n+1;
     else
        D1x1=1./(RN/x1-D0x1)-RN/x1;
        D1x2=1./(RN/x2-D0x2)-RN/x2;
        CHIx2=(2*RN-1)*CHI1x2/x2-CHI0x2;
        CHIy2=(2*RN-1)*CHI1y2/y2-CHI0y2;
        CHIPx2=CHI1x2-RN*CHIx2/x2;
        CHIPy2=CHI1y2-RN*CHIy2/y2;
        ANCAP=refrel*D1x1-D1x2;
        ANCAP=ANCAP/(refrel*D1x1*CHIx2-CHIPx2);
        ANCAP=ANCAP/(CHIx2*D1x2-CHIPx2);
        BRACK=ANCAP/(CHIy2*D1y2-CHIPy2);
        BNCAP=refrel*D1x2-D1x1;
        BNCAP=BNCAP/(refrel*CHIPx2-D1x1*CHIx2);
        BNCAP=BNCAP/(CHIx2*D1x2-CHIPx2);
        CRACK=BNCAP/(CHIy2*D1y2-CHIPy2);
        AMESS1=BRACK*CHIPy2;
        AMESS2=BRACK*CHIy2;
        AMESS3=CRACK*CHIPy2;
        AMESS4=CRACK*CHIy2;
        if (abs(AMESS1)>del*abs(D1y2) | abs(AMESS2)>del | abs(AMESS3)>del*abs(D1y2) | abs(AMESS4)>del)
            DNBAR=D1y2-BRACK*CHIPy2;
            DNBAR=DNBAR/(1.-BRACK*CHIy2);
            GNBAR=D1y2-CRACK*CHIPy2;
            GNBAR=GNBAR/(1.-CRACK*CHIy2);
            AN=(DNBAR/refrel2+RN/y)*PSIy-PSI1y;
            AN=AN/((DNBAR/refrel2+RN/y)*XIy-XI1y);
            BN=(refrel2*GNBAR+RN/y)*PSIy-PSI1y;
            BN=BN/((refrel2*GNBAR+RN/y)*XIy-XI1y);
            Qscat=Qscat+(2.*RN+1)*(abs(AN)^2+abs(BN)^2);
            XBACK=XBACK+(2*RN+1)*(-1)^n*(AN-BN);
            Qext=Qext+(2*RN+1)*(real(AN)+real(BN));
            PSI0y=PSI1y;
            PSI1y=PSIy;
            CHI0y=CHI1y;
            CHI1y=CHIy;
            XI1y=PSI1y-sqrt(-1)*CHI1y;
            CHI0x2=CHI1x2;
            CHI1x2=CHIx2;
            CHI0y2=CHI1y2;
            CHI1y2=CHIy2;
            D0x1=D1x1;
            D0x2=D1x2;
            D0y2=D1y2;
            n=n+1;
        else
            BRACK=0+sqrt(-1)*0;
            CRACK=0+sqrt(-1)*0;
            Iflag=1;
        end
    end   
end
Qscat=(2./y^2)*Qscat;
Qext=(2./y^2)*Qext;
%Note: Qback is not Qbb, but the radar back scattering.
Qback=XBACK*conj(XBACK)/y^2;
