www.pudn.com > FDTD3D_Main.zip > FDTD3D_Main.m, change:2013-01-18,size:6101b


function [Ex,Ey,Ez]=FDTD3D_Main(handles) 
 
global SimRunStop 
 
ie = get(handles.xslider,'Value');       %number of grid cells in x-direction 
je = get(handles.yslider,'Value');       %number of grid cells in y-direction 
ke = get(handles.zslider,'Value');       %number of grid cells in z-direction 
 
ib=ie+1;      
jb=je+1;    
kb=ke+1; 
 
%*********************************************************************** 
%     All Domains Fields Ini. 
%*********************************************************************** 
 
Ex=zeros(ie,jb,kb); 
Ey=zeros(ib,je,kb); 
Ez=zeros(ib,jb,ke); 
Hx=zeros(ib,je,ke); 
Hy=zeros(ie,jb,ke); 
Hz=zeros(ie,je,kb); 
 
%*********************************************************************** 
%     Fundamental constants 
%*********************************************************************** 
 
param.cc=2.99792458e8;            %speed of light in free space 
 
param.muz=4.0*pi*1.0e-7;          %permeability of free space 
param.epsz = 1.0/(param.cc*param.cc*param.muz); %permittivity of free space 
 
%*********************************************************************** 
%     Grid parameters 
%*********************************************************************** 
 
 
param.is = get(handles.xsource,'Value');       %location of z-directed current source 
param.js = get(handles.ysource,'Value');       %location of z-directed current source 
 
param.kobs = floor(ke/2);      %Surface of observation 
 
param.dx = get(handles.CellSize,'Value');;          %space increment of cubic lattice 
param.dt = param.dx/(2.0*param.cc);    %time step 
 
param.nmax = get(handles.TimeStep,'Value');;          %total number of time steps 
 
%*********************************************************************** 
%     Differentiated Gaussian pulse excitation 
%*********************************************************************** 
 
param.rtau=get(handles.tau,'Value')*100e-12; 
param.tau=param.rtau/param.dt; 
param.ndelay=3*param.tau; 
param.Amp = get(handles.sourceamp,'Value')*10e11; 
param.srcconst=-param.dt*param.Amp; 
 
%*********************************************************************** 
%     Material parameters 
%*********************************************************************** 
 
param.eps= get(handles.epsilon,'Value'); 
param.sig= get(handles.sigma,'Value');        
 
%*********************************************************************** 
%     Updating coefficients 
%*********************************************************************** 
 
param.ca=(1.0-(param.dt*param.sig)/(2.0*param.epsz*param.eps))/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps)); 
param.cb=(param.dt/param.epsz/param.eps/param.dx)/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps)); 
param.da=1.0; 
param.db=param.dt/param.muz/param.dx; 
 
%*********************************************************************** 
%     Calling FDTD Algorithm 
%*********************************************************************** 
 
ex=zeros(ib,jb,kb); 
ey=zeros(ib,jb,kb); 
ez=zeros(ib,jb,kb); 
hx=zeros(ib,jb,kb); 
hy=zeros(ib,jb,kb); 
hz=zeros(ib,jb,kb); 
 
[X,Y,Z] = meshgrid(1:ib,1:jb,1:kb); % Grid coordinates 
 
Psim = zeros(param.nmax,1); 
Panl = zeros(param.nmax,1); 
 
       
if ((p.ip == 1)&(p.jp == 0)) 
    x = ceil(ie/p.PN)     
    param.a = [1:x-1:ie-x]; 
    param.b = [x+1:x-1:ie]; 
    param.c = [1,1]; 
    param.d = [je,je]; 
    m2 = 1; 
    for n=1:1:param.nmax 
        for m1=1:1:p.PN                         
                [ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
                [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
        end                 
        [Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n); 
        field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); 
        Dyn_FFT 
        pause(0.01); 
    end     
elseif ((p.ip == 0)&(p.jp == 1)) 
    y = ceil(je/p.PN); 
    param.c = [1:y-1:je-y]; 
    param.d = [y+1:y-1:je]; 
    param.a = [1,1]; 
    param.b = [ie,ie]; 
    m1 = 1; 
    for n=1:1:param.nmax 
        for m2=1:1:p.PN  
                [ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
                [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
        end 
        [Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n); 
        field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); 
        pause(0.01); 
    end 
elseif ((p.ip == 1)&(p.jp == 1))     
    x = ceil(ie/p.PN);     
    param.a = [1:x-1:ie-x]; 
    param.b = [x+1:x-1:ie]; 
    y = ceil(je/p.PN); 
    param.c = [1:y-1:je-y]; 
    param.d = [y+1:y-1:je]; 
 
    for n=1:1:param.nmax 
        for m2=1:1:p.PN  
            for m1=1:1:p.PN  
                [ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
                [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
            end 
        end 
        [Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n); 
        field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); 
        pause(0.01); 
    end     
else         
    param.a = 1; 
    param.b = ie; 
    param.c = 1; 
    param.d = je; 
    m1 = 1;m2=1; 
    for n=1:1:param.nmax 
                [ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
                [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); 
                SimRunStop = get(handles.Stop_sim,'value');    
                    if SimRunStop == 1 
                        h = warndlg('Simulation Run is Stopped by User !','!! Warning !!'); 
                        waitfor(h); 
                        break; 
                    end                           
                [Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n); 
                if n>=2 
                    Panl(n)=Panl(n) + Panl(n-1); 
                end 
                field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); 
                pause(0.01);                 
    end     
end