%callbh calculates the size parameter (x) and relative refractive index (refrel)
%for a given sphere refractive index, medium refractive index, radius, and free
%space wavelength. It then calls bhmie.m, the subroutine that computes amplitude
%scattering matrix ellements and efficiencies.
%This routine was written by E. Boss 

clear all

refmed=1;   %input('(real) refractive index of surrounding medium  ');
refre=1.05; %input('(real) refractive index of sphere  ');
refim=0.01; %input('(imaginary) refractive index of sphere  ');

refrel=(refre+sqrt(-1)*refim)/refmed
wavel=.75*.44;%input('wavelength [same units as the radius] in medium (i.e. .75*.65 for 650nm)  ');

Dmin=0.7;%minimal size in population
Dmax=100;%maximal size in population
Nint=35;%number of size classes

fac=(Dmax/Dmin)^(1/(Nint-1));
for i=1:Nint
    D(i)=Dmin*fac^(i-1);
end

x=pi*D*refmed/wavel
nang=91;%input('number of angles between 0 and 90 deg (e.g. 91 for output every degree). needs ot be an odd integer.  ');
[mm nn]=size(x)
nang=91;
dang=pi/2/(nang-1);


%loop over all sizes:
for k=1:nn
    [S1,S2,Qscat(k),Qext(k),Qback]=bhmie(x(k),refrel,nang);
    k
    nan=2*nang-1;
    for j=1:nan
        S11(k,j)=0.5*(abs(S2(j))^2+abs(S1(j))^2);
        ang(j)=dang*(j-1)*180/pi;
    end
    %angles and angles intervals in radians:
    sin_fac=sin(ang*pi/180);
    dtheta=(ang(2)-ang(1))*pi/180;

    %the phase function is proportional to S11, however we have to normalized it:

    ntgrnd=S11(k,:).*sin_fac;

    S11_integral=3/8*(ntgrnd(1)+ntgrnd(nan))+7/6*(ntgrnd(2)+ntgrnd(nan-1))+23/24*(ntgrnd(3)+ntgrnd(nan-2));
    S11_integral=(S11_integral+sum(ntgrnd(4:nan-3)))*2*pi*dtheta;
    S11(k,:)=S11(k,:)/S11_integral; %phase function

    %calculating the back scattering efficiency:
    pid2=(nan+1)/2;
    back_integral=3/8*(ntgrnd(pid2)+ntgrnd(nan))+7/6*(ntgrnd(pid2+1)+ntgrnd(nan-1))+23/24*(ntgrnd(pid2+2)+ntgrnd(nan-2));
    back_integral=(back_integral+sum(ntgrnd(pid2+3:nan-3)))*2*pi*dtheta;
    Qbb(k)=back_integral/S11_integral*Qscat(k);
end
%save data into output file used later for population
save phyto_data


ro=2*x*(refre-1);

figure(1)
subplot(2,1,1)
   semilogx(D,Qscat,'-.',D,Qext,'-+',D,Qext-Qscat,'-o')
   xlabel('D')
   title('Qa-o Qb-. and Qc-+')
subplot(2,1,2)
	semilogx(D,Qbb,'-.')
   title('Qbb')
   xlabel('D')
figure(2)
	semilogx(D,Qbb./Qscat,'-.')
   title('Qbb/Qb')
   xlabel('D')
