[已答复] mie散射:带有下标的赋值维度不匹配

[复制链接]
123MATLAB呀 发表于 2022-7-27 13:03:26
本帖最后由 123MATLAB呀 于 2022-7-27 16:20 编辑

clc,clear,close;

lamda=1550E-3;  %%  波长:μm
nn=1.811+0.1623*log(lamda)+0.027*(log(lamda)).^2+0.0417*(log(lamda)).^3;  %%  复折射率实部
kk=0.5821+0.1213*log(lamda)+0.2309*(log(lamda)).^2-0.01*(log(lamda)).^3;  %%  复折射率虚部
m=nn+kk*1i;  %%  复折射率

a=0.005+(0.5-0.005)*rand(1,1);  %%  气溶胶粒子半径,随机产生:0.005-0.5μm
x=(2*pi*a)./lamda;  %%  气溶胶粒径参数
theta=rand(1,1)*pi;  %%  散射角
d=0:0.4:4; %%  传输距离
L=11;
syms r
r=0.005:0.0495:0.5; %%  半径
nd=0:4E6:4E7;  %%  粒子数浓度
sigma=1.5; %%  平均标准差
rg=0.06;  %%  平均半径
Nr=nd.*exp(-((log(r)-log(rg))./(sqrt(2*pi).*log(sigma))).^2)./(sqrt(2*pi).*log(sigma));  %%  尺度分布


%---------------------递推公式参数------------------------%
%%% 勒让德函数:π(x)-pin(x)、τ(x)-tao(x)
pin=zeros(L,1);  %%%  初始值
tao=zeros(L,1);
pin(1)=1;
pin(2)=3*cos(theta);
tao(1)=cos(theta);
tao(2)=6*cos(theta)*cos(theta)-3;
    for n=3:L
        pin(n)=((2*n-1)/(n-1))*cos(theta)* pin(n-1)-(n/(n-1))*pin(n-2);
        tao(n)=n*cos(theta)*pin(n)-(n+1)*pin(n-1);
    end
  %%%  黎卡提—贝塞尔函数
  %%%  变量为x的psi、ksi、psi导数、导数/变量为mx的psi、ksi、psi导数、ksi导数
    psi=zeros(L,1);  %%  ψ(x)
    PSI=zeros(L,1);  %%  ψ'(x)导数

    psi(1)=sin(x)/x-cos(x);
    psi(2)=(2/x^2)*sin(x)-(2*cos(x))/x-sin(x);
    PSI(1)=(-1/x^2)*sin(x)+(1/x)*cos(x)+sin(x);
    PSI(2)=(1/x^3)*sin(x)-(1/x^2)*cos(x)-cos(x);
    for n=3:L
        psi(n)=2*(n-1)/x*psi(n-1)-psi(n-2);
        PSI(n)=(-n/x)*psi(n)-psi(n-1);
    end
%%%  黎卡提—汉开尔函数
    ksi=zeros(L,1);  %%  ξ(x)
    KSI=zeros(L,1);  %%  ξ'(x)
    ksi(1)=(1/x)*(sin(x)+1i*cos(x))-(cos(x)-1i*sin(x));
    ksi(2)=(2/x)*ksi(1)-(sin(x)+1i*cos(x));
    KSI(1)=(-1/x)*ksi(1)+ sin(x)+1i*cos(x);
    KSI(2)=(-2/x)*ksi(2)+ ksi(1);
     for n=3:L
        ksi(n)=2*(n-1)/x*ksi(n-1)-ksi(n-2);
        KSI(n)=(-n/x)*ksi(n)-ksi(n-1);
     end

%%%  ψ(mx)、ψ'(mx)
    psim=zeros(L,1);  %%  ψ(mx)
    PSIm=zeros(L,1);  %%  ψ'(mx)导数
    psim(1)=sin(m*x)/x-cos(m*x);
    psim(2)=(2/(m*x)^2)*sin(m*x)-(2*cos(m*x))/(m*x)-sin(m*x);
    PSIm(1)=(-1/(m*x)^2)*sin(m*x)+(1/m*x)*cos(m*x)+sin(m*x);
    PSIm(2)=(1/(m*x)^3)*sin(m*x)-(1/(m*x)^2)*cos(m*x)-cos(m*x);
    for n=3:L
        psim(n)=2*(n-1)/x*psi(n-1)-psi(n-2);
        PSIm(n)=(-n/x)*psi(n)-psi(n-1);
    end

%%%  Mie散射系数:an,bn、消光效率因子Qext、消光系数Kext、链路衰减Qlac
     for n=1:L
            a(n)=(   PSIm(n)*psi(n)-m*psim(n)*PSI(n)   )/...
                (   PSIm(n)*ksi(n)-m*psim(n)*KSI(n)  );
            b(n)=(  m*PSIm(n)*psi(n)-psim(n)*PSI(n)     )/...
                (   m*PSIm(n)*ksi(n)-psim(n)*KSI(n)  );
            Qext(n)=(2/x^2)*(2*n+1)*( real(a(n)+b(n)) );  %%  消光效率因子
               fun=@(r)(r.^2.*Qext(n).*Nr);
               Kext(n)=integral(fun,0.005,0.5);
               vpa(Kext(n));
               Qlac(n)=10*d*Kext(n)*log10(exp(1));        
     end
     surf(Qlac,d,nd);

其中, Kext(n)=integral(fun,0.005,0.5,'ArrayValued',true);显示带有下标的赋值维度不匹配,不知道哪里出现错误,求大神解答


2 条回复


