1970阅读
6回复

[求助]RCS计算程序

UID ：11863

注册：2008-05-03
登录：2021-10-24
发帖：35
等级：仿真新人

0楼发表于: 2009-12-21 21:36:20

编了好久的RCS的程序，结果总是有问题，有没有高手看下程序有什么问题，编的程序是按照葛老书附录写的，如下：
clear all;
Media=zeros(4,15);%各介质的电磁参数
cc=2.99792458e8;
muz=4.0*pi*1.0e-7;
mur=1;
epsz=1.0/(cc*cc*muz);
freq=3.0e10;
WaveLength=cc/freq;
WL=40;%数值波长，每波长的网格数
i=0;
j=0;
Iflag=1;
IflagI=1;
dx=WaveLength/WL;
dt=WaveLength/(WL*2*3e8);
ie=100;
je=100;
ib=ie+1;
jb=je+1;
%定义电磁场变量
Ez=zeros(ib,jb);%FDTD采样点Ez
Hx=zeros(ib,je);%FDTD采样点Hx
Hy=zeros(ie,jb);%FDTD采样点Hy

Ein=zeros(1600,1);%入射波E采样点
Hin=zeros(1600,1);%入射波H采样点

EBin=zeros(4);%入射波吸收边界条件临时变量

ECB=zeros(ib,4);%连接边界处E入射波的分量
HCB=zeros(ib,4);%连接边界处H入射波的分量

EAB=zeros(ib,4);%吸收边界条件临时变量
EAC=zeros(4,2,2);%吸收边界条件角点临时变量

FE1=zeros(11,11);%FDTD迭代系数
FE2=zeros(11,11);
FH1=zeros(11,11);
FH2=zeros(11,11);

Eini=0;

TEMflag=1;%1代表TM，2代表TE

Mx=0;
My=0;
Jz=0;
Step=0;

T1=0;%存储临时变量
T2=0;

CT3=0;
CT4=0;
I=0;
J=0;
N=0;
II=0;
III=0;
MediaNumber=2;%目标所包含介质总数目
rcs=zeros(181,1);
%%%%%%%%%%%%%%%%%%初始化变量
TEMFlag=1;

Isource=3;

Imin=1;
Imax=80;
Jmin=1;
Jmax=80;
Itmin=10;%连接边界
Itmax=50;
Jtmin=10;
Jtmax=50;
iomax=53;%输出边界，用于近远场的变换
iomin=7;
jomax=53;
jomin=7;
length=iomax-iomin+1;
outdataez=zeros(4,length);
outdatahx=zeros(4,length);
outdatahy=zeros(4,length);
outdataezr=zeros(4,length);
outdatahxr=zeros(4,length);
outdatahyr=zeros(4,length);
ob=ones(Imax,Jmax);%存放网格介质的编号
obeps=ones(Imax,Jmax);
Phi=0;
Phi=Phi*pi/180;
cp=cos(Phi);
sp=sin(Phi);
TimeStop=1200;
IncidentStart=1;
IncidentEnd=1420;

Isbias=Isource;
%%%%%%%%%%%%%%常量初始化
MU0=pi*0.0000004;
EPS0=8.85e-12;
SMU0=sqrt(MU0);
SEPS0=sqrt(EPS0);
Z=sqrt(MU0/EPS0);
k=2*pi/WaveLength;

Media(1,1)=1;%背景为真空epsilon_r
Media(2,1)=1;%mu_r
Media(3,1)=0;%导电率
Media(4,1)=0;%磁导率

MediaNo=1;%只有一种介质
Media(1,2)=1;%epsilon_r
Media(2,2)=1;%mu_r
Media(3,2)=1e7;%导电率
% Media(3,2)=0;%导电率
Media(4,2)=0;%磁导率

for i=Imin:Imax
    for j=Jmin:Jmax
        if ((i-(Imin+Imax)/2)^2+(j-(Jmin+Jmax)/2)^2)<100
            ob(i,j)=2;
        end
    end
end

for i=1:MediaNumber
    Media(3,i)=-Media(3,i)*Z/k;
    Media(4,i)=-Media(4,i)/(Z*k);
end

