才发现没贴全
@熬夜修头仙子= an / ((d(n)/drefrl+en/dx)xi - xi1);
bn = (drefrld(n)+en/dx)psi - drefrlpsi1;
bn = bn / (drefrld(n)xi + en/dxxi1);
%** 计算散射幅度
s1(n) = (anpsi0 - an1psi)/ (bnpsi0 - bn1psi);
s2(n) = (anchi0 - an1chi)/ (bnchi0 - bn1chi);
%
%*** 计算消光和散射效率因子, 以及相移
if (n > 1)
qsca = qsca + (2.en+1.) (abs(an)*abs(an)+abs(bn)abs(bn));
gsca = gsca + ((2.en+1.)/ (en (en+1.))) (real(an)*real(bn)+imag(an)imag(bn));
end
psi1 = psi;
psi0 = psi1;
chi1 = chi;
chi0 = chi1;
xi1 = psi1 - chi1i;
%
% 对于入射角度theta计算π_n和τ_n,需要递归计算
if (n <= nang)
pi(n) = pi1(n-1);
pi1(n-1) = ((2.*n-1.)/(n-1.))amu(n)pi(n);
if (n > 1)
tau(n-1) = (n-1.)pi(n-1) - npi0(n-1);
end
end
end
%
% 计算反向散射系数和不对称参数
qsca = (2./ (xx))qsca;
gsca = (4./ (xx))gsca/qsca;
qext = (4./ (xx))real(s1(1));
qback = (4./ (xx))abs(s1(nmx))^2;
%
% 计算S1和S2的幅度相对值
s1 = s1/ (drefrldrefrl);
s2 = s2/ (drefrldrefrl);
end