Home > atmlab > mie > phasefunbeta0.m

phasefunbeta0

PURPOSE ^

"beta0" of phase function in Mie Theory (s. Eq. 9 of

SYNOPSIS ^

function result = phasefunbeta0(m,x,teta0,factor);

DESCRIPTION ^

 "beta0" of phase function in Mie Theory (s. Eq. 9 of
 Meador and Weaver, 1980 J. Atm. Sci. 37, pp. 630-643).
 beta0 is a function of incidence angle teta;
 the whole function is computed, and for a given
 value of teta0, the function value is given as beta00

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

phasefunbeta0.m

SOURCE CODE ^

0001 function result = phasefunbeta0(m,x,teta0,factor);
0002 
0003 % "beta0" of phase function in Mie Theory (s. Eq. 9 of
0004 % Meador and Weaver, 1980 J. Atm. Sci. 37, pp. 630-643).
0005 % beta0 is a function of incidence angle teta;
0006 % the whole function is computed, and for a given
0007 % value of teta0, the function value is given as beta00
0008 
0009 % Input:
0010 % m: refractive index, x: size parameter
0011 % teta0: incidence angle (rad)
0012 % factor: integer >= 1, best for about factor=10.
0013 
0014 % Output: beta00, mu0=cos(teta2), beta0(teta2),
0015 % for 0<=teta2<=pi/2.     C. M�zler, July 2003
0016 
0017 dn0=0.7;    % Teta-Versatz der St�zwerte: 0 bis < 1
0018 nj=ceil((x+1.5)*2.3);
0019 nsteps=factor*nj;
0020 n2=2*nsteps;
0021 dteta=pi/n2;
0022 m1=real(m); m2=imag(m);
0023 nx=(1:n2);
0024 teta=(nx-dn0)*dteta;
0025 u=cos(teta); s=sin(teta);
0026     for j = 1:n2, 
0027         a(:,j)=Mie_S12(m,x,u(j));
0028         SL(j)= real(a(1,j)'*a(1,j));
0029         SR(j)= real(a(2,j)'*a(2,j));
0030         A(j)=(SL(j)+SR(j)).*s(j);
0031     end;
0032 Q=mie(m,x);
0033 Qext=Q(1); Qsca=Q(2); asy=Q(5); w=Qsca/Qext;
0034 AA=dteta*A./Qsca/x.^2;
0035 
0036 L=[];
0037 for jj=1:nj,
0038     x0=legendre(jj-1,u);
0039     x01=x0(1,:);
0040     L=[L,x01'];    % Legendre Polynom(teta), Order jj-1 in Column jj
0041     gj(jj)=x01*AA'; % gj are the coefficients gl in Mead-Weav
0042     gj2(jj)=gj(jj)*(2*jj-1);
0043 end;
0044 for j1=1:n2,
0045     for j2=1:n2,
0046         ab(j1,j2)=gj2.*L(j1,:)*L(j2,:)';
0047     end;
0048 end;
0049 ptest=dteta*s*ab(:,1)   % ptest should be 2
0050 nx2=(1:nsteps);
0051 teta2=teta(1:nsteps);
0052 s2=s(1:nsteps);
0053 c2=cos(teta2);
0054 L2=L(1:nsteps,:);
0055 Q=0.5*dteta*s2*L2;
0056 for j=1:nsteps,
0057     beta0(j)=1-gj2.*L(j,:)*Q';  % beta0(teta0), MeadWeav 1980, Eq. 9
0058 end;
0059 % treatment of boundary values at teta=0 and teta=pi/2:
0060 teta2=[teta2,pi/2]; beta0=[beta0,0.5]; c2=[c2,0]; s2=[s2,1];
0061 nsteps=nsteps+1;
0062 plot(c2,beta0,'k-'), xlabel('mu0'),ylabel('beta0'),
0063 title(sprintf('x=%g, m=%g+%gi',x,m1,m2));
0064 % interpolation to find value beta00=beta0(teta0):
0065 t2=min(find(teta2>teta0));
0066 t1=max(find(teta2<=teta0));
0067 if t2==1,
0068     beta00=beta0(1);
0069 elseif t1==nsteps,
0070     beta00=0.5;
0071 elseif t1==nsteps-1,
0072     dt=teta0-pi/2;
0073     beta01=beta0(nsteps-1);
0074     beta02=0.5; dteta=pi/2-teta2(nsteps);
0075     beta00=beta02+dt*(beta02-beta01)/(pi/2-teta2(nsteps-1));
0076 else,
0077     dt=teta0-teta2(t1);
0078     beta01=beta0(t1);
0079     beta02=beta0(t2);
0080     beta00=beta01+dt*(beta02-beta01)/dteta;
0081 end;
0082 result.beta0=beta00;
0083 result.mu0=c2;
0084 result.beta0fun=beta0;

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