Home > atmlab > scattering > mie_phase_mat.m

mie_phase_mat

PURPOSE ^

MIE_PHASE_MAT Mie phase matrix

SYNOPSIS ^

function [y] = mie_phase_mat(m,lambda,r,theta,stokes_dim)

DESCRIPTION ^

 MIE_PHASE_MAT   Mie phase matrix

   Returns the phase matrix for Mie particles. Multiple particle sizes and
   scattering angles are handled.

   The returned phase matrix is not normalised. That is, the upper left
   element of the matrix, y(1,:,1,1), is related to the scattering
   cross-section as:

     int_{4pi}[ y(1,:,1,1) ] domega = y_sca

   where the integration is performed over the complete sphere, omega is solid
   angle and y_sca is the scattering cross-section.

 FORMAT  y = mie_phase_mat(m,lambda,r,theta[,stokes_dim])

 OUT     y          Mie phase matrix. Size(r,theta,stokes_dim,stokes_dim)

 IN      m          complex refractive index (imaginary part should 
                    be positive)   (scalar)
         lambda     wavelength [m] (scalar)
         theta      scattering angle [deg], (vector)
         r          particle radius [m] (vector)   
 OPT     stokes_dim Stokes dimension.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

mie_phase_mat.m

SOURCE CODE ^

0001 % MIE_PHASE_MAT   Mie phase matrix
0002 %
0003 %   Returns the phase matrix for Mie particles. Multiple particle sizes and
0004 %   scattering angles are handled.
0005 %
0006 %   The returned phase matrix is not normalised. That is, the upper left
0007 %   element of the matrix, y(1,:,1,1), is related to the scattering
0008 %   cross-section as:
0009 %
0010 %     int_{4pi}[ y(1,:,1,1) ] domega = y_sca
0011 %
0012 %   where the integration is performed over the complete sphere, omega is solid
0013 %   angle and y_sca is the scattering cross-section.
0014 %
0015 % FORMAT  y = mie_phase_mat(m,lambda,r,theta[,stokes_dim])
0016 %
0017 % OUT     y          Mie phase matrix. Size(r,theta,stokes_dim,stokes_dim)
0018 %
0019 % IN      m          complex refractive index (imaginary part should
0020 %                    be positive)   (scalar)
0021 %         lambda     wavelength [m] (scalar)
0022 %         theta      scattering angle [deg], (vector)
0023 %         r          particle radius [m] (vector)
0024 % OPT     stokes_dim Stokes dimension.
0025 
0026 % History: 2010-02-02 Created by Bengt Rydberg
0027 
0028 function [y] = mie_phase_mat(m,lambda,r,theta,stokes_dim)
0029 %
0030 if nargin < 5
0031   stokes_dim = 1;
0032 end
0033 %                                                                           %&%
0034 rqre_datatype( m, @istensor0 );                                             %&%
0035 rqre_datatype( lambda, @istensor0 );                                        %&%
0036 rqre_alltypes( r, {@isnumeric,@isvector} );                                 %&%
0037 rqre_alltypes( theta, {@isnumeric,@isvector} );                             %&%
0038 rqre_alltypes( stokes_dim, {@istensor0,@iswhole} );                         %&%
0039 rqre_in_range( stokes_dim, 1 , 4 );                                         %&%
0040 if imag( m ) < 0                                                            %&%
0041   error( 'Imaginary part of the refractive index must be positive.' );      %&%
0042 end                                                                         %&%
0043 
0044 DEG2RAD = constants( 'DEG2RAD' );
0045 
0046 k = 2 * pi / lambda;
0047 x = k * r;
0048 
0049 y = zeros( length(r), length(theta), stokes_dim, stokes_dim );
0050 
0051 
0052 %loop over particle sizes
0053 for i=1:length(x)
0054   
0055   %loop over scattering angle
0056   for l=1:length(theta)
0057     u          = cos(DEG2RAD*theta(l));
0058     S          = mie_S12(m,x(i),u);
0059     S11        = 1/2 * ( abs(S(1))^2 + abs(S(2))^2 )/k^2;
0060     y(i,l,1,1) = S11;
0061     %
0062     if stokes_dim > 1        
0063       S12        = 1/2 * ( abs(S(2))^2 - abs(S(1))^2 )/k^2;
0064       y(i,l,2,2) = S11;
0065       y(i,l,1,2) = S12;
0066       y(i,l,2,1) = S12;
0067       %
0068       if stokes_dim > 2
0069         S33        = 1/2 * ( S(2)*conj(S(1)) + S(1)*conj(S(2)))/k^2;
0070         y(i,l,3,3) = S33;
0071         %
0072         if stokes_dim > 3
0073           S34        = -1/2 *imag(-S(2)*conj(S(1)) + S(1)*conj(S(2))) /k^2;
0074           y(i,l,4,4) = S33;
0075           y(i,l,3,4) = S34;
0076           y(i,l,4,3) = -S34;
0077         end
0078       end
0079     end
0080   end
0081 end
0082

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