Home > atmlab > mie > mie_teta.m

mie_teta

PURPOSE ^

Computation of Mie Efficienicies and gi coefficients of Legendre

SYNOPSIS ^

function result = Mie_teta(m, x, tetadlim)

DESCRIPTION ^

 Computation of Mie Efficienicies and gi coefficients of Legendre 
 Polynomial decomposition for complex refractive-index ratio m=m'+im",
 size parameters x=k0*a, scattered intensity with and without 
 diffraction peak. 
 tetadlim (deg): Optional value of integration limit 
 for beam and diffraction efficiencies (default 180�)
 Output.etab: etab(180�), etab0(180�), etab(tetadlim)
 Output.Q: Qext, Qsca, Qabs, Qb, Qd, Qdlim, asy
 Output.gi and Output.g0i Legendre Coefficients
 C. M�zler, April 2004.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

mie_teta.m

SOURCE CODE ^

0001 function result = Mie_teta(m, x, tetadlim)
0002 
0003 % Computation of Mie Efficienicies and gi coefficients of Legendre
0004 % Polynomial decomposition for complex refractive-index ratio m=m'+im",
0005 % size parameters x=k0*a, scattered intensity with and without
0006 % diffraction peak.
0007 % tetadlim (deg): Optional value of integration limit
0008 % for beam and diffraction efficiencies (default 180�)
0009 % Output.etab: etab(180�), etab0(180�), etab(tetadlim)
0010 % Output.Q: Qext, Qsca, Qabs, Qb, Qd, Qdlim, asy
0011 % Output.gi and Output.g0i Legendre Coefficients
0012 % C. M�zler, April 2004.
0013 
0014 nsteps=round(23*x);     % number of angular steps
0015 if nargin==2,
0016     tetadlim=180;       % default value 180�
0017 end;
0018 nx=(1:nsteps); dteta=pi/nsteps;
0019 Q=mie(m,x);  Qext=Q(1); Qsca=Q(2); Qabs=Q(3); Qb=Q(4); asy=Q(5);
0020 nmax=round(2+x+4*x^(1/3));
0021 ab=mie_ab(m,x); an=ab(1,:); bn=ab(2,:);
0022 teta=(nx-0.5).*dteta; tetad=teta*180/pi;
0023 u=cos(teta); s=sin(teta);
0024 px=pi*x^2;
0025 st=pi*s*dteta/Qsca;      % Constant factor of angular integrands
0026 for j = 1:nsteps, 
0027     pt=mie_pt(u(j),nmax);
0028     pin =pt(1,:);
0029     tin =pt(2,:);
0030     n=(1:nmax);
0031     n2=(2*n+1)./(n.*(n+1));
0032     pin=n2.*pin; tin=n2.*tin;
0033     S1=(an*pin'+bn*tin'); 
0034     S2=(an*tin'+bn*pin');
0035     xs=x.*s(j);
0036     if abs(xs)<0.00002   % Diffraction pattern according to BH, p. 110
0037         S3=x.*x*0.25.*(1+u(j));    % avoiding division by zero
0038     else
0039         S3=x.*x*0.5.*(1+u(j)).*besselj(1,xs)./xs;    
0040     end;
0041     S4=S1-S3;
0042     S5=S2-S3;
0043     SR(j)= real(S1'*S1)/px;
0044     SL(j)= real(S2'*S2)/px;
0045     SR0(j)=real(S4'*S4)/px;        
0046     SL0(j)=real(S5'*S5)/px;        
0047 end;
0048 z=st.*(SL+SR);       % Integrand
0049 z0=st.*(SL0+SR0);    % Integrand
0050 etabdlim=cumsum(z);  % Beam Efficiency with limited angle
0051 nj=11;               
0052 for jj=1:nj,         % Phase fct. decomposition in Legendre Polynomials
0053     xa=legendre(jj-1,u);  % Legendre Function
0054     x0=xa(1,:);      % Legendre Polynomial
0055     gi(jj)=x0*z';    % Beam Eff., asymmetry factor and higher gi's
0056     g0i(jj)=x0*z0';  % same, but diffraction signal removed
0057 end;
0058 etab=gi(1);     gi=gi/etab;     gi=gi(2:nj);
0059 etab0=g0i(1);   g0i=g0i/etab0;  g0i=g0i(2:nj);
0060 Qd=Qsca*(etab-etab0);      % Qd = diffraction efficiency
0061 n=max(find(tetad<tetadlim)); 
0062 Qdlim=Qsca*etabdlim(n);    % Qd = limited diffraction efficiency
0063 
0064 result.etab=[etab, etab0, etabdlim(n)];
0065 result.Q=[Qext, Qsca, Qabs, Qb, Qd, Qdlim, asy];
0066 result.gi=gi;
0067 result.g0i=g0i;

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