Home > atmlab > mie > mie_phasefunction.m

mie_phasefunction

PURPOSE ^

Computation of Mie Phase Function p=p1+p2 (unpolarised)

SYNOPSIS ^

function result = mie_phasefunction(m, x, u, peak)

DESCRIPTION ^

 Computation of Mie Phase Function p=p1+p2 (unpolarised)
 with normalisation to 'one' when integrated 
 over all directions/(4*pi), see Chandrasekhar 1960
 Eq. (28), with complex refractive index m=m'+im", 
 size parameter x=k0*a, and u=cos(scattering angle),
 where k0=vacuum wave number, a=sphere radius;
 u=cos(teta), teta=scattering angle
 peak=0 if diffraction signal is to be subtracted
 s. p. 111-114, Bohren and Huffman (1983) BEWI:TDD122
 C. M�zler, July 2003, revised April 2004.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

mie_phasefunction.m

SOURCE CODE ^

0001 function result = mie_phasefunction(m, x, u, peak)
0002 
0003 % Computation of Mie Phase Function p=p1+p2 (unpolarised)
0004 % with normalisation to 'one' when integrated
0005 % over all directions/(4*pi), see Chandrasekhar 1960
0006 % Eq. (28), with complex refractive index m=m'+im",
0007 % size parameter x=k0*a, and u=cos(scattering angle),
0008 % where k0=vacuum wave number, a=sphere radius;
0009 % u=cos(teta), teta=scattering angle
0010 % peak=0 if diffraction signal is to be subtracted
0011 % s. p. 111-114, Bohren and Huffman (1983) BEWI:TDD122
0012 % C. M�zler, July 2003, revised April 2004.
0013 
0014 nmax=round(2+x+4*x^(1/3));
0015 ab=Mie_ab(m,x);
0016 an=ab(1,:);
0017 bn=ab(2,:);
0018 
0019 pt=Mie_pt(u,nmax);
0020 pin =pt(1,:);
0021 tin =pt(2,:);
0022 
0023 n=(1:nmax);
0024 n2=(2*n+1)./(n.*(n+1));
0025 pin=n2.*pin;
0026 tin=n2.*tin;
0027 S1=(an*pin'+bn*tin');
0028 S2=(an*tin'+bn*pin');
0029 if peak==0,
0030     % Computation of diffraction pattern S according to BH, p. 110
0031     xs=x.*sqrt(1-u.*u);
0032     if abs(xs)<0.0001
0033         S=x.*x*0.25.*(1+u);            % avoiding division by zero
0034     else
0035         S=x.*x*0.5.*(1+u).*besselj(1,xs)./xs;    
0036     end;
0037     S1=S1-S;
0038     S2=S2-S;
0039 end;
0040 Q=mie(m,x);
0041 Qext=Q(1); Qsca=Q(2); asy=Q(5); w0=Qsca/Qext;
0042 p=2*(S1'*S1+S2'*S2)/(Qsca*x^2);  
0043 % Qsca to be exchanged by Qext above if normalisation to
0044 % single-scattering albedo, w, is required
0045 result=[p,w0,asy];

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005