1587阅读
1回复

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

UID ：41252

注册：2009-09-09
登录：2010-07-20
发帖：49
等级：仿真一级

0楼发表于: 2009-10-13 21:37:21

最近编了一段球的散射的程序，大家帮我看下我的这段频域三维远场外推的程序对不对
ein0=cmplx(e0r,e0i)      原点处的入射电场
hin0=ein0/cexp(cmplx(0,0.5*del_t*2*pi*f))    入射电场滞后半个时间步
do i=imin,imax
do j=jmin,jmax
do k=kmin,kmax
exf(i,j,k)=cmplx(exr(i,j,k),exi(i,j,k))/ein0           电场归一
eyf(i,j,k)=cmplx(eyr(i,j,k),eyi(i,j,k))/ein0
ezf(i,j,k)=cmplx(ezr(i,j,k),ezi(i,j,k))/ein0
hxf(i,j,k)=cmplx(hxr(i,j,k),hxi(i,j,k))/hin0          磁场归一
hyf(i,j,k)=cmplx(hyr(i,j,k),hyi(i,j,k))/hin0
hzf(i,j,k)=cmplx(hzr(i,j,k),hzi(i,j,k))/hin0
end do
end do
end do


c     外推边界的等效电磁流
c     左平面
do j=jomin,jomax-1
do k=komin,komax-1
jz(iomin,j,k)=-0.25*(hyf(iomin,j,k)+hyf(iomin,j+1,k)
     &     +hyf(iomin-1,j,k)+hyf(iomin-1,j+1,k))
      jy(iomin,j,k)=0.25*(hzf(iomin,j,k)+hzf(iomin,j,k+1)
     &     +hzf(iomin-1,j,k)+hzf(iomin-1,j,k+1))
jmz(iomin,j,k)=0.5*(eyf(iomin,j,k)+eyf(iomin,j,k+1))
jmy(iomin,j,k)=-0.5*(ezf(iomin,j,k)+eyf(iomin,j+1,k))
kr(iomin,j,k)=iomin*sin(thida)*cos(phi)
     & +(j+0.5)*sin(thida)*sin(phi)+(k+0.5)*cos(phi)
end do
end do

c 右平面
      do j=jomin,jomax-1
do k=komin,komax-1
jz(iomax,j,k)=0.25*(hyf(iomax,j,k)+hyf(iomax,j+1,k)
     &     +hyf(iomax-1,j,k)+hyf(iomax-1,j+1,k))
      jy(iomax,j,k)=-0.25*(hzf(iomax,j,k)+hzf(iomax,j,k+1)
     &     +hzf(iomax-1,j,k)+hzf(iomax-1,j,k+1))
jmz(iomax,j,k)=-0.5*(eyf(iomax,j,k)+eyf(iomax,j,k+1))
jmy(iomax,j,k)=0.5*(ezf(iomax,j,k)+eyf(iomax,j+1,k))
kr(iomax,j,k)=iomax*sin(thida)*cos(phi)
     & +(j+0.5)*sin(thida)*sin(phi)+(k+0.5)*cos(phi)
end do
end do
C     前平面
      do i=iomin,iomax-1
do k=komin,komax-1
jz(i,jomin,k)=0.25*(hxf(i,jomin,k)+hxf(i+1,jomin,k)
     &     +hxf(i,jomin-1,k)+hxf(i+1,jomin-1,k))
jx(i,jomin,k)=-0.25*(hzf(i,jomin,k)+hzf(i,jomin,k+1)
     &     +hzf(i,jomin-1,k)+hzf(i,jomin-1,k+1))
      jmz(i,jomin,k)=-0.5*(exf(i,jomin,k)+exf(i,jomin,k+1))
      jmx(i,jomin,k)=0.5*(ezf(i,jomin,k)+ezf(i+1,jomin,k))
kr(i,jomin,k)=(i+0.5)*sin(thida)*cos(phi)
     & +jomin*sin(thida)*sin(phi)+(k+0.5)*cos(phi)
end do
end do
C     后平面
      do i=iomin,iomax-1
do k=komin,komax-1
jz(i,jomax,k)=-0.25*(hxf(i,jomax,k)+hxf(i+1,jomax,k)
     &     +hxf(i,jomax-1,k)+hxf(i+1,jomax-1,k))
jx(i,jomax,k)=0.25*(hzf(i,jomax,k)+hzf(i,jomax,k+1)
     &     +hzf(i,jomax-1,k)+hzf(i,jomax-1,k+1))
      jmz(i,jomax,k)=0.5*(exf(i,jomax,k)+exf(i,jomax,k+1))
      jmx(i,jomax,k)=-0.5*(ezf(i,jomax,k)+ezf(i+1,jomax,k))
kr(i,jomax,k)=(i+0.5)*sin(thida)*cos(phi)
     & +jomax*sin(thida)*sin(phi)+(k+0.5)*cos(phi)
end do
end do
c     下平面
      do i=iomin,iomax-1
do j=jomin,jomax-1
jx(i,j,komin)=0.25*(hyf(i,j,komin)+hyf(i,j+1,komin)
     &     +hyf(i,j,komin-1)+hyf(i,j+1,komin-1))
jy(i,j,komin)=-0.25*(hxf(i,j,komin)+hxf(i+1,j,komin)
     &     +hxf(i,j,komin-1)+hxf(i+1,j,komin-1))
      jmx(i,j,komin)=-0.5*(eyf(i,j,komin)+eyf(i+1,j,komin))
jmy(i,j,komin)=0.5*(exf(i,j,komin)+exf(i,j+1,komin))
kr(i,j,komin)=(i+0.5)*sin(thida)*cos(phi)
     & +(j+0.5)*sin(thida)*sin(phi)+komin*cos(phi)
end do
end do
c     上平面
      do i=iomin,iomax-1
do j=jomin,jomax-1
jx(i,j,komax)=-0.25*(hyf(i,j,komax)+hyf(i,j+1,komax)
     &     +hyf(i,j,komax-1)+hyf(i,j+1,komax-1))
jy(i,j,komax)=0.25*(hxf(i,j,komax)+hxf(i+1,j,komax)
     &     +hxf(i,j,komax-1)+hxf(i+1,j,komax-1))
      jmx(i,j,komax)=0.5*(eyf(i,j,komax)+eyf(i+1,j,komax))
jmy(i,j,komax)=-0.5*(exf(i,j,komax)+exf(i,j+1,komax))
kr(i,j,komax)=(i+0.5)*sin(thida)*cos(phi)
     & +(j+0.5)*sin(thida)*sin(phi)+komax*cos(phi)
end do
end do

C 由于非外推边界面上的电磁流取为0，所以为了方便，取区域内所有电磁流的和
do i=iomin,iomax
do j=jomin,jomax
do k=komin,komax
fx=fx+jx(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))        3个方向的电流矩，KK为波矢
      fy=fy+jy(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))
      fz=fz+jz(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))
      fmx=fmx+jmx(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))       3个方向的磁流矩
      fmy=fmy+jmy(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))
      fmz=fmz+jmz(i,j,k)*exp(cmplx(0,kr(i,j,k)*kk))
end do
end do
end do

  m1=z0*(fx*cos(thida)*cos(phi)+fy*cos(thida)*sin(phi)          THIDA方向远场公式的一部分
     &   -fz*sin(thida))+(-fmx*sin(phi)+fmy*cos(phi))
t1=(abs(m1))**2

m2=z0*(fx*sin(phi)-fy ..

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

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

1楼发表于: 2009-10-15 14:07:34

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

逆流而上

发帖回复