%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, U. of Maine 

clear all

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

refrel=(refre+sqrt(-1)*refim)/refmed

rad=input('radius of sphere  ');
wavel=input('wavelength [same units as the radius] in medium (i.e. .75*.65 for 650nm)  ');

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

nfile=input('name of output file  (e.g. test.dat)  ','s');
fid = fopen(nfile,'w')
fprintf(fid,'real(refrel) imag(refrel)  wavelength radius size-parameter\n');
fprintf(fid,'%8.3f   %8.3f    %8.3f   %8.3f   %8.3f \n',real(refrel),imag(refrel),wavel,rad,x);

%call the function that perform the mie calculation:

[S1,S2,Qscat,Qext,Qback]=bhmie(x,refrel,nang); %note: Qback is not Qbb!
nan=2*nang-1;
for j=1:nan
   S11(j)=0.5*(abs(S2(j))^2+abs(S1(j))^2);
   ang(j)=dang*(j-1)*180/pi;
end

%the phase function is proportional to S11, however we have to normalized it (using a high order formula from Press et al's numerical recipes):

sin_fac=sin(ang*pi/180);
dtheta=(ang(2)-ang(1))*pi/180;
ntgrnd=S11.*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=S11/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=back_integral/S11_integral*Qscat;

fprintf(fid,'Qc  	Qb 		  Qa 	    Qbb\n');
fprintf(fid,'%6.2e  %6.2e  %6.2e  %6.2e\n',Qext,Qscat,Qext-Qscat,Qbb);
fprintf(fid,'Angle  phase-function 	\n');
for j=1:nan
   fprintf(fid,'%6.2f  %6.2e  \n',ang(j),S11(j));
end
fclose(fid);


