Home > atmlab > mie > mie_phasefunplot.m

mie_phasefunplot

PURPOSE ^

Plot of Mie phasefunction for unpolarised radiation

SYNOPSIS ^

function result = mie_phasefunplot(m, x, nsteps, peak)

DESCRIPTION ^

 Plot of Mie phasefunction for unpolarised radiation 
 with normalisation to 'one' when integrated 
 over all directions/(4*pi), see Chandrasekhar 1960
 Eq. (28), for complex refractive-index ratio m=m'+im", 
 size parameters x=k0*a, 
 according to Bohren and Huffman (1983) BEWI:TDD122
 INPUT: 
 m, x: as usual
 nsteps: number of angles to be considered
 peak=1 for normal case, peak=0 without forward peak
 C. M�zler, July 2003, revised March 2004 (peak).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

mie_phasefunplot.m

SOURCE CODE ^

0001 function result = mie_phasefunplot(m, x, nsteps, peak)
0002 
0003 % Plot of Mie phasefunction for unpolarised radiation
0004 % with normalisation to 'one' when integrated
0005 % over all directions/(4*pi), see Chandrasekhar 1960
0006 % Eq. (28), for complex refractive-index ratio m=m'+im",
0007 % size parameters x=k0*a,
0008 % according to Bohren and Huffman (1983) BEWI:TDD122
0009 % INPUT:
0010 % m, x: as usual
0011 % nsteps: number of angles to be considered
0012 % peak=1 for normal case, peak=0 without forward peak
0013 % C. M�zler, July 2003, revised March 2004 (peak).
0014 
0015 dteta=pi/nsteps;
0016 m1=real(m); m2=imag(m);
0017 nx=(1:nsteps);
0018 teta=(nx-0.5).*dteta;
0019 u=cos(teta);
0020     for j = 1:nsteps, 
0021         if peak==1
0022             a(:,j)=mie_S12(m,x,u(j));
0023         elseif peak==0
0024             a(:,j)=mie_S12nopeak(m,x,u(j));
0025         end;
0026         SL(j)= real(a(1,j)'*a(1,j));
0027         SR(j)= real(a(2,j)'*a(2,j));
0028     end;
0029 tetad=teta*180/pi;
0030 s=sin(teta);
0031 Q=mie(m,x); Qext=Q(1); Qsca=Q(2); Qabs=Q(3); asy=Q(5); w=Qsca/Qext;
0032 p=(SL+SR)./Qsca/x.^2;
0033 eta=dteta*(p*s'-(p(1)*s(1)+p(nsteps)*s(nsteps))*0.5);
0034 asy0=dteta*(p*(u.*s)'-(p(1)*s(1)-p(nsteps)*s(nsteps))*0.5)/eta;
0035 Qsca0=Qsca*eta;Qext0=Qabs+Qsca0;
0036 w0=Qsca0/Qext0;
0037 pdB=10*log10(p);
0038 % Qsca to be exchanged by Qext=Q(1) above if normalisation
0039 % to single-scattering albedo is required.
0040 plot(tetad,pdB,'r-')
0041 title(sprintf('Phase Function: m=%g+%gi, x=%g, w=%g, g=%g',m1,m2,x,w0,asy0))
0042 xlabel('Scattering Angle (deg)'),ylabel('Phase Function (dB)')
0043 result=[teta;p]'; % Phase function

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