mie_tetascanallold

 Computation and plot of Mie Power Scattering and diffraction functions 
 for complex refractive-index ratio m=m'+im", size parameters x=k0*a, 
 according to Bohren and Huffman (1983) BEWI:TDD122
 1) polar diagram linear or in dB scale with respect to minimum, with
 SL in upper semicircle, SR in lower semicircle and 3 cartesian diagrams
 2) same for SL0 and SR0 without diffraction pattern, 
 3) scattered intensity S (lin or log scale), and degree of polarisation
 4) scattered intensity without diffraction peak S0 (lin or log scale),
 5) beam efficiencies of S and S0.
 nsteps: number of scattering angles
 nsmooth: number of values to be averaged in polarisation, S and S0
 type:= 'log' or 'lin' for logarithmic or linear plots
 C. M�zler, April 2004.

0001 function result = Mie_tetascanall(m, x, nsteps, nsmooth, type)
0002 
0003 % Computation and plot of Mie Power Scattering and diffraction functions
0004 % for complex refractive-index ratio m=m'+im", size parameters x=k0*a,
0005 % according to Bohren and Huffman (1983) BEWI:TDD122
0006 % 1) polar diagram linear or in dB scale with respect to minimum, with
0007 % SL in upper semicircle, SR in lower semicircle and 3 cartesian diagrams
0008 % 2) same for SL0 and SR0 without diffraction pattern,
0009 % 3) scattered intensity S (lin or log scale), and degree of polarisation
0010 % 4) scattered intensity without diffraction peak S0 (lin or log scale),
0011 % 5) beam efficiencies of S and S0.
0012 % nsteps: number of scattering angles
0013 % nsmooth: number of values to be averaged in polarisation, S and S0
0014 % type:= 'log' or 'lin' for logarithmic or linear plots
0015 % C. M�zler, April 2004.
0016 
0017 nstart=min(nsteps,round(nsteps*pi/x)+nsmooth);
0018 m1=real(m); m2=imag(m);
0019 nx=(1:nsteps); dteta=pi/nsteps;
0020 Q=mie(m,x);  Qext=Q(1); Qsca=Q(2); Qabs=Q(3); Qb=Q(4); asy=Q(5);
0021 nmax=round(2+x+4*x^(1/3));
0022 ab=mie_ab(m,x); an=ab(1,:); bn=ab(2,:);
0023 teta=(nx-0.5).*dteta; tetad=teta*180/pi;
0024 u=cos(teta); s=sin(teta);
0025 px=pi*x^2;
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'); S2=(an*tin'+bn*pin');
0034     xs=x.*s(j);
0035     if abs(xs)<0.00002          % Diffraction pattern according to BH, p. 110
0036         S3=x.*x*0.25.*(1+u(j));    % avoiding division by zero
0037     else
0038         S3=x.*x*0.5.*(1+u(j)).*besselj(1,xs)./xs;    
0039     end;
0040     S4=S1-S3;
0041     S5=S2-S3;
0042     SL(j)= real(S1'*S1)/px;
0043     SR(j)= real(S2'*S2)/px;
0044     SD(j)= real(S3'*S3)/px;        
0045     SL0(j)=real(S4'*S4)/px;        
0046     SR0(j)=real(S5'*S5)/px;        
0047 end;
0048 st=2*pi*s;
0049 SSL=st.*SL;
0050 SSR=st.*SR;
0051 SSLb=st.*SL0;
0052 SSRb=st.*SR0;
0053 z=0.5*dteta*cumsum(SSL+SSR)/Qsca;    % full beam efficiency
0054 z0=0.5*dteta*cumsum(SSLb+SSRb)/Qsca; % diffraction-free beam efficiency
0055 g=0.5*dteta*(SSL+SSR)*u'/Qsca;       % full asymmetry factor
0056 g0=0.5*dteta*(SSLb+SSRb)*u'/Qsca;    % diffraction-free asymmetry factor
0057 
0058 figure;
0059 xmin=min(tetad/180);
0060 semilogx(tetad/180,z,'r-',tetad/180,z0,'k--',tetad/180,z-z0,'b:'),
0061 %title(sprintf('Cumulative Fraction of Scattered Power: m=%g+%gi, x=%g',m1,m2,x)),
0062 xlabel('Maximum Scattering Angle/180�'),
0063 axis([xmin, 1, 0, 1.1]);
0064 
0065 S=SL+SR; S0=SL0+SR0;            % Intensity
0066 Ss=smooth(S,nsmooth); S0s=smooth(S0,nsmooth);
0067 dS=(SL-SR)./S;                  % Degree of polarisation
0068 dSs=smooth(dS,nsmooth);
0069 testa=tetad(nstart:nsteps);
0070 
0071 figure;
0072 if type=='lin'                   % linear plots
0073     y=[teta teta+pi;SL SR(nsteps:-1:1)]'; 
0074     polar(y(:,1),y(:,2)),
0075     title(sprintf('Mie Scattering Diagram: m=%g+%gi, x=%g',m1,m2,x)),
0076     xlabel('Scattering Angle'),
0077 figure;
0078     y=[teta teta+pi;SL0 SR0(nsteps:-1:1)]'; 
0079     polar(y(:,1),y(:,2)),
0080     title(sprintf('No-Peak Scattering Diagram: m=%g+%gi, x=%g',m1,m2,x)),
0081     xlabel('Scattering Angle'),
0082 figure;
0083 subplot(2,1,1);
0084     plot(tetad,Ss,'k-')
0085     title(sprintf('Mie Angular Scattering: m=%g+%gi, x=%g',m1,m2,x)),
0086     xlabel('Scattering Angle'),
0087     ylabel('S');
0088 subplot(2,1,2);
0089     plot(tetad,dSs,'k-')
0090     xlabel('Scattering Angle'),
0091     ylabel('Polarisation Degree ');
0092 figure;
0093     plot(tetad,S0s,'k-')
0094     title(sprintf('No-Peak Angular Scattering: m=%g+%gi, x=%g',m1,m2,x)),
0095     xlabel('Scattering Angle'),
0096     ylabel('S0');
0097 elseif type=='log',              % logar. plots
0098     y=[teta teta+pi;10*log10(SL) 10*log10(SR(nsteps:-1:1))]';
0099     ymin=min(y(:,2));  % Minimum for normalisation of log-polar plot
0100     y(:,2)=y(:,2)-ymin;                  
0101     polar(y(:,1),y(:,2)),
0102     title(sprintf('Mie Scattering Diagram: m=%g+%gi, x=%g, min(dB)= %g',m1,m2,x,ymin)),
0103     xlabel('Scattering Angle (deg)');
0104     y=[teta teta+pi;10*log10(SL0) 10*log10(SR0(nsteps:-1:1))]';
0105     ymin=min(y(:,2));  % Minimum for normalisation of log-polar plot
0106     y(:,2)=y(:,2)-ymin;                  
0107 figure;
0108     polar(y(:,1),y(:,2)),
0109     title(sprintf('No-Peak Scattering Diagram: m=%g+%gi, x=%g, min(dB)= %g',m1,m2,x,ymin)),
0110     xlabel('Scattering Angle (deg)');
0111 figure;
0112 subplot(2,1,1);
0113     semilogy(tetad,Ss,'k-'),    
0114     title(sprintf('Mie Angular Scattering: m=%g+%gi, x=%g',m1,m2,x)),
0115     xlabel('Scattering Angle'),
0116     ylabel('S');
0117 subplot(2,1,2);
0118     plot(tetad,dSs,'k-'),
0119     xlabel('Scattering Angle'),
0120     ylabel('Polarisation Degree ');
0121 figure;
0122     semilogy(testa,Ss(nstart:nsteps),'r:',tetad,S0s,'k-'),    
0123     title(sprintf('No-Peak Angular Scattering: m=%g+%gi, x=%g',m1,m2,x)),
0124     xlabel('Scattering Angle'),
0125     ylabel('S0');
0126 end;
0127 etab=z(nsteps); etab0=z0(nsteps); Qd=Qsca*(1-etab0);
0128 result.s=[tetad',SR',SL',SR0',SL0',SD',dS',z',z0']; 
0129 result.Q=[Qext,Qsca,Qabs,Qb,Qd];
0130 result.etab=[etab,etab0]; 
0131 result.g=[asy, g, g0];