mie散射：带有下标的赋值维度不匹配

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);显示带有下标的赋值维度不匹配，不知道哪里出现错误，求大神解答

仅供参考

1. 
   
2. clc,clear,close;
   
3. 
   
4. lamda=1550E-3;  %%  波长：μm
   
5. nn=1.811+0.1623*log(lamda)+0.027*(log(lamda)).^2+0.0417*(log(lamda)).^3;  %%  复折射率实部
   
6. kk=0.5821+0.1213*log(lamda)+0.2309*(log(lamda)).^2-0.01*(log(lamda)).^3;  %%  复折射率虚部
   
7. m=nn+kk*1i;  %%  复折射率
   
8. 
   
9. a=0.005+(0.5-0.005)*rand(1,1);  %%  气溶胶粒子半径，随机产生：0.005-0.5μm
   
10. x=(2*pi*a)./lamda;  %%  气溶胶粒径参数
    
11. theta=rand(1,1)*pi;  %%  散射角
    
12. d=0:0.4:4; %%  传输距离
    
13. L=8;
    
14. syms r
    
15. r=0.005:0.0495:0.5; %%  半径
    
16. nd=0:4E6:4E7;  %%  粒子数浓度
    
17. sigma=1.5; %%  平均标准差
    
18. rg=0.06;  %%  平均半径
    
19. Nr=nd.*exp(-((log(r)-log(rg))./(sqrt(2*pi).*log(sigma))).^2)./(sqrt(2*pi).*log(sigma));  %%  尺度分布
    
20. 
    
21. 
    
22. %---------------------递推公式参数------------------------%
    
23. %%% 勒让德函数：π(x)-pin(x)、τ(x)-tao(x)
    
24. pin=zeros(L,1);  %%%  初始值
    
25. tao=zeros(L,1);
    
26. pin(1)=1;
    
27. pin(2)=3*cos(theta);
    
28. tao(1)=cos(theta);
    
29. tao(2)=6*cos(theta)*cos(theta)-3;
    
30. for n=3:L
    
31.     pin(n)=((2*n-1)/(n-1))*cos(theta)* pin(n-1)-(n/(n-1))*pin(n-2);
    
32.     tao(n)=n*cos(theta)*pin(n)-(n+1)*pin(n-1);
    
33. end
    
34. %%%  黎卡提—贝塞尔函数
    
35. %%%  变量为x的psi、ksi、psi导数、导数/变量为mx的psi、ksi、psi导数、ksi导数
    
36. psi=zeros(L,1);  %%  ψ(x)
    
37. PSI=zeros(L,1);  %%  ψ'(x)导数
    
38. 
    
39. psi(1)=sin(x)/x-cos(x);
    
40. psi(2)=(2/x^2)*sin(x)-(2*cos(x))/x-sin(x);
    
41. PSI(1)=(-1/x^2)*sin(x)+(1/x)*cos(x)+sin(x);
    
42. PSI(2)=(1/x^3)*sin(x)-(1/x^2)*cos(x)-cos(x);
    
43. for n=3:L
    
44.     psi(n)=2*(n-1)/x*psi(n-1)-psi(n-2);
    
45.     PSI(n)=(-n/x)*psi(n)-psi(n-1);
    
46. end
    
47. %%%  黎卡提—汉开尔函数
    
48. ksi=zeros(L,1);  %%  ξ(x)
    
49. KSI=zeros(L,1);  %%  ξ'(x)
    
50. ksi(1)=(1/x)*(sin(x)+1i*cos(x))-(cos(x)-1i*sin(x));
    
51. ksi(2)=(2/x)*ksi(1)-(sin(x)+1i*cos(x));
    
52. KSI(1)=(-1/x)*ksi(1)+ sin(x)+1i*cos(x);
    
53. KSI(2)=(-2/x)*ksi(2)+ ksi(1);
    
54. for n=3:L
    
55.     ksi(n)=2*(n-1)/x*ksi(n-1)-ksi(n-2);
    
56.     KSI(n)=(-n/x)*ksi(n)-ksi(n-1);
    
57. end
    
58. 
    
59. %%%  ψ(mx)、ψ'(mx)
    
60. psim=zeros(L,1);  %%  ψ(mx)
    
61. PSIm=zeros(L,1);  %%  ψ'(mx)导数
    
62. psim(1)=sin(m*x)/x-cos(m*x);
    
63. psim(2)=(2/(m*x)^2)*sin(m*x)-(2*cos(m*x))/(m*x)-sin(m*x);
    
64. PSIm(1)=(-1/(m*x)^2)*sin(m*x)+(1/m*x)*cos(m*x)+sin(m*x);
    
65. PSIm(2)=(1/(m*x)^3)*sin(m*x)-(1/(m*x)^2)*cos(m*x)-cos(m*x);
    
66. for n=3:L
    
67.     psim(n)=2*(n-1)/x*psi(n-1)-psi(n-2);
    
68.     PSIm(n)=(-n/x)*psi(n)-psi(n-1);
    
69. end
    
70. 
    
71. %%%  Mie散射系数:an,bn、消光效率因子Qext、消光系数Kext、链路衰减Qlac
    
72. for n=1:L
    
73.     a(n)=(   PSIm(n)*psi(n)-m*psim(n)*PSI(n)   )/...
    
74.         (   PSIm(n)*ksi(n)-m*psim(n)*KSI(n)  );
    
75.     b(n)=(  m*PSIm(n)*psi(n)-psim(n)*PSI(n)     )/...
    
76.         (   m*PSIm(n)*ksi(n)-psim(n)*KSI(n)  );
    
77.     Qext(n)=(2/x^2)*(2*n+1)*( imag(a(n))^2+imag(b(n))^2   );  %%  消光效率因子
    
78.     fun=@(rr) (rr.^2.*Qext(n).*Nr(n));
    
79.     Kext(n)=integral(fun,0.005,0.5);
    
80.     vpa(Kext(n));
    
81.     Qlac(n)=10*d(n)*4/pi*Kext(n)/log(10);
    
82. end
    
83. surf(Qlac,d,nd);

复制代码

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

你好，这个最后算出Qlac是1*11，surf无法画三维