登录 | 注册
开启辅助访问
设为首页加入收藏关于我们

头雁微网

头雁微网»论坛 › 技术应用 › 情报信息 › 几个用Matlab和C写的关于FDTD程序
1234下一页
返回列表 发新帖
收起左侧
楼主: gaohy2006 - 

几个用Matlab和C写的关于FDTD程序

  [复制链接]
发表于 2008-6-10 15:39:14  | 显示全部楼层
%COXIAL
clear
  
%***********************************************************************
%     Fundamental constants
%***********************************************************************
  
cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
epsz=1.0/(cc*cc*muz);       %permittivity of free space
ii=sqrt(-1)

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=34;      %number of grid cells in x-direction  r=17mm  rout=34mm
je=20;
ke=15;       %number of grid cells in z-direction  l=7.5mm
  
ib=ie+1;
jb=je+1;
kb=ke+1;

dtheta=2*pi/je;
dx=0.0005;          %space increment of cubic lattice 0.5mm
dt=dx/(2.0*cc);     %time step
nmax=2000;          %total number of time steps
cb=dt/epsz;
db=dt/muz/dx;

er=zeros(ie,jb,kb);
et=zeros(ib,jb,kb);
ez=zeros(ib,jb,ke);
hr=zeros(ib,jb,ke);
ht=zeros(ie,jb,ke);
hz=zeros(ie,jb,kb);
r=zeros(ib,jb,kb);
ez0=zeros(ib,jb);
Nf=1000
fmax=100e+9
df=fmax/Nf

f0=30e+9


ndecay=1/(2*fmax*dt);
n0=1/(2*f0*dt);
t=3.123/dx/ie/2;
%compute the r(i)
for i=1:ib
    r(i,:,:)=dx*ie+dx*(i-1);
end
    for i=1:ib
        ez0(i,:)=besselj(0,t*(ie+i-1)*dx)*bessely(0,t*dx*2*ie)-bessely(0,t*(ie+i-1)*dx)*besselj(0,t*2*dx*ie);
    end
   
    %***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************
  
for n=1:nmax
  
%***********************************************************************
%     Update electric fields
%***********************************************************************
er(1:ie,2:jb,2:ke)=er(1:ie,2:jb,2:ke)+...
                    dt/dtheta/epsz/(r(1:ie,2:jb,2:ke)+dx/2).*(hz(1:ie,2:jb,2:ke)-hz(1:ie,1:je,2:ke))-...
                    dt/epsz/dx*(ht(1:ie,2:jb,2:ke)-ht(1:ie,2:jb,1:ke-1));

ez(2:ie,2:jb,1:ke)=ez(2:ie,2:jb,1:ke)+...
                   dt/epsz/dx*(ht(2:ie,2:jb,1:ke)-ht(1:ie-1,2:jb,1:ke))+dt/epsz/2*(ht(2:ie,2:jb,1:ke)+ht(1:ie-1,2:jb,1:ke))./r(2:ie,2:jb,1:ke)-...
                   dt/dtheta/epsz*(hr(2:ie,2:jb,1:ke)-hr(2:ie,1:je,1:ke))./r(2:ie,2:jb,1:ke);

et(2:ie,1:jb,2:ke)=et(2:ie,1:jb,2:ke)+dt/epsz/dx*(hr(2:ie,1:jb,2:ke)-hr(2:ie,1:jb,1:ke-1)-hz(2:ie,1:jb,2:ke)+hz(1:ie-1,1:jb,2:ke));
                  

   %add excite**************************depend upon the mode which compute.
         ez(5,5,3)=exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(5,5,3);
        %ez(:,:,3)=ez0(:,:)*exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(:,:,3);
   %************* periodic  boundary condition   for theta direction *****************************

