902阅读
8回复

[求助]谁能帮我看看这个程序？FDTD matlab

UID ：83117

注册：2011-09-22
登录：2012-06-12
发帖：60
等级：仿真一级

0楼发表于: 2012-05-28 13:47:08

哪位高手帮我看看这个程序有什么错误，为什么出不了图形？请各位高手帮忙啊！是关于同轴电缆的FDTD的~~~大家帮帮忙~~~~~感激不尽

clear;
clc;
close all;
%***********************************************************************
%    基本常数
%***********************************************************************
pi = 4*atan(1);
cc=2.99792458e8;            %真空中的波速
mu0=4.0*pi*1.0e-7;          %真空中的磁导系数
eps0=1.0/(cc*cc*mu0);       %真空中的介电系数
Fr=500000000;               %光波频率500MHz
Lamda=cc/Fr;                %波长
k0=2*pi/Lamda;              %波数
bat0=k0*sqrt(eps0);         %传播常数
v=2*pi*Fr/bat0;             %波速
%***********************************************************************
%     同轴电缆参数，参考王均宏论文上的模型
%***********************************************************************
b=8;                        %内导体半径
a=20;                       %外导体半径
c=22.5;                     %护套(绝缘）厚度2.5mm
L=100;                      %电缆长度
%***********************************************************************
%     网格参数
%***********************************************************************
dr=0.5;                    %r方向上的栅距
dfi=10;                    %φ方向上的栅距
dz=0.5;                    %z方向上的栅距
dt=0.5/(2*v);              %时间步
ds=v*dt;                   %空间步
tmax=100;
t0=4;
Nrr=floor((a-b)/dr);        %就近取整，r方向上内导体的外边界到外导体的内边界所占网格数
Nr=floor((c-b)/dr);         %r方向上内导体的外边界到护套外层所占网格数
Nfi=2*pi/dfi;                %φ方向上的网格数
Nz=floor(L/dz);             %z方向上的网格数
Nrb=Nr+1;
Nfib=Nfi+1;
Nzb=Nz+1;
%***********************************************************************
%     求相对介质系数epsr和r(i)
%***********************************************************************
for i=1:1:Nrr             %内导体和外导体之间的介电常数为1.26
    for j=1:1:Nfi
        for k=1:1:Nz
            epsr=1.26*eps0;
        end
    end
end
for i=Nrr+1:1:Nr
    for j=1:1:Nfi          %绝缘层的介电常数为2.3
        for k=1:1:Nz
            epsr=2.3*eps0;
        end
    end
for i=1:Nrb
    r(i,:,:)=dr*i+a; %%%%%%%??????????????????????????????????????????????????????????????????????????
end
%***********************************************************************
%     入射场初始化
%***********************************************************************
Er=zeros(Nr,Nfib,Nzb);
Ef=zeros(Nrb,Nfi,Nzb);
Ez=zeros(Nrb,Nfib,Nz);
Hr=zeros(Nrb,Nfi,Nz);
Hf=zeros(Nr,Nfib,Nz);
Hz=zeros(Nr,Nfi,Nzb);
%***********************************************************************
%     FDTD迭代
%***********************************************************************
for Nt=1:1:tmax

%%%%%%%%%%%%迭代系数%%%%%%%%%%%%
a0=dt/mu0;
a1=a0/dr;
a2=a0/dr/dfi;
a3=a0/dz;
b0=dt/eps;
b1=b0/dr;
b2=b0/dr/dfi;
b3=b0/dz;
%%%%%%第一步求H%%%%%%%%%%%
Hr(2:Nr,1:Nfi,1:Nz)=Hr(2:Nr,1:Nfi,1:Nz)-...
                    a2/r(2:Nr,1:Nfi,1:Nz)*(Ez(2:Nr,2:Nfib,1:Nz)-Ez(2:Nr,1:Nfi,1:Nz))+...
                    a3*(Ef(2:Nr,1:Nfi,2:Nzb)-Ef(2:Nr,1:Nfi,1:Nz));
Hf(1:Nr,2:Nfi,1:Nz)=Hf(1:Nr,2:Nfi,1:Nz)-...
                    a3*(Er(1:Nr,2:Nfi,1:Nz)-Er(1:Nr,2:Nfi,2:Nzb))+...
                    a1*(Ez(2:Nrb,2:Nfi,1:Nz)-Ez(1:Nr,2:Nfi,1:Nz));
Hz(1:Nr,1:Nfi,2:Nz)=Hz(1:Nr,1:Nfi,2:Nz)-...
                    a1/(r(1:Nr,1:Nfi,2:Nz)+ds/2)*((r(1:Nr,1:Nfi,2:Nz)+ds)*Ef(1:Nr,1:Nfi,2:Nz)-r(1:Nr,1:Nfi,2:Nz)*Ef(2:Nbr,1:Nfi,2:Nz))+...
                    a2/(r(1:Nr,1:Nfi,2:Nz)+ds/2)*(Er(1:Nr,2:Nfib,2:Nz)-Er(1:Nr,1:Nfi,2:Nz));

%%%%%%第一步求E%%%%%%%%%%%
Er(1:Nr,2:Nfb,2:Nz)=Er(1:Nr,2:Nfi,2:Nz)+...
                      b2/(r(1:Nr,1:Nfi,2:Nz)+ds/2)*(Hz(1:Nr,2:Nfi,2:Nz)-Hz(1:Nr,1:Nfi-1,2:Nz))-...
                      b3*(Hf(1:Nr,2:Nfi,2:Nz)-Hf(1:Nr,2:Nfi,1:Nz-1));
Ef(2:Nr,1:Nfi,2:Nz)=Ef(2:Nr,1:Nfi,2:Nz)+...
                     b3*(Hr(2:Nr,1:Nfi,2:Nz)-Hr(2:Nr,1:Nfi,1:Nz-1))-...
                     b1*(Hz(2:Nr,1:Nfi,2:Nz)-Hz(1:Nr-1,1:Nfi,2:Nz));

Ez(2:Nr,2:Nfi,1:Nz)=Ez(2:Nr,2:Nfi,1:Nz)+...
                     b1/r(2:Nr,2:Nfi,1:Nz)*((r(2:Nr,2:Nfi,1:Nz)+ds/2)*Hf(2:Nr,2:Nfi,1:Nz)-(r(2:Nr,2:Nfi,1:Nz)-ds/2)*Hf(1:Nr-1,2:Nfi,1:Nz))-...
                     b2/r(2:Nr,2:Nfi,1:Nz)*(Hr(2:Nr,2:Nfi,1:Nz)-Hr(2:Nr,1:Nfi-1,1:Nz));
%***********************************************************************
%     导体处电场切向为0
%***********************************************************************
Ef(Nr,1:Nfi,1:Nz)=0;
Ez(Nr,1:Nfi,1:Nz)=0;
%***********************************************************************
%     激 ..

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