Home > atmlab > mie > cloud_phasefunction.m

cloud_phasefunction

PURPOSE ^

Computation of Phase Function pm (unpolarised)

SYNOPSIS ^

function result = cloud_phasefunction(lamda, rc, mu, nsize)

DESCRIPTION ^

 Computation of Phase Function pm (unpolarised)
 of clouds with normalisation to 'one' when integrated 
 over all directions/(4*pi), see Chandrasekhar 1960
 Eq. (28), 
 Input: complex refractive index m=m'+im", 
 mode radius rc, wavelength lamda, both in micron, 
 mu=cos(scattering angle), and nsize=number of drop sizes  
 Output: pm,  
 asym=asymmetry parameter, and wm=single-scattering albedo
 s. p. 111-114, Bohren and Huffman (1983) BEWI:TDD122
 C. M�zler, July 2003

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

cloud_phasefunction.m

SOURCE CODE ^

0001 function result = cloud_phasefunction(lamda, rc, mu, nsize)
0002 
0003 % Computation of Phase Function pm (unpolarised)
0004 % of clouds with normalisation to 'one' when integrated
0005 % over all directions/(4*pi), see Chandrasekhar 1960
0006 % Eq. (28),
0007 % Input: complex refractive index m=m'+im",
0008 % mode radius rc, wavelength lamda, both in micron,
0009 % mu=cos(scattering angle), and nsize=number of drop sizes
0010 % Output: pm,
0011 % asym=asymmetry parameter, and wm=single-scattering albedo
0012 % s. p. 111-114, Bohren and Huffman (1983) BEWI:TDD122
0013 % C. M�zler, July 2003
0014 
0015 dr0=0.5*rc;
0016 dr=1.5*rc/nsize;
0017 mv=1; alfa=6; gam=1;   % Cloud parameters according to Ulaby et al.1981
0018 m=nwater(lamda); 
0019 nx=(1:nsize)';
0020 r=(nx-1)*dr+dr0;
0021 rm=1e-6*r; rcm=1e-6*rc; drm=dr*1e-6;   % transformations to m
0022 xj=2*pi*r./lamda;       % size parameter
0023 sigmag=pi*rm.*rm;      % geometric cross section in m2
0024 NMP=modgammaulaby(rm,rcm,mv,alfa,gam); 
0025 for j=1:nsize,        % computations for increasing drop size
0026     x=xj(j);
0027     nmax=round(2+x+4*x^(1/3));
0028     ab=mie_ab(m,x);
0029     an=ab(1,:);
0030     bn=ab(2,:);
0031     pt=mie_pt(mu,nmax);
0032     pin =pt(1,:);
0033     tin =pt(2,:);
0034     n=(1:nmax);
0035     n2=(2*n+1)./(n.*(n+1));
0036     pin=n2.*pin;
0037     tin=n2.*tin;
0038     S1=(an*pin'+bn*tin');
0039     S2=(an*tin'+bn*pin');
0040     Q=mie(m,x);
0041     Qext(j)=NMP(j)*Q(1); Qsca=Q(2); asy(j)=NMP(j)*Q(5); 
0042     w(j)=NMP(j)*Qsca/Q(1);
0043     p(j)=2*NMP(j)*(S1'*S1+S2'*S2)/(Qsca*x^2);  
0044 end;
0045 sNMP=sum(NMP);
0046 pm=sum(p)/sNMP;
0047 asym=sum(asy)/sNMP;
0048 wm=sum(w)/sNMP;
0049 Qextm=sum(Qext)/sNMP;
0050 % Qsca to be exchanged by Qext above if normalisation to
0051 % single-scattering albedo, w, is required
0052 result=[pm,wm,asym,Qextm];

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