er(1:ie,1,2:ke)=er(1:ie,1,2:ke)+...
                    dt/dtheta/epsz/(r(1:ie,1,2:ke)+dx/2).*(hz(1:ie,1,2:ke)-hz(1:ie,je,2:ke))-...
                   dt/epsz/dx*(ht(1:ie,1,2:ke)-ht(1:ie,1,1:ke-1));

ez(2:ie,1,1:ke)=ez(2:ie,1,1:ke)+...
                   dt/epsz/dx*(ht(2:ie,1,1:ke)-ht(1:ie-1,1,1:ke))+dt/epsz/2*(ht(2:ie,1,1:ke)+ht(1:ie-1,1,1:ke))./r(2:ie,1,1:ke)-...
                   dt/dtheta/epsz*(hr(2:ie,1,1:ke)-hr(2:ie,je,1:ke))./r(2:ie,1,1:ke);


%***********************************************************************
%     Update magnetic fields
%***********************************************************************
hr(2:ie,1:je,1:ke)=hr(2:ie,1:je,1:ke)-...
                    dt/muz/dtheta*(ez(2:ie,2:jb,1:ke)-ez(2:ie,1:je,1:ke))./r(2:ie,1:je,1:ke)+dt/muz/dx*(et(2:ie,1:je,2:kb)-et(2:ie,1:je,1:ke));
  
ht(1:ie,1:jb,1:ke)=ht(1:ie,1:jb,1:ke)+...
                   dt/muz/dx*(er(1:ie,1:jb,1:ke)-er(1:ie,1:jb,2:kb)+ez(2:ib,1:jb,1:ke)-ez(1:ie,1:jb,1:ke));

hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)-dt/muz/2/(r(1:ie,1:je,2:ke)+dx/2).*(et(1:ie,1:je,2:ke)+et(2:ib,1:je,2:ke))-...
                     dt/muz/dx*(et(2:ib,1:je,2:ke)-et(1:ie,1:je,2:ke))+dt/muz/dtheta/(r(1:ie,1:je,2:ke)+dx/2).*(er(1:ie,2:jb,2:ke)-er(1:ie,1:je,2:ke));  
                     
%************* periodic  boundary condition   for theta direction *****************************
hr(2:ie,jb,1:ke)=hr(2:ie,jb,1:ke)-...
                    dt/muz/dtheta*(ez(2:ie,2,1:ke)-ez(2:ie,jb,1:ke))./r(2:ie,jb,1:ke)+dt/muz/dx*(et(2:ie,jb,2:kb)-et(2:ie,jb,1:ke));
               
hz(1:ie,jb,2:ke)=hz(1:ie,jb,2:ke)-dt/muz/2/(r(1:ie,jb,2:ke)+dx/2).*(et(1:ie,jb,2:ke)+et(2:ib,jb,2:ke))-...
                     dt/muz/dx*(et(2:ib,jb,2:ke)-et(1:ie,jb,2:ke))+dt/muz/dtheta/(r(1:ie,jb,2:ke)+dx/2).*(er(1:ie,2,2:ke)-er(1:ie,jb,2:ke));
               
%***********************************************************************
%% restroe the imformation at the picked point
%***********************************************************************
gVer1(n)=ez(10,9,10);

end

ff=abs(gVer1);
%plot(ff)
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-6-10 15:39:35  | 显示全部楼层
%COXIAL
clear
  
%***********************************************************************
%     Fundamental constants
%***********************************************************************
  
cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
epsz=1.0/(cc*cc*muz);       %permittivity of free space
ii=sqrt(-1)

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=34;      %number of grid cells in x-direction  r=17mm  rout=34mm
je=20;
ke=15;       %number of grid cells in z-direction  l=7.5mm
  
ib=ie+1;
jb=je+1;
kb=ke+1;

dtheta=2*pi/je;
dx=0.0005;          %space increment of cubic lattice 0.5mm
dt=dx/(2.0*cc);     %time step
nmax=2000;          %total number of time steps
cb=dt/epsz;
db=dt/muz/dx;

