[讨论]球散射的外推程序，请大家帮忙看看

最近编了一段球的散射的程序，大家帮我看下我的这段频域三维远场外推的程序对不对 Lnl(2xD  
 ein0=cmplx(e0r,e0i)      原点处的入射电场 @l5"nBs<_:  
 hin0=ein0/cexp(cmplx(0,0.5*del_t*2*pi*f))    入射电场滞后半个时间步 (UD@q>c  
 do i=imin,imax 8U"v6S~A%Q  
 do j=jmin,jmax vzAaxk%  
 do k=kmin,kmax :uS\3toj  
 exf(i,j,k)=cmplx(exr(i,j,k),exi(i,j,k))/ein0           电场归一 oG?Xk%7&\  
 eyf(i,j,k)=cmplx(eyr(i,j,k),eyi(i,j,k))/ein0 l}|%5.5-  
 ezf(i,j,k)=cmplx(ezr(i,j,k),ezi(i,j,k))/ein0 |kg7LP3(8,  
 hxf(i,j,k)=cmplx(hxr(i,j,k),hxi(i,j,k))/hin0           磁场归一 Y;M|D'y+  
 hyf(i,j,k)=cmplx(hyr(i,j,k),hyi(i,j,k))/hin0 r.&Vw|*>  
 hzf(i,j,k)=cmplx(hzr(i,j,k),hzi(i,j,k))/hin0 !;v|'I  
 end do ?pmHFlx  
 end do [-K&R  
 end do ZgcMv,=  
       Hyl%mJ  
`UyG_;  
c     外推边界的等效电磁流 9d659iC  
c     左平面 e\l7Iu  
 do j=jomin,jomax-1 Ykw*&opz  
 do k=komin,komax-1 >Eto( y"q  
 jz(iomin,j,k)=-0.25*(hyf(iomin,j,k)+hyf(iomin,j+1,k) 1<@W6@]  
     &     +hyf(iomin-1,j,k)+hyf(iomin-1,j+1,k)) :6 R\OeH+  
      jy(iomin,j,k)=0.25*(hzf(iomin,j,k)+hzf(iomin,j,k+1) O"9\
5(
w  
     &     +hzf(iomin-1,j,k)+hzf(iomin-1,j,k+1)) s WvBv  
 jmz(iomin,j,k)=0.5*(eyf(iomin,j,k)+eyf(iomin,j,k+1)) zT]8KA   
 jmy(iomin,j,k)=-0.5*(ezf(iomin,j,k)+eyf(iomin,j+1,k)) s?}e^/"v  
 kr(iomin,j,k)=iomin*sin(thida)*cos(phi) H[$"+&q  
     &  +(j+0.5)*sin(thida)*sin(phi)+(k+0.5)*cos(phi) ;7V%#-  
 end do R4cM%l_#W  
 end do Y\k#*\'Y~  
  prF%.(G2)  
c 右平面 ^A/k)x6  
      do j=jomin,jomax-1 I-*S&SiXjI  
 do k=komin,komax-1 n0 {i&[I~+  
 jz(iomax,j,k)=0.25*(hyf(iomax,j,k)+hyf(iomax,j+1,k) %)W2H^  
     &     +hyf(iomax-1,j,k)+hyf(iomax-1,j+1,k)) } 9Eg=%0v  
      jy(iomax,j,k)=-0.25*(hzf(iomax,j,k)+hzf(iomax,j,k+1)  skViMo  
     &     +hzf(iomax-1,j,k)+hzf(iomax-1,j,k+1)) 19)i*\+  
 jmz(iomax,j,k)=-0.5*(eyf(iomax,j,k)+eyf(iomax,j,k+1)) L|xbR#v  
 jmy(iomax,j,k)=0.5*(ezf(iomax,j,k)+eyf(iomax,j+1,k)) xGg )Y#  
 kr(iomax,j,k)=iomax*sin(thida)*cos(phi) }@+0/W?\.  
     &  +(j+0.5)*sin(thida)*sin(phi)+(k+0.5)*cos(phi) }B+C~@j  
 end do lT?v^\(H  
 end do      ;bib/  
C     前平面 `iNSr?N.  
      do i=iomin,iomax-1 }!r|1$,kL  
 do k=komin,komax-1 %s|Ely)  
 jz(i,jomin,k)=0.25*(hxf(i,jomin,k)+hxf(i+1,jomin,k)  s};{ZAtE  
     &     +hxf(i,jomin-1,k)+hxf(i+1,jomin-1,k)) O
m\vMd@!  
 jx(i,jomin,k)=-0.25*(hzf(i,jomin,k)+hzf(i,jomin,k+1) @o _}g !9=  
     &     +hzf(i,jomin-1,k)+hzf(i,jomin-1,k+1)) hx%v+/  
      jmz(i,jomin,k)=-0.5*(exf(i,jomin,k)+exf(i,jomin,k+1)) Rtl"Ub@HV  
      jmx(i,jomin,k)=0.5*(ezf(i,jomin,k)+ezf(i+1,jomin,k)) =s2*H8]  
 kr(i,jomin,k)=(i+0.5)*sin(thida)*cos(phi) osAd1<EI
