miecoated.m

来自「MieMatlabMaetzler.zip 非常全功能的Mie散射程序 ma」· M 代码 · 共 57 行

function result = Miecoated(m1,m2,x,y,opt)

% Mie Efficiencies of coated spheres for given complex refractive-index
% ratios m1=m1'+im1", m2= m2'+im2" of kernel and coating, resp.,
% and size parameters x=k0*a, y=k0*b where k0= wave number in ambient 
% medium, a,b= inner,outer sphere radius, using complex Mie Coefficients
% an and bn for n=1 to nmax,
% s. Bohren and Huffman (1983) BEWI:TDD122, p. 181-185,483.
% Result: m1',m1", m2',m2", x,y, efficiencies for extinction (qext), 
% scattering (qsca), absorption (qabs), backscattering (qb), 
% asymmetry parameter (asy=<costeta>) and (qratio=qb/qsca).
% opt selects the function "Miecoated_ab.." for an and bn, n=1 to nmax.
% Note that 0<=x<=y;   C. M鋞zler, August 2002.

if x==y                     % To reduce computing time
    result=mie(m1,y);
elseif x==0                 % To avoid a singularity at x=0
    result=mie(m2,y);
elseif m1==m2
    result=mie(m1,y);       % To reduce computing time
elseif x>0                  % This is the normal situation
    nmax=round(2+y+4*y.^(1/3));
    n1=nmax-1;
    n=(1:nmax);cn=2*n+1; c1n=n.*(n+2)./(n+1); c2n=cn./n./(n+1);
    y2=y.*y;
if opt==1
    f=miecoated_ab1(m1,m2,x,y);
elseif opt==2
    f=miecoated_ab2(m1,m2,x,y);
elseif opt==3
    f=miecoated_ab3(m1,m2,x,y);
end;
    anp=(real(f(1,:))); anpp=(imag(f(1,:)));
    bnp=(real(f(2,:))); bnpp=(imag(f(2,:)));
    g1(1:4,nmax)=[0; 0; 0; 0]; % displaced numbers used for
    g1(1,1:n1)=anp(2:nmax);    % asymmetry parameter, p. 120
    g1(2,1:n1)=anpp(2:nmax);
    g1(3,1:n1)=bnp(2:nmax);
    g1(4,1:n1)=bnpp(2:nmax);   
    dn=cn.*(anp+bnp);
    q=sum(dn);
    qext=2*q./y2;
    en=cn.*(anp.*anp+anpp.*anpp+bnp.*bnp+bnpp.*bnpp);
    q=sum(en);
    qsca=2*q./y2;
    qabs=qext-qsca;
    fn=(f(1,:)-f(2,:)).*cn;
    gn=(-1).^n;
    f(3,:)=fn.*gn;
    q=sum(f(3,:));
    qb=q*q'./y2;
    asy1=c1n.*(anp.*g1(1,:)+anpp.*g1(2,:)+bnp.*g1(3,:)+bnpp.*g1(4,:));
    asy2=c2n.*(anp.*bnp+anpp.*bnpp);
    asy=4/y2*sum(asy1+asy2)/qsca;
    qratio=qb/qsca;
    result=[qext qsca qabs qb asy qratio];
end;