er=zeros(ie,jb,kb);
et=zeros(ib,jb,kb);
ez=zeros(ib,jb,ke);
hr=zeros(ib,jb,ke);
ht=zeros(ie,jb,ke);
hz=zeros(ie,jb,kb);
r=zeros(ib,jb,kb);
ez0=zeros(ib,jb);
Nf=1000
fmax=100e+9
df=fmax/Nf

f0=30e+9


ndecay=1/(2*fmax*dt);
n0=1/(2*f0*dt);
t=3.123/dx/ie/2;
%compute the r(i)
for i=1:ib
    r(i,:,:)=dx*ie+dx*(i-1);
end
    for i=1:ib
        ez0(i,:)=besselj(0,t*(ie+i-1)*dx)*bessely(0,t*dx*2*ie)-bessely(0,t*(ie+i-1)*dx)*besselj(0,t*2*dx*ie);
    end
   
    %***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************
  
for n=1:nmax
  
%***********************************************************************
%     Update electric fields
%***********************************************************************
er(1:ie,2:jb,2:ke)=er(1:ie,2:jb,2:ke)+...
                    dt/dtheta/epsz/(r(1:ie,2:jb,2:ke)+dx/2).*(hz(1:ie,2:jb,2:ke)-hz(1:ie,1:je,2:ke))-...
                    dt/epsz/dx*(ht(1:ie,2:jb,2:ke)-ht(1:ie,2:jb,1:ke-1));

ez(2:ie,2:jb,1:ke)=ez(2:ie,2:jb,1:ke)+...
                   dt/epsz/dx*(ht(2:ie,2:jb,1:ke)-ht(1:ie-1,2:jb,1:ke))+dt/epsz/2*(ht(2:ie,2:jb,1:ke)+ht(1:ie-1,2:jb,1:ke))./r(2:ie,2:jb,1:ke)-...
                   dt/dtheta/epsz*(hr(2:ie,2:jb,1:ke)-hr(2:ie,1:je,1:ke))./r(2:ie,2:jb,1:ke);

et(2:ie,1:jb,2:ke)=et(2:ie,1:jb,2:ke)+dt/epsz/dx*(hr(2:ie,1:jb,2:ke)-hr(2:ie,1:jb,1:ke-1)-hz(2:ie,1:jb,2:ke)+hz(1:ie-1,1:jb,2:ke));
                  

   %add excite**************************depend upon the mode which compute.
         ez(5,5,3)=exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(5,5,3);
        %ez(:,:,3)=ez0(:,:)*exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(:,:,3);
   %************* periodic  boundary condition   for theta direction *****************************

er(1:ie,1,2:ke)=er(1:ie,1,2:ke)+...
                    dt/dtheta/epsz/(r(1:ie,1,2:ke)+dx/2).*(hz(1:ie,1,2:ke)-hz(1:ie,je,2:ke))-...
                   dt/epsz/dx*(ht(1:ie,1,2:ke)-ht(1:ie,1,1:ke-1));

ez(2:ie,1,1:ke)=ez(2:ie,1,1:ke)+...
                   dt/epsz/dx*(ht(2:ie,1,1:ke)-ht(1:ie-1,1,1:ke))+dt/epsz/2*(ht(2:ie,1,1:ke)+ht(1:ie-1,1,1:ke))./r(2:ie,1,1:ke)-...
                   dt/dtheta/epsz*(hr(2:ie,1,1:ke)-hr(2:ie,je,1:ke))./r(2:ie,1,1:ke);


%***********************************************************************
%     Update magnetic fields
%***********************************************************************
hr(2:ie,1:je,1:ke)=hr(2:ie,1:je,1:ke)-...
                    dt/muz/dtheta*(ez(2:ie,2:jb,1:ke)-ez(2:ie,1:je,1:ke))./r(2:ie,1:je,1:ke)+dt/muz/dx*(et(2:ie,1:je,2:kb)-et(2:ie,1:je,1:ke));
  