if TEMFlag==1 %计算FDTD迭代公式系数
    for i=1:MediaNumber
        for j=1:MediaNumber
            FE1(i,j)=((Media(1,i)+Media(1,j))+pi*0.5*(Media(3,i)+Media(3,j))/WL)/((Media(1,i)+Media(1,j))-pi*0.5*(Media(3,i)+Media(3,j))/WL);
            FE2(i,j)=1./((Media(1,i)+Media(1,j))-pi*0.5*(Media(3,i)+Media(3,j))/WL);
            FH1(i,j)=((Media(2,i)+Media(2,j))+pi*0.5*(Media(4,i)+Media(4,j))/WL)/((Media(2,i)+Media(2,j))-pi*0.5*(Media(4,i)+Media(4,j))/WL);
            FH2(i,j)=1./((Media(2,i)+Media(2,j))-pi*0.5*(Media(4,i)+Media(4,j))/WL);
        end
    end
else
    for i=1:MediaNumber
        for j=1:MediaNumber
            FH1(i,j)=((Media(1,i)+Media(1,j))+pi*0.5*(Media(3,i)+Media(3,j))/WL)/((Media(1,i)+Media(1,j))-pi*0.5*(Media(3,i)+Media(3,j))/WL);
            FH2(i,j)=-1./((Media(1,i)+Media(1,j))-pi*0.5*(Media(3,i)+Media(3,j))/WL);
            FE1(i,j)=((Media(2,i)+Media(2,j))+pi*0.5*(Media(4,i)+Media(4,j))/WL)/((Media(2,i)+Media(2,j))-pi*0.5*(Media(4,i)+Media(4,j))/WL);
            FE2(i,j)=-1./((Media(2,i)+Media(2,j))-pi*0.5*(Media(4,i)+Media(4,j))/WL);
        end
    end
end
FH=FH2(1,1);
FE=FE2(1,1);
OMIGA=pi/WL;

