一维FDTD色散媒质.docx

% Frequency-dependent media % More on One-dimensional simulation 2_3 % 1D FDTD simulation %KE=200; % space long %ex=zeros(1,KE+1);hy=zeros(1,KE+1);dx=zeros(1,KE+1);ix=zeros(1,KE+1);sx=zeros(1,KE+1);ga=zeros(1,KE+1);gb=zeros(1,KE+1);gc=zeros(1,KE+1);mag=zeros(1,KE+1);arg=zeros(1,3);% Real and imaginary parts of the Fourier Transform %real_pt=zeros(3,KE+1);imag_pt=zeros(3,KE+1);% Amplitude and phase of the Fourier Transform %ampn=zeros(3,KE+1);phasen=zeros(3,KE+1);% Amplitude and phase of input pulse %amp_in=zeros(1,3);phase_in=zeros(1,3);% Fourier Trans. of input pulse %real_in=zeros(1,3);imag_in=zeros(1,3);% Parameters %ex_low_m1=0;ex_low_m2=0;ex_high_m1=0;ex_high_m2=0; T=0; NSTEP=500; % time step %kc=KE/2; % Medium interface %t0=50 ; % Gaussian pulse parameter %spread=10; % Gaussian pulse parameter %ddx=0.01; % set the cell size to 1cm %dt=ddx/(2*3e8); % calculate the time step %freq_in=7*1e8; % signal frequency %epsilon=2; % relative dielectric constant %kstart=100; % position of the dielectric medium %epsz=8.85419e-12; % absolute dielectric constant %sigma=0.01; % conductivity %chil=2;tau=1e-9; % signal frequency %freq=50.e6,200.e6,500.e6; for i=1:3arg(i)=2*pi*freq(i)*dt;end% initialize to free space %for n=1:KE+1 ga(n)=1; gb(n)=0; gc(n)=0;enddel_exp=exp(-dt/tau);for n=kstart:KE+1ga(n)=1/(epsilon+(sigma*dt/epsz)+chil*dt/tau);gb(n)=sigma*dt/epsz;gc(n)=chil*dt/tau;end for n=1:NSTEP T=T+1;% calculate dx field % for k=2:KE+1 dx(k)=dx(k)+0.5*(hy(k-1)-hy(k); end % insert source % pulse=exp(-0.5*(t0-T)2/spread2); % sin(2*pi*freq_in*dt*T); % dx(5)=dx(5)+pulse; % calculate Ex from dx % for k=2:KE+1 ex(k)=ga(k)*(dx(k)-ix(k)-del_exp*sx(k); ix(k)=ix(k)+gb(k)*ex(k); sx(k)=del_exp*sx(k)+gc(k)*ex(k); end% calculate the Fourier transform of Ex % for k=2:KE+1 for i=1:3 real_pt(i,k)=real_pt(i,k)+cos(arg(i)*T)*ex(k); imag_pt(i,k)=imag_pt(i,k)-sin(arg(i)*T)*ex(k); %ampn(i,k)=sqrt( real_pt(i,k)2+imag_pt(i,k)2); %phasen(i,k)=atan2(imag_pt(i,k),real_pt(i,k); endend% Fourier Transform of the input pulse % if T100 for i=1:3 real_in(i)=real_in(i)+cos(arg(i)*T)*ex(10); imag_in(i)=imag_in(i)-sin(arg(i)*T)*ex(10); end end % Calculate the amplitude and phase of each pulse % % Amplitude and phase of the input pulse % for i=1:3 amp_in(i)=sqrt( real_in(i)2+imag_in(i)2); phase_in(i)=atan2(imag_in(i),real_in(i); for k=2:KE+1 ampn(i,k)=(1.0/amp_in(i)*sqrt( real_pt(i,k)2+imag_pt(i,k)2); phasen(i,k)=atan2(imag_pt(i,k),real_pt(i,k)-phase_in(i); endend% absorbing boundary condition % ex(1)=ex_low_m2; ex_low_m2=ex_low_m1; ex_low_m1=ex(2); ex(KE) = ex_high_m2; ex_high_m2 = ex_high_m1; ex_high_m1 = ex(KE-1);% calculate H field % for k=1:KE hy(k)=hy(k)+0.5*(ex(k)-ex(k+1); end% display using matlab % m=moviein(300); i=1:KE+1; subplot(2,2,1); plot(i,ex); axis(0 200 -1.5 1.5) title(电场 t= num2str(n) ; subplot(2,2,2); plot(i,ampn(1,:); axis(0 200 -1.5 1.5) title(ampn Freq.Domain at 100 MHz t= num2str(n) ; subplot(2,2,3); plot(i,ampn(2,:); axis(0 200 -2 2) title(Amplitude Freq.Domain at 200 MHz t= num2