ht(1:ie,1:jb,1:ke)=ht(1:ie,1:jb,1:ke)+...
                   dt/muz/dx*(er(1:ie,1:jb,1:ke)-er(1:ie,1:jb,2:kb)+ez(2:ib,1:jb,1:ke)-ez(1:ie,1:jb,1:ke));

hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)-dt/muz/2/(r(1:ie,1:je,2:ke)+dx/2).*(et(1:ie,1:je,2:ke)+et(2:ib,1:je,2:ke))-...
                     dt/muz/dx*(et(2:ib,1:je,2:ke)-et(1:ie,1:je,2:ke))+dt/muz/dtheta/(r(1:ie,1:je,2:ke)+dx/2).*(er(1:ie,2:jb,2:ke)-er(1:ie,1:je,2:ke));  
                     
%************* periodic  boundary condition   for theta direction *****************************
hr(2:ie,jb,1:ke)=hr(2:ie,jb,1:ke)-...
                    dt/muz/dtheta*(ez(2:ie,2,1:ke)-ez(2:ie,jb,1:ke))./r(2:ie,jb,1:ke)+dt/muz/dx*(et(2:ie,jb,2:kb)-et(2:ie,jb,1:ke));
               
hz(1:ie,jb,2:ke)=hz(1:ie,jb,2:ke)-dt/muz/2/(r(1:ie,jb,2:ke)+dx/2).*(et(1:ie,jb,2:ke)+et(2:ib,jb,2:ke))-...
                     dt/muz/dx*(et(2:ib,jb,2:ke)-et(1:ie,jb,2:ke))+dt/muz/dtheta/(r(1:ie,jb,2:ke)+dx/2).*(er(1:ie,2,2:ke)-er(1:ie,jb,2:ke));
               
%***********************************************************************
%% restroe the imformation at the picked point
%***********************************************************************
gVer1(n)=ez(10,9,10);

end

ff=abs(gVer1);
%plot(ff)
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-6-10 15:40:36  | 显示全部楼层
我把楼主的程序贴出来，大家需要下载的可以到一楼下载
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-6-10 17:38:25  | 显示全部楼层
辛苦了，不错
给一部分懒人提供方便，^_^
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-6-10 17:42:13  | 显示全部楼层

回复 沙发 的帖子

如果把数组改成动态数组最好。
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

skygardon 该用户已被删除
发表于 2008-6-10 22:57:22  | 显示全部楼层
提示: 作者被禁止或删除 内容自动屏蔽
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-6-11 08:12:08  | 显示全部楼层

回复 16楼 的帖子

编程也没什么，多磨几次，编程就提高了
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-11-24 13:28:41  | 显示全部楼层
希望能如楼上所说多练就可以掌握门道。
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-11-24 15:12:46  | 显示全部楼层
：23de         我要编程 呵呵 支持
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

发表于 2008-11-30 20:13:42  | 显示全部楼层
上面的代码作为参考很不错，但不能完全用他的代码作为自己的研究，这样更累
发现有不少学生原封不动地抄袭上面的代码作为作业提交，这样对自己的学习很不负责任，当以后真正做研究要用这部分的知识时，会后悔当初没有真正来学
以己之微·网博天下：博览微网之术·创造成功之路！
回复 支持 反对

举报

下一页 »
1234下一页
返回列表 发新帖

发表回复
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 注册

本版积分规则

返回列表 客服中心 搜索
关于我们
关于我们
关注我们
联系我们
帮助中心
资讯中心
企业生态
社区论坛
服务支持
资源下载
售后服务
推广服务
关注我们
官方微博
官方空间
官方微信
诚信网站 可信网站 Powered by Discuz! X3.5 © 2008-2024 头雁微网 版权所有 蜀ICP备2023016130号-1 客服QQ: 68870544  技术支持: 克米设计
快速回复 返回顶部 返回列表