% subplot(3,1,1),pcolor(Hx');
% shading flat;
% caxis([-80.0 80.0]);
% axis([1 60 1 60]);
% colorbar;
% axis image;
% axis off;
% title('Hx at time step = 0');
%
% subplot(3,1,2),pcolor(Hy');
% shading flat;
% caxis([-80.0 80.0]);
% axis([1 60 1 60]);
% colorbar;
% axis image;
% axis off;
% title('Hy at time step = 0');
%
% subplot(3,1,3),pcolor(Ez');
% shading flat;
% caxis([-0.2 0.2]);
% axis([1 60 1 60]);
% colorbar;
% axis image;
% axis off;
% title('Ez at time step = 0');
%
% rect=get(gcf,'Position');
% rect(1:2)=[0 0];
%
%
%
% M=moviein(TimeStop/4,gcf,rect);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          主程序
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for N=1:TimeStop+WL/2
    %%%%%%%%%%产生入射波
    for i=IncidentStart:IncidentEnd-1
        Hin(i)=Hin(i)-FH*(Ein(i+1)-Ein(i));
    end
    for i=IncidentStart+1:IncidentEnd-1
        Ein(i)=Ein(i)-FE*(Hin(i)-Hin(i-1));
    end
    %%%加入源
    Ein(Isource)=sin(OMIGA*N);
    if N<WL
        Ein(Isource)=Ein(Isource)*0.5*(1-cos(OMIGA*N));%开关函数，升余弦函数
    end
    Ein(IncidentStart)=EBin(4);
    EBin(4)=EBin(3);
    EBin(3)=Ein(IncidentStart+1);
    Ein(IncidentEnd)=EBin(2);
    EBin(2)=EBin(1);
    EBin(1)=Ein(IncidentEnd-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       连接边界处入射场各分量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=Itmin:Itmax
    T1=I*cp+Jtmin*sp+Isbias;
    II=round(T1);
    III=II+1;
    T1=III-T1;
    ECB(I,1)=T1*(Ein(II)-Ein(III))+Ein(III);%%%%%%%%%底部连接边界处入射电场分量

    T1=I*cp+Jtmax*sp+Isbias;
    II=round(T1);
    III=II+1;
    T1=III-T1;
    ECB(I,3)=T1*(Ein(II)-Ein(III))+Ein(III);%%%%%%%%%%%%%顶部连接边界处入射电场分量
end
for J=Jtmin:Jtmax
    T1=J*sp+Itmin*cp+Isbias;
    II=round(T1);
    III=II+1;
    T1=III-T1;
    ECB(J,4)=T1*(Ein(II)-Ein(III))+Ein(III);%%%%%%%%%%%%%左边连接边界处入射电场分量

    T1=J*sp+Itmax*cp+Isbias;
    II=round(T1);
    III=II+1;
    T1=III-T1;
    ECB(J,2)=T1*(Ein(II)-Ein(III))+Ein(III);%%%%%%%%%%%%%右边连接边界处入射电场分量
end

for I=Itmin:Itmax
    T1=I*cp+(Jtmin-0.5)*sp+Isbias;
    II=round(T1-0.5);
    III=II+1;
    T1=III+0.5-T1;
    HCB(I,1)=T1*(Hin(II)-Hin(III))+Hin(III);
    HCB(I,1)=HCB(I,1)*sp;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%底部连接边界入射波磁场分量Hx

    T1=I*cp+(Jtmax+0.5)*sp+Isbias;
    II=round(T1-0.5);
    III=II+1;
    T1=III+0.5-T1;
    HCB(I,3)=T1*(Hin(II)-Hin(III))+Hin(III);
    HCB(I,3)=HCB(I,3)*sp;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%顶部连接边界入射波磁场分量Hx
end

for J=Jtmin:Jtmax
    T1=J*sp+(Itmin-0.5)*cp+Isbias;
    II=round(T1-0.5);
    III=II+1;
    T1=III+0.5-T1;
    HCB(J,4)=T1*(Hin(II)-Hin(III))+Hin(III);
    HCB(J,4)=-HCB(J,4)*cp;%%%%%%%%%%%%%%%%%%%%左边连接边界入射波磁场分量

    T1=J*sp+(Itmax+0.5)*cp+Isbias;
    II=round(T1-0.5);
    III=II+1;
    T1=III+0.5-T1;
    HCB(J,2)=T1*(Hin(II)-Hin(III))+Hin(III);
    HCB(J,2)=-HCB(J,2)*cp;%%%%%%%%%%%%%%%%%%%%左边连接边界入射波磁场分量

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     Ez的FDTD迭代公式
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=Imin+1:Imax-1
    for J=Jmin+1:Jmax-1
        T1=Hy(I,J)-Hy(I-1,J)-Hx(I,J)+Hx(I,J-1);
        Ez(I,J)=FE1(ob(I,J),ob(I,J))*Ez(I,J)+FE2(ob(I,J),ob(I,J))*T1;
        %Dz(I,J)=Dz(I,J)+(1/(2))*T1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     连接边界处引入入射波Ez分量         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=Itmin+1:Itmax-1
    Ez(I,Jtmax)=Ez(I,Jtmax)-FE*HCB(I,3);
    Ez(I,Jtmin)=Ez(I,Jtmin)+FE*HCB(I,1);
end

for J=Jtmin+1:Jtmax-1
    Ez(Itmax,J)=Ez(Itmax,J)+FE*HCB(J,2);
    Ez(Itmin,J)=Ez(Itmin,J)-FE*HCB(J,4);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         连接边界角点引入入射波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T1=HCB(Jtmax,2)-HCB(Itmax,3);
Ez(Itmax,Jtmax)=Ez(Itmax,Jtmax)+FE*T1;

T1=-HCB(Jtmin,4)+HCB(Itmin,1);
Ez(Itmin,Jtmin)=Ez(Itmin,Jtmin)+FE*T1;

T1=HCB(Jtmin,2)+HCB(Itmax,1);
Ez(Itmax,Jtmin)=Ez(Itmax,Jtmin)+FE*T1;

T1=-HCB(Itmin,3)-HCB(Jtmax,4);
Ez(Itmin,Jtmax)=Ez(Itmin,Jtmax)+FE*T1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         吸收边界条件Ez
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=Imin+1:Imax-1
    T1=Hy(I,Jmin)-Hy(I-1,Jmin);
    T1=T1+Hy(I,Jmin+1)-Hy(I-1,Jmin+1);
    T1=T1*TEMFlag;
    T2=Ez(I,Jmin+1)-Ez(I,Jmin);
    Ez(I,Jmin)=EAB(I,1)-1/3*T2;
    Ez(I,Jmin)=Ez(I,Jmin)+1/6*T1;
    EAB(I,1)=Ez(I,Jmin+1);

    T1=Hy(I,Jmax)-Hy(I-1,Jmax);
    T1=T1+Hy(I,Jmax-1)-Hy(I-1,Jmax-1);
    T1=T1*TEMFlag;
    T2=Ez(I,Jmax-1)-Ez(I,Jmax);
    Ez(I,Jmax)=EAB(I,3)-1/3*T2;
    Ez(I,Jmax)=Ez(I,Jmax)+1/6*T1;
    EAB(I,3)=Ez(I,Jmax-1);
end

for J=Jmin+1:Jmax-1
    T1=Hx(Imin,J)-Hx(Imin,J-1);
    T1=T1+Hx(Imin+1,J)-Hx(Imin+1,J-1);
    T1=T1*TEMFlag;
    T2=Ez(Imin+1,J)-Ez(Imin,J);
    Ez(Imin,J)=EAB(J,4)-1/3*T2;
    Ez(Imin,J)=Ez(Imin,J)-1/6*T1;
    EAB(J,4)=Ez(Imin+1,J);

    T1=Hx(Imax,J)-Hx(Imax,J-1);
    T1=T1+Hx(Imax-1,J)-Hx(Imax-1,J-1);
    T1=T1*TEMFlag;
    T2=Ez(Imax-1,J)-Ez(Imax,J);
    Ez(Imax,J)=EAB(J,2)-1/3*T2;
    Ez(Imax,J)=Ez(Imax,J)-1/6*T1;
    EAB(J,2)=Ez(Imax-1,J);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         吸收边界角点Ez
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
IflagI=1-IflagI;%%交换存储标志
Iflag=IflagI+1;
Ez(Imin,Jmin)=0.2928932*EAC(1,1,Iflag)+0.7071068*EAC(1,2,Iflag);
EAC(1,1,Iflag)=Ez(Imin,Jmin);
EAC(1,2,Iflag)=Ez(Imin+1,Jmin+1);

Ez(Imax,Jmin)=.2928932*EAC(2,1,Iflag)+0.7071068*EAC(2,2,Iflag);
EAC(2,1,Iflag)=Ez(Imax,Jmin);
EAC(2,2,Iflag)=Ez(Imax-1,Jmin+1);

Ez(Imin,Jmax)=0.2928932*EAC(3,1,Iflag)+0.7071068*EAC(3,2,Iflag);
EAC(3,1,Iflag)=Ez(Imin,Jmax);
EAC(3,2,Iflag)=Ez(Imin+1,Jmax-1);

Ez(Imax,Jmax)=0.2928932*EAC(4,1,Iflag)+0.7071068*EAC(4,2,Iflag);
EAC(4,1,Iflag)=Ez(Imax,Jmax);
EAC(4,2,Iflag)=Ez(Imax-1,Jmax-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      Hx Hy 的FDTD迭代公式
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=Imin:Imax
    for J=Jmin:Jmax-1
        T1=Ez(I,J+1)-Ez(I,J);

        Hx(I,J)=FH1(ob(I,J),ob(I,J+1))*Hx(I,J)-FH2(ob(I,J),ob(I,J+1))*T1;
    end
end

for I=Imin:Imax-1
    for J=Jmin:Jmax
        T1=Ez(I+1,J)-Ez(I,J);
        Hy(I,J)=FH1(ob(I,J),ob(I+1,J))*Hy(I,J)+FH2(ob(I,J),ob(I+1,J))*T1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         连接边界引入入射波Hx Hy分量
for I=Itmin:Itmax
    Hx(I,Jtmin-1)=Hx(I,Jtmin-1)+FH*ECB(I,1);
    Hx(I,Jtmax)=Hx(I,Jtmax)-FH*ECB(I,3);
end

for J=Jtmin:Jtmax
    Hy(Itmin-1,J)=Hy(Itmin-1,J)-FH*ECB(J,4);
    Hy(Itmax,J)=Hy(Itmax,J)+FH*ECB(J,2);
end

%输出边界上的虚部
if N==TimeStop
   %output Ez imaginary part
   for i=iomin:iomax
       outdataez(1,i)=0.5*(Ez(i,jomin)+Ez(i,jomin-1));%down
       outdatahx(1,i)=Hx(i,jomin-1)/Z;                %down
       outdataez(2,i)=0.5*(Ez(i,jomax)+Ez(i,jomax+1));%up
       outdatahx(2,i)=Hx(i,jomax)/Z;                  %up
   end

   for j=jomin:jomax
       outdataez(3,j)=0.5*(Ez(iomin,j)+Ez(iomin-1,j));%left
       outdatahy(1,j)=Hy(iomin-1,j)/Z;                %left
       outdataez(4,j)=0.5*(Ez(iomax,j)+Ez(iomax+1,j));%right
       outdatahy(2,j)=Hy(iomax,j)/Z;                  %right
   end
   Eini=Ein(Isource);
end

% if mod(N,4)==0;
%
% timestep=int2str(N);
%
% subplot(3,1,1),pcolor(Hx');
% shading flat;
% caxis([-1.3 1.3]);
% axis([1 60 1 60]);
% colorbar;
% axis image;
% axis off;
% title(['Hx at time step = ',timestep]);
%
% subplot(3,1,2),pcolor(Hy');
% shading flat;
% caxis([-1.3 1.3]);
% axis([1 60 1 60]);
% colorbar;
% axis image;
% axis off;
% title(['Hy at time step = ',timestep]);
%
% subplot(3,1,3),pcolor(Ez');
% shading flat;
% caxis([-1.6 1.6]);
% axis([1 60 1 60]);
% colorbar;
% axis image;
% axis off;
% title(['Ez at time step = ',timestep]);
%
% nn=N/4;
% M(:,nn)=getframe(gcf,rect);
%
% end;

end
Einr=Ein(Isource);

%输出结果组成复数
%实部
   for i=iomin:iomax
       outdataezr(1,i)=0.5*(Ez(i,jomin)+Ez(i,jomin-1));%down
       outdatahxr(1,i)=Hx(i,jomin-1)/Z;                %down
       outdataezr(2,i)=0.5*(Ez(i,jomax)+Ez(i,jomax+1));%up
       outdatahxr(2,i)=Hx(i,jomax)/Z;                  %up
   end

   for j=jomin:jomax
       outdataezr(3,j)=0.5*(Ez(iomin,j)+Ez(iomin-1,j));%left
       outdatahyr(1,j)=Hy(iomin-1,j)/Z;                %left
       outdataezr(4,j)=0.5*(Ez(iomax,j)+Ez(iomax+1,j));%right
       outdatahyr(2,j)=Hy(iomax,j)/Z;                  %right
   end
CT1=1/(Einr+sqrt(-1)*Eini);%入射波电场在原点处的值的倒数
CT2=exp(-0.5*sqrt(-1)*OMIGA)*CT1;%CT2比CT1倒数滞后半个时间步
%对电场和磁场归一化
   for i=iomin:iomax
       outdataez(1,i)=CT1*(outdataezr(1,i)+sqrt(-1)*outdataez(1,i));%down
       outdatahx(1,i)=CT2*(outdatahxr(1,i)+sqrt(-1)*outdatahx(1,i));%down
       outdataez(2,i)=CT1*(outdataezr(2,i)+sqrt(-1)*outdataez(2,i));%up
       outdatahx(2,i)=CT2*(outdatahxr(2,i)+sqrt(-1)*outdatahx(2,i));%up
   end

   for j=jomin:jomax ..

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

离线casey

有失才有得

UID ：2188

注册：2007-04-25
登录：2013-10-17
发帖：1813
等级：荣誉管理员

1楼发表于: 2009-12-22 12:40:52

看别人code，最头疼了.....

Casey

离线kerbcurb

UID ：5263

注册：2007-10-02
登录：2024-02-10
发帖：417
等级：积极交流五级

2楼发表于: 2009-12-24 10:50:37

你应该把问题说清楚，从程序看，好像2D方柱的散射吧？似乎是后面的部分有问题。
我感觉这一句：CT1=sqrt(-1)*(i*cos(phi)+(jomin-0.5)*sin(phi))*rk; 似乎有问题，
没说错的话，你是想用这个这个表示入射波，我记得入射波一般表示成：x *cos(phi) + y * sin(phi)，x, y的取值范围是边宽，但是这里 i 的取值范围是iomin:iomax，这里是不是有问题？

代码太长，前面的代码没有仔细看，有没有错误不知道。

不知道说的对不对，很少做FDTD

共1条评分

cc81372365

rf币 +1

积极参与论坛交流，欢迎继续参与本贴交流！

2010-12-28

离线xiaobestlc

UID ：38453

注册：2009-07-31
登录：2011-03-20
发帖：29
等级：仿真一级

3楼发表于: 2010-01-14 18:35:48

2D方柱我用C做了个，楼主需要可以联系一起讨论

离线沔水长清

UID ：94167

注册：2012-05-17
登录：2017-03-13
发帖：126
等级：仿真二级

4楼发表于: 2012-07-13 11:10:10

好长啊。。。。。。。。。

离线smlqfouk

UID ：95065

注册：2012-06-03
登录：2019-03-04
发帖：29
等级：仿真新人

5楼发表于: 2012-07-19 12:01:59

值得收藏...

張曉舟游本昌楊海彪 avast-nod32-antivirus.org/

离线qazwsx

UID ：18226

注册：2008-09-28
登录：2017-02-16
发帖：13
等级：仿真新人

6楼发表于: 2017-02-15 14:47:45

这和································

发帖回复