请帮我运行下列matlab的程序,并给出每一行代码的注释

% 基于米氏散射理论计算
% 输入:
% x - 大小参数 =2pi*λ/半径
% refrel - 复合形式的反射指数(相对折射率),例如:1.5+0.02*i;
% nang - S1 和 S2 函数在 0 到 pi/2 范围内的角度
% 输出:
% S1, S2 - 相位函数有关因子
% Qext - 消光效率 Qc 因子
% Qsca - 散射效率 Qb 因子
% Qback 反向散射效率
% gsca- 不对称参数
[s1,s2,qext,qsca,qback,gsca]=bhmie(8.42,1.022+i*1.119e-02,90)
%PSI归一化并画出其随角度变化的图像
PSI0=sum((abs(s1).*abs(s1)+abs(s2).*abs(s2))/2)
PSI=((abs(s1).*abs(s1)+abs(s2).*abs(s2))/2)/PSI0
plot(1:179,PSI)
xlabel('Ψ[deg]')
ylabel('β')
function [s1,s2,qext,qsca,qback,gsca]=bhmie(x,refrel,nang)
% 引入计算过程中的参量
mxnang=1000;
nmxx=150000;
s1=zeros(1,2*mxnang-1);
s2=zeros(1,2*mxnang-1);
d=zeros(1,nmxx);
amu=zeros(1,mxnang);
pi=zeros(1,mxnang);
pi0=zeros(1,mxnang);
pi1=zeros(1,mxnang);
tau=zeros(1,mxnang);
if (nang > mxnang)
disp('error: nang > mxnang in bhmie')
return
end
if (nang < 2)
nang = 2;
end
pii = 4.*atan(1.);
dx = x;
drefrl = refrel;
y = x*drefrl;
ymod = abs(y);
% 求最接近的整数
% 从 NMX 向下计算的对数导数
xstop = x + 4.*x^0.3333 + 2.;
nmx = max(xstop,ymod) + 15;
nmx=fix(nmx);
nstop = xstop;
if (nmx > nmxx)
'error: nmx > nmxx=', nmxx, ' for |m|x=', ymod
return
end
%
% 为了计算散射强度引入散射角
dang = 0.;
if (nang > 1)
dang = .5*pii/ (nang-1);
end
for j=1: nang
theta = (j-1)*dang;
amu(j) = cos(theta);
end
%π0=0,π1=1
for j=1: nang
pi0(j) = 0.;
pi1(j) = 1.;
end
nn = 2*nang - 1;
%
%通过向下递归计算对数导数 D(j),从初始值(0,0)开始,j=nmx
nn = nmx - 1;
for n=1: nn
en = nmx - n + 1;
d(nmx-n) = (en/y) - (1./ (d(nmx-n+1)+en/y));
end
%
% 通过实参数 X 由向上递归计算Riccati-Bessel函数Ψn和ξn
psi0 = cos(dx);
psi1 = sin(dx);
chi0 = -sin(dx);
chi1 = cos(dx);
xi1 = psi1-chi1*i;
qsca = 0.;
gsca = 0.;
p = -1;
for n=1: nstop
en = n;
fn = (2.*en+1.)/ (en* (en+1.));
% 对于给定的n, psi = psi_n chi = chi_n
% psi1 = psi_{n-1} chi1 = chi_{n-1}
% psi0 = psi_{n-2} chi0 = chi_{n-2}
% 计算 psi_n and chi_n
psi = (2.*en-1.)*psi1/dx - psi0;
chi = (2.*en-1.)*chi1/dx - chi0;
xi = psi-chi*i;
%
%*** 存储之前的an和bn以供使用
% 用来计算g=<cos(theta)>
if (n > 1)
an1 = an;
bn1 = bn;
end
%
%*** 计算an和bn:
an = (d(n)/drefrl+en/dx)*psi - psi1;
an