C  
     & +jomin*sin(thida)*sin(phi)+(k+0.5)*cos(phi) k$R-#f;  
 end do PiIpnoM  
 end do #OD/$f_  
C     后平面 4F'LBS]=0  
      do i=iomin,iomax-1 $a"Oc   
 do k=komin,komax-1 AjMh,@  
 jz(i,jomax,k)=-0.25*(hxf(i,jomax,k)+hxf(i+1,jomax,k)  )9`qG:b'  
     &     +hxf(i,jomax-1,k)+hxf(i+1,jomax-1,k)) E,U+o $  
 jx(i,jomax,k)=0.25*(hzf(i,jomax,k)+hzf(i,jomax,k+1) 9}<ile7^  
     &     +hzf(i,jomax-1,k)+hzf(i,jomax-1,k+1)) !)0;&
e5  
      jmz(i,jomax,k)=0.5*(exf(i,jomax,k)+exf(i,jomax,k+1)) ;:g@zAV  
      jmx(i,jomax,k)=-0.5*(ezf(i,jomax,k)+ezf(i+1,jomax,k)) OKR "4n:  
 kr(i,jomax,k)=(i+0.5)*sin(thida)*cos(phi) q dBrQC  
     & +jomax*sin(thida)*sin(phi)+(k+0.5)*cos(phi) pJ"qu,w  
 end do .j0$J\:i  
 end do  ]Ie
 0S~  
c     下平面 *zvx$
yJ?  
      do i=iomin,iomax-1 l'.VKh\C  
 do j=jomin,jomax-1 :q%M_  
 jx(i,j,komin)=0.25*(hyf(i,j,komin)+hyf(i,j+1,komin) WlC:l  
     &     +hyf(i,j,komin-1)+hyf(i,j+1,komin-1)) ]MitOkX  
 jy(i,j,komin)=-0.25*(hxf(i,j,komin)+hxf(i+1,j,komin) g7`LEF <A  
     &     +hxf(i,j,komin-1)+hxf(i+1,j,komin-1)) _op}1   
      jmx(i,j,komin)=-0.5*(eyf(i,j,komin)+eyf(i+1,j,komin))   K`zdc`/  
 jmy(i,j,komin)=0.5*(exf(i,j,komin)+exf(i,j+1,komin)) 6Y?|w3f   
 kr(i,j,komin)=(i+0.5)*sin(thida)*cos(phi) Hk.TM2{w  
     & +(j+0.5)*sin(thida)*sin(phi)+komin*cos(phi) ~KX/ Ai  
 end do )U:m:cr<  
 end do h2]P]@nW;W  
c     上平面 | rtD.,m   
      do i=iomin,iomax-1 {@{']Y  
 do j=jomin,jomax-1 dOH&  
 jx(i,j,komax)=-0.25*(hyf(i,j,komax)+hyf(i,j+1,komax) agDM~=#F  
     &     +hyf(i,j,komax-1)+hyf(i,j+1,komax-1)) :>f )g  
 jy(i,j,komax)=0.25*(hxf(i,j,komax)+hxf(i+1,j,komax) :KP@RZm  
     &     +hxf(i,j,komax-1)+hxf(i+1,j,komax-1)) {qJ1ko)$  
      jmx(i,j,komax)=0.5*(eyf(i,j,komax)+eyf(i+1,j,komax)) #a,PZDaE  
 jmy(i,j,komax)=-0.5*(exf(i,j,komax)+exf(i,j+1,komax)) BG]#o|KW  
 kr(i,j,komax)=(i+0.5)*sin(thida)*cos(phi) 9-a0:bP  
     & +(j+0.5)*sin(thida)*sin(phi)+komax*cos(phi) f+)L#>Gl?  
 end do nT$SfGFj8  
 end do :bq8N@P/  
      ~-Qw.EdC  
C 由于非外推边界面上的电磁流取为0，所以为了方便，取区域内所有电磁流的和 s)D;a-F  
 do i=iomin,iomax ,m|h<faZL
 
 do j=jomin,jomax
mL{6L?  
 do k=komin,komax KBc1{adDx@  
 fx=fx+jx(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))        3个方向的电流矩，KK为波矢 fxHH;hRfv  
      fy=fy+jy(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk)) (v
JNHY M  
      fz=fz+jz(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk)) O-hAFKx  
      fmx=fmx+jmx(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))       3个方向的磁流矩 G}raA%  
      fmy=fmy+jmy(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk)) Vv=. -&'  
      fmz=fmz+jmz(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk)) i3mcx)d@H  
 end do PB*&aYLU  
 end do 8dIgjQX|  
 end do 0znR0%~  
    .g<DD)`  
  m1=z0*(fx*cos(thida)*cos(phi)+fy*cos(thida)*sin(phi)          THIDA方向远场公式的一部分 qt"m  
     &    -fz*sin(thida))+(-fmx*sin(phi)+fmy*cos(phi)) .|fHy  
 t1=(abs(m1))**2     
A]oV"`f  
                 16(QR-  
 m2=z0*(fx*sin(phi)-fy .. p6Gy,C.  
"`e{/7I  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

还是想办法找个已知的例子验证吧。 'e'cb>GnA  
裸看code，有点痛苦。