| 查看: 401 | 回复: 2 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
求助，程序怎么改？运行有错误。
|
||
|
%光振幅函数S1和S2程序% function result=Mie_S12(m,x,u) m=1.5; x=8; u=0.5; nmax=round(2+x+4*x^(1/3)); n=(1:nmax);nu=(n+0.5);z=m.^x;m2=m.*m; sqx=sqrt(0.5*pi./x);sqz=sqrt(0.5*pi./z); bx=besselj(nu,x).*sqx; bz=besselj(nu,z).*sqz; yx=bessely(nu,x).*sqx; hx=bx+i*yx; blx=[sin(x)/x,bx(1:nmax-1)]; blz=[sin(z)/z,bz(1:nmax-1)]; ylx=[-cos(x)/x,yx(1:nmax-1)]; hlx=blx+i*ylx; ax=x.*blx-n.*bx; az=z.*blz-n.*bz; ahx=x.*hlx-n.*hx; an=(m2.*bz.*ax-bx.*az)./(m2.*bz.*ahx-hx.*az); bn=(bz.*ax-bx.*az)./(bz.*ahx-hx.*az); %-1<=u<=1,u=cos(theta) p(1)=1; t(1)=u; p(2)=3*u; t(2)=3*cos(2*(acos(u))); for n1=3:nmax, p1=(2*n1-1)./(n1-1).*p(n1-1).*u; p2=n1./(n1-1).*p(n1-2); p(n1)=p1-p2; t1=n1*u.*p(n1); t2=(n1+1).*p(n1-1); t(n1)=t1-t2; end; ab=Mie_ab(m,x); an=ab(1,  ;bn=ab(2, ;pt=Mie_pt(u,nmax); pin=pt(1, ;tin=pt(2, ;n=(1:nmax); n2=(2*n+1)./(n.*(n+1)); pin=n2.*pin; tin=n2.*tin; S1=(an*pin'+bn*tin'); S2=(an*tin'+bn*pin'); result=[S1;S2]; |
回复此楼» 猜你喜欢
小木虫上这么多卖论文的，真有人买论文么？感觉没必要啊
已经有4人回复
-
为什么网上很多人说本周 12号出结果
已经有4人回复
-
基金中了
已经有15人回复
-
帮忙看看fileCode
已经有8人回复
-
FileCode能看出啥？
已经有15人回复
-
filecode
已经有8人回复
-
确定了，国自然21号放榜
已经有7人回复
-
综述论文作为代表作会不会影响评审专家的印象分？
已经有11人回复
-
关于代码变化问题，想知道的进来
已经有24人回复
-
静等基金结果
已经有16人回复
3楼2017-02-25 15:12:37
2楼2015-12-29 23:48:49
100