一维等离子体仿真求解反射系数和透射.txt

本程序教我们一维等离子体的仿真求解反射系数及其投射问题

clc
clear all
%***********************************************************************
%     Fundamental constants
%***********************************************************************
cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
eps0=1.0/(cc*cc*muz);       %permittivity of free space
freq0=80e9;                 %frequency of source excitation
freqp=28.7e9;               %frequency of plasma resonance                
lambda=cc/freq0;            %wavelength of source excitation
omega0=2.0*pi*freq0;
wp=2.0*pi*freqp;
vc=2*pi*20e9;
%***********************************************************************
%     Grid parameters
%***********************************************************************
Imax=800;                     %number of grid cells in x-direction(electric)
k=51;                         %number of frequency selected
dx=lambda/40.0;               %space increment of 1-D lattice
dt=dx/cc/2;                   %time step
spread=136;                   %width of derivative Gauss pulse 
N0=68;                        %delay                       

%***********************************************************************
%     Updating coefficients for space region with nonpermeable media
%***********************************************************************
cgb=dt/dx/muz;
cga=dt/dx/eps0;
epsw=1;
gama0=wp.^2./vc*dt-(wp./vc).^2*(1-exp(-vc*dt));
dgama0=-(wp./vc).^2*(1-exp(-vc*dt)).^2;
dgama1=-(wp./vc).^2*exp(-vc*dt)*(1-exp(-vc*dt)).^2;
%***********************************************************************
%     Field arrays
%***********************************************************************
Ex1(1:Imax+1)=0.0;
csy1(1:Imax+1)=0.0;
Hy1(1:Imax)=0.0;
Ex2(1:Imax+1)=0.0;
csy2(1:Imax+1)=0.0;
Hy2(1:Imax)=0.0;
real_pt1(1:k)=0;
imag_pt1(1:k)=0;
amp1(1:k)=0;
phase1(1:k)=0;
real_pt2(1:k)=0;
imag_pt2(1:k)=0;
amp2(1:k)=0;
phase2(1:k)=0;
real_pt3(1:k)=0;
imag_pt3(1:k)=0;
amp3(1:k)=0;
phase3(1:k)=0;
high1_m1=0.0;
high1_m2=0.0;
low1_m1=0.0;
low1_m2=0.0;
high2_m1=0.0;
high2_m2=0.0;
low2_m1=0.0;
low2_m2=0.0;
x=linspace(dx,Imax*dx,Imax);
%***********************************************************************
%     Frequency for the fourier transform
%***********************************************************************
for m=1:k
    freq(m)=0e9+(m-1)*20e9/10;
    omegadt(m)=2*pi*freq(m)*dt;
end
%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************
for n=1:9600
%***********************************************************************
%    Update electric fields
%***********************************************************************  
  pulse1=(n-N0)/spread*exp(-4*pi*(n-N0).^2/spread.^2);
  Ex1(2)=pulse1;
  for i=2:Imax
      Ex1(i)=Ex1(i)+(Hy1(i-1)-Hy1(i))*cga;      
   end 
  %Calculate the fourier transform
    for m=1:k
       real_pt1(m)=real_pt1(m)+cos(omegadt(m)*n)*Ex1(299);
       imag_pt1(m)=imag_pt1(m)-sin(omegadt(m)*n)*Ex1(299);
       amp1(m)=sqrt(real_pt1(m).^2+imag_pt1(m).^2);
       phase1(m)=atan2(imag_pt1(m),real_pt1(m));
   end       
%***********************************************************************
%     Absorbing Boundary
%***********************************************************************
   Ex1(Imax+1)=high1_m2;
   high1_m2=high1_m1;
   high1_m1=Ex1(Imax); 
   Ex1(1)=low1_m2;
   low1_m2=low1_m1;
   low1_m1=Ex1(2);
%***********************************************************************
%    Update Magnetic fields
%***********************************************************************
  for j=1:Imax
       Hy1(j)=Hy1(j)-(Ex1(j+1)-Ex1(j))*cgb;
  end
%end
%***********************************************************************
%   the Second TIME-STEPPING LOOP(plasma has been added up)
%***********************************************************************
%***********************************************************************
%    Update electric fields
%***********************************************************************  
 Ex2(2)=pulse1;
 for i=1:Imax+1
     if i>1&i<300       
         Ex2(i)=Ex2(i)+(Hy2(i-1)-Hy2(i))*cga;      
     elseif i>299&i<501 % plasma occuppied
         Ex2(i)=(epsw+dgama0)./(epsw+gama0)*Ex2(i)+1./(epsw+gama0)*csy2(i)+...
         dt./eps0*(epsw+gama0)*((Hy2(i-1)-Hy2(i))./dx);
         csy2(i)=Ex2(i)*dgama1+exp(-vc*dt)*csy2(i);   
     elseif i>500&i<801
         Ex2(i)=Ex2(i)+(Hy2(i-1)-Hy2(i))*cga;   
     end  
 end
 %Calculate the fourier transform
 for m=1:k
     real_pt2(m)=real_pt2(m)+cos(omegadt(m)*n)*Ex2(501);
     imag_pt2(m)=imag_pt2(m)-sin(omegadt(m)*n)*Ex2(501);
     amp2(m)=sqrt(real_pt2(m).^2+imag_pt2(m).^2);
     phase2(m)=atan2(imag_pt2(m),real_pt2(m))*180/pi; 
 end 
   %***********************************************************************
%     Absorbing Boundary
%***********************************************************************
   Ex2(Imax+1)=high2_m2;
   high2_m2=high2_m1;
   high2_m1=Ex2(Imax); 
   Ex2(1)=low2_m2;
   low2_m2=low2_m1;
   low2_m1=Ex2(2);
%***********************************************************************
%    Update Magnetic fields
%***********************************************************************
 for j=1:Imax
     Hy2(j)=Hy2(j)-(Ex2(j+1)-Ex2(j))*cgb;
 end
 for m=1:k
     real_pt3(m)=real_pt3(m)+cos(omegadt(m)*n)*(Ex2(299)-Ex1(299));
     imag_pt3(m)=imag_pt3(m)-sin(omegadt(m)*n)*(Ex2(299)-Ex1(299));
     amp3(m)=sqrt(real_pt3(m).^2+imag_pt3(m).^2);
     phase3(m)=atan2(imag_pt3(m),real_pt3(m))*180/pi;
  end 
end
a=(amp3./amp1);
a2=(amp2)./amp1;
s11=20*log10(a);
s21=20*log10(a2);
figure(1)
plot(freq,amp1,'-og',freq,amp3,'-vb',freq,s11,'-*r');  
figure(2)
plot(freq,s21,'-ro') 
figure(3)
plot(freq,phase3,'-ro') 
figure(4)
plot(freq,phase2,'-go')