20141303 发表于 2022-7-27 14:30:06
仅供参考

  1. clc,clear,close;

  2. lamda=1550E-3;  %%  波长:μm
  3. nn=1.811+0.1623*log(lamda)+0.027*(log(lamda)).^2+0.0417*(log(lamda)).^3;  %%  复折射率实部
  4. kk=0.5821+0.1213*log(lamda)+0.2309*(log(lamda)).^2-0.01*(log(lamda)).^3;  %%  复折射率虚部
  5. m=nn+kk*1i;  %%  复折射率

  6. a=0.005+(0.5-0.005)*rand(1,1);  %%  气溶胶粒子半径,随机产生:0.005-0.5μm
  7. x=(2*pi*a)./lamda;  %%  气溶胶粒径参数
  8. theta=rand(1,1)*pi;  %%  散射角
  9. d=0:0.4:4; %%  传输距离
  10. L=8;
  11. syms r
  12. r=0.005:0.0495:0.5; %%  半径
  13. nd=0:4E6:4E7;  %%  粒子数浓度
  14. sigma=1.5; %%  平均标准差
  15. rg=0.06;  %%  平均半径
  16. Nr=nd.*exp(-((log(r)-log(rg))./(sqrt(2*pi).*log(sigma))).^2)./(sqrt(2*pi).*log(sigma));  %%  尺度分布


  17. %---------------------递推公式参数------------------------%
  18. %%% 勒让德函数:π(x)-pin(x)、τ(x)-tao(x)
  19. pin=zeros(L,1);  %%%  初始值
  20. tao=zeros(L,1);
  21. pin(1)=1;
  22. pin(2)=3*cos(theta);
  23. tao(1)=cos(theta);
  24. tao(2)=6*cos(theta)*cos(theta)-3;
  25. for n=3:L
  26.     pin(n)=((2*n-1)/(n-1))*cos(theta)* pin(n-1)-(n/(n-1))*pin(n-2);
  27.     tao(n)=n*cos(theta)*pin(n)-(n+1)*pin(n-1);
  28. end
  29. %%%  黎卡提—贝塞尔函数
  30. %%%  变量为x的psi、ksi、psi导数、导数/变量为mx的psi、ksi、psi导数、ksi导数
  31. psi=zeros(L,1);  %%  ψ(x)
  32. PSI=zeros(L,1);  %%  ψ'(x)导数

  33. psi(1)=sin(x)/x-cos(x);
  34. psi(2)=(2/x^2)*sin(x)-(2*cos(x))/x-sin(x);
  35. PSI(1)=(-1/x^2)*sin(x)+(1/x)*cos(x)+sin(x);
  36. PSI(2)=(1/x^3)*sin(x)-(1/x^2)*cos(x)-cos(x);
  37. for n=3:L
  38.     psi(n)=2*(n-1)/x*psi(n-1)-psi(n-2);
  39.     PSI(n)=(-n/x)*psi(n)-psi(n-1);
  40. end
  41. %%%  黎卡提—汉开尔函数
  42. ksi=zeros(L,1);  %%  ξ(x)
  43. KSI=zeros(L,1);  %%  ξ'(x)
  44. ksi(1)=(1/x)*(sin(x)+1i*cos(x))-(cos(x)-1i*sin(x));
  45. ksi(2)=(2/x)*ksi(1)-(sin(x)+1i*cos(x));
  46. KSI(1)=(-1/x)*ksi(1)+ sin(x)+1i*cos(x);
  47. KSI(2)=(-2/x)*ksi(2)+ ksi(1);
  48. for n=3:L
  49.     ksi(n)=2*(n-1)/x*ksi(n-1)-ksi(n-2);
  50.     KSI(n)=(-n/x)*ksi(n)-ksi(n-1);
  51. end

  52. %%%  ψ(mx)、ψ'(mx)
  53. psim=zeros(L,1);  %%  ψ(mx)
  54. PSIm=zeros(L,1);  %%  ψ'(mx)导数
  55. psim(1)=sin(m*x)/x-cos(m*x);
  56. psim(2)=(2/(m*x)^2)*sin(m*x)-(2*cos(m*x))/(m*x)-sin(m*x);
  57. PSIm(1)=(-1/(m*x)^2)*sin(m*x)+(1/m*x)*cos(m*x)+sin(m*x);
  58. PSIm(2)=(1/(m*x)^3)*sin(m*x)-(1/(m*x)^2)*cos(m*x)-cos(m*x);
  59. for n=3:L
  60.     psim(n)=2*(n-1)/x*psi(n-1)-psi(n-2);
  61.     PSIm(n)=(-n/x)*psi(n)-psi(n-1);
  62. end

  63. %%%  Mie散射系数:an,bn、消光效率因子Qext、消光系数Kext、链路衰减Qlac
  64. for n=1:L
  65.     a(n)=(   PSIm(n)*psi(n)-m*psim(n)*PSI(n)   )/...
  66.         (   PSIm(n)*ksi(n)-m*psim(n)*KSI(n)  );
  67.     b(n)=(  m*PSIm(n)*psi(n)-psim(n)*PSI(n)     )/...
  68.         (   m*PSIm(n)*ksi(n)-psim(n)*KSI(n)  );
  69.     Qext(n)=(2/x^2)*(2*n+1)*( imag(a(n))^2+imag(b(n))^2   );  %%  消光效率因子
  70.     fun=@(rr) (rr.^2.*Qext(n).*Nr(n));
  71.     Kext(n)=integral(fun,0.005,0.5);
  72.     vpa(Kext(n));
  73.     Qlac(n)=10*d(n)*4/pi*Kext(n)/log(10);
  74. end
  75. surf(Qlac,d,nd);
复制代码

123MATLAB呀 发表于 4 天前

你好,这个最后算出Qlac是1*11,surf无法画三维
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表