通信基于matlab FDTD法研究移动通信终端电磁辐射对人体的影响【含Matlab源码 761期】.zip

%*********************************************************************** % 3-D FDTD code with PML absorbing boundary conditions %*********************************************************************** function fdtd3D clear clc %*********************************************************************** % 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 etaz=sqrt(muz/epsz); freq=9e8; lambda=cc/freq; omega=2.0*pi*freq; %*********************************************************************** % Grid parameters %*********************************************************************** ie=80; %number of grid cells in x-direction if mod(ie,2)==0 ie=ie+1; end %Force it to be a odd integer! je=80; %number of grid cells in y-direction if mod(je,2)==0 je=je+1; end %Force it to be a odd integer! ke=80; %number of grid cells in z-direction if mod(ke,2)==0 ke=ke+1; end %Force it to be a odd integer! ib=ie+1; jb=je+1; kb=ke+1; is=45; %location of hard source in x axis js=25; %location of hard source in y axis ks=50; %location of hard source in z axis ds=lambda/20; %space increment of square lattice dt=ds/(2.0*cc); %time step nmax=800; iebc=8; %thickness of left and right PML region jebc=8; %thickness of front and back PML region kebc=8; %thickness of bottom and top PML region ibbc=iebc+1; jbbc=jebc+1; kbbc=kebc+1; iefbc=ie+2*iebc; jefbc=je+2*jebc; kefbc=ke+2*kebc; ibfbc=iefbc+1; jbfbc=jefbc+1; kbfbc=kefbc+1; rmax=1.0e-7; orderbc=2; %*********************************************************************** % Material parameters %*********************************************************************** media=2; eps=[1.0 1.0]; sig=[0.0 1.0e+7]; mur=[1.0 1.0]; sim=[0.0 0.0]; ds = ds/sqrt(eps(1)*mur(1)); %*********************************************************************** % Wave excitation %*********************************************************************** source=zeros(1,nmax); for n=1:nmax source(n)=15.5*sin(omega*n*dt); end %*********************************************************************** % Field arrays %*********************************************************************** ex=zeros(ie,jb,kb); %fields in main grid 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); exybcf=zeros(iefbc,jebc,kbfbc); %fields in front PML region exzbcf=zeros(iefbc,jebc,kbfbc); eyzbcf=zeros(ibfbc,jebc,kbfbc); eyxbcf=zeros(ibfbc,jebc,kbfbc); ezxbcf=zeros(ibfbc,jebc,kefbc); ezybcf=zeros(ibfbc,jebc,kefbc); hxybcf=zeros(ibfbc,jebc,kefbc); hxzbcf=zeros(ibfbc,jebc,kefbc); hyzbcf=zeros(iefbc,jebc,kefbc); hyxbcf=zeros(iefbc,jebc,kefbc); hzxbcf=zeros(iefbc,jebc,kbfbc); hzybcf=zeros(iefbc,jebc,kbfbc); exybcb=zeros(iefbc,jebc,kbfbc); %fields in back PML region exzbcb=zeros(iefbc,jebc,kbfbc); eyzbcb=zeros(ibfbc,jebc,kbfbc); eyxbcb=zeros(ibfbc,jebc,kbfbc); ezxbcb=zeros(ibfbc,jebc,kefbc); ezybcb=zeros(ibfbc,jebc,kefbc); hxybcb=zeros(ibfbc,jebc,kefbc); hxzbcb=zeros(ibfbc,jebc,kefbc); hyzbcb=zeros(iefbc,jebc,kefbc); hyxbcb=zeros(iefbc,jebc,kefbc); hzxbcb=zeros(iefbc,jebc,kbfbc); hzybcb=zeros(iefbc,jebc,kbfbc); exybcl=zeros(iebc,jb,kbfbc); %fields in left PML region exzbcl=zeros(iebc,jb,kbfbc); eyzbcl=zeros(iebc,je,kbfbc); eyxbcl=zeros(iebc,je,kbfbc); ezxbcl=zeros(iebc,jb,kefbc); ezybcl=zeros(iebc,jb,kefbc); hxybcl=zeros(iebc,je,kefbc); hxzbcl=zeros(iebc,je,kefbc); hyzbcl=zeros(iebc,jb,kefbc); hyxbcl=zeros(iebc,jb,kefbc); hzxbcl=zeros(iebc,je,kbfbc); hzybcl=zeros(iebc,je,kbfbc); exybcr=zeros(iebc,jb,kbfbc); %fields in right PML region exzbcr=zeros(iebc,jb,kbfbc); eyzbcr=zeros(iebc,je,kbfbc); eyxbcr=zeros(iebc,je,kbfbc); ezxbcr=zeros(iebc,jb,kefbc); ezybcr=zeros(iebc,jb,kefbc); hxybcr=zeros(iebc,je,kefbc); hxzbcr=zeros(iebc,je,kefbc); hyzbcr=zeros(iebc,jb,kefbc); hyxbcr=zeros(iebc,jb,kefbc); hzxbcr=zeros(iebc,je,kbfbc); hzybcr=zeros(iebc,je,kbfbc); exybcd=zeros(ie,jb,kebc); %fields in bottom PML region exzbcd=zeros(ie,jb,kebc); eyzbcd=zeros(ib,je,kebc); eyxbcd=zeros(ib,je,kebc); ezxbcd=zeros(ib,jb,kebc); ezybcd=zeros(ib,jb,kebc); hxybcd=zeros(ib,je,kebc); hxzbcd=zeros(ib,je,kebc); hyzbcd=zeros(ie,jb,kebc); hyxbcd=zeros(ie,jb,kebc); hzybcd=zeros(ie,je,kebc); hzxbcd=zeros(ie,je,kebc); exybct=zeros(ie,jb,kebc); %fields in top PML region exzbct=zeros(ie,jb,kebc); eyzbct=zeros(ib,je,kebc); eyxbct=zeros(ib,je,kebc); ezxbct=zeros(ib,jb,kebc); ezybct=zeros(ib,jb,kebc); hxybct=zeros(ib,je,kebc); hxzbct=zeros(ib,je,kebc); hyzbct=zeros(ie,jb,kebc); hyxbct=zeros(ie,jb,kebc); hzxbct=zeros(ie,je,kebc); hzybct=zeros(ie,je,kebc); %*********************************************************************** % Updating coefficients %*********************************************************************** for i=1:media eaf=dt*sig(i)/(2.0*epsz*eps(i)); ca(i)=(1.0-eaf)/(1.0+eaf); cb(i)=dt/epsz/eps(i)/ds/(1.0+eaf); haf=dt*sim(i)/(2.0*muz*mur(i)); da(i)=(1.0-haf)/(1.0+haf); db(i)=dt/muz/mur(i)/ds/(1.0+haf); end %*********************************************************************** % Geometry specification (main grid) %*********************************************************************** % Initialize entire main grid to free space caex(1:ie,1:jb,1:kb)=ca(1); cbex(1:ie,1:jb,1:kb)=cb(1); caey(1:ib,1:je,1:kb)=ca(1); cbey(1:ib,1:je,1:kb)=cb(1); caez(1:ib,1:jb,1:ke)=ca(1); cbez(1:ib,1:jb,1:ke)=cb(1); dahx(1:ib,1:je,1:ke)=da(1); dbhx(1:ib,1:je,1:ke)=db(1); dahy(1:ie,1:jb,1:ke)=da(1); dbhy(1:ie,1:jb,1:ke)=db(1); dahz(1:ie,1:je,1:kb)=da(1); dbhz(1:ie,1:je,1:kb)=db(1); % % Add cylinder %*********************************************************************** % Fill the PML regions %*********************************************************************** delbc=iebc*ds; sigmam = -(orderbc+1)*log(rmax)/(2*etaz*sqrt(mur(1)/eps(1))*delbc); bcfactor = sigmam/(ds*(orderbc+1)*(delbc^orderbc)); sigmams = sigmam*((muz*mur(1))/(epsz*eps(1))); bcfactors = sigmams/(ds*(orderbc+1)*(delbc^orderbc)); % FRONT region % front face caexybcf(1:iefbc,1,1:kbfbc)=1.0; cbexybcf(1:iefbc,1,1:kbfbc)=0.0; caexzbcf(1:iefbc,1,1:kbfbc)=1.0; cbexzbcf(1:iefbc,1,1:kbfbc)=0.0; caezxbcf(1:ibfbc,1,1:kefbc)=1.0; cbezxbcf(1:ibfbc,1,1:kefbc)=0.0; caezybcf(1:ibfbc,1,1:kefbc)=1.0; cbezybcf(1:ibfbc,1,1:kefbc)=0.0; dahyxbcf(1:iefbc,1,1:kefbc)=1.0; dbhyxbcf(1:iefbc,1,1:kefbc)=0.0; dahyzbcf(1:iefbc,1,1:kefbc)=1.0; dbhyzbcf(1:iefbc,1,1:kefbc)=0.0; % left face caeyxbcf(1,1:jebc,1:kbfbc)=1.0; cbeyxbcf(1,1:jebc,1:kbfbc)=0.0; caeyzbcf(1,1:jebc,1:kbfbc)=1.0; cbeyzbcf(1,1:jebc,1:kbfbc)=0.0; caezxbcf(1,1:jebc,1:kefbc)=1.0; cbezxbcf(1,1:jebc,1:kefbc)=0.0; caezybcf(1,1:jebc,1:kefbc)=1.0; cbezybcf(1,1:jebc,1:kefbc)=0.0; dahxybcf(1,1:jebc,1:kebc)=1.0; dbhxybcf(1,1:jebc,1:kebc)=0.0; dahxzbcf(1,1:jebc,1:kebc)=1.0; dbhxzbcf(1,1:jebc,1:kebc)=0.0; % right face caeyxbcf(ibfbc,1:jebc,1:kbfbc)=1.0; cbeyxbcf(ibfbc,1:jebc,1:kbfbc)=0.0; caeyzbcf(ibfbc,1:jebc,1:kbfbc)=1.0; cbeyzbcf(ibfbc,1:jebc,1:kbfbc)=0.0; caezxbcf(ibfbc,1:jebc,1:kefbc)=1.0; cbezxbcf(ibfbc,1:jebc,1:kefbc)=0.0; caezybcf(ibfbc,1:jebc,1:kefbc)=1.0; cbezybcf(ibfbc,1:jebc,1:kefbc)=0.0; dahxybcf(ibfbc,1:jebc,1:kebc)=1.0; dbhxybcf(ibfbc,1:jebc,1:kebc)=0.0; dahxzbcf(ibfbc,1:jebc,1:kebc)=1.0; dbhxzbcf(ibfbc,1:jebc,1:kebc)=0.0; % bottom face caexybcf(1:iefbc,1:jebc,1)=1.0; cbexybcf(1:iefbc,1:je