1 D FDTD for a plane harmonic wave entering a glass sheet with 2.25 relative permittivity

A plane harmonic wave at a wavelength of λ = 650nm (red light) travels in air along a unique direction. It then impinges on a glass sheet with a relative permittivity of 2.25 in a thickness of d at a normal angle. The medium behind the glass sheet is again air. Using FDTD, show the 1D E/M field distributions of the plane harmonic wave when d = λ, λ/2, λ/4 inside a computation window that includes the entire air-glass-air system. The starting and ending point of the computation window can be chosen at places, e.g., 5λ away from the glass sheet, i.e., the computation window size can be chosen as 5λ (air) + d (glass) + 5λ (air)

I have the code below but I don't understand how to introduce the glass thing in the code. I will be very thankful if someone could help me out.

%Clearing variables in memory and Matlab command screen
clear all;
clc;

% Grid Dimension in x (xdim) direction
xdim=200;

%Total no of time steps
time_tot=350;

%Position of the source (center of the domain)
xsource=3;

%Courant stability factor
S=1;

% Parameters of free space (permittivity and permeability and speed of
% light) are all not 1 and are given real values
epsilon0=(1/(36*pi))*1e-9;
mu0=4*pi*1e-7;
c=3e+8;

% Spatial grid step length (spatial grid step = 1 micron)
delta=1e-6;
% Temporal gris step obtained using Courant condition
deltat=S*delta/c;

% Initialization of field vectors
Ez=zeros(1,xdim);
Hy=zeros(1,xdim);

% Initialization of permittivity and permeability vectors
epsilon=epsilon0*ones(1,xdim);
mu=mu0*ones(1,xdim);

% Initializing electric and magnetic field vectors
sigma=4e-4*ones(1,xdim);
sigma_star=4e-4*ones(1,xdim);

% The user can give a frequency of his choice for sinusoidal (if sine=1 above) waves in Hz 
frequency=1.5e+13;

%Multiplication factor vectors for H vector update to avoid being calculated many times 
%in the time update loop so as to increase computation speed
A=((mu-0.5*deltat*sigma_star)./(mu+0.5*deltat*sigma_star)); 
B=(deltat/delta)./(mu+0.5*deltat*sigma_star);

%Multiplication factor vectors for E vector update to avoid being calculated many times 
%in the time update loop so as to increase computation speed                          
C=((epsilon-0.5*deltat*sigma)./(epsilon+0.5*deltat*sigma)); 
D=(deltat/delta)./(epsilon+0.5*deltat*sigma);                     


% Update loop begins
for n=1:1:time_tot
    
    % Setting time dependent boundaries to update only relevant parts of the 
    % vector where the wave has reached to avoid unnecessary updates.
    if n<xsource-2
        n1=xsource-n-1;
    else
        n1=1;
    end
    if n<xdim-1-xsource
        n2=xsource+n;
    else
        n2=xdim-1;
    end
    
    %Vector Update instead of for loop for Hy field incorporating magnetic
    %conductivity
    Hy(n1:n2)=A(n1:n2).*Hy(n1:n2)+B(n1:n2).*(Ez(n1+1:n2+1)-Ez(n1:n2));
    
    %Vector Update instead of for loop for Ez field incorporating magnetic
    %conductivity
    Ez(n1+1:n2+1)=C(n1+1:n2+1).*Ez(n1+1:n2+1)+D(n1+1:n2+1).*(Hy(n1+1:n2+1)-Hy(n1:n2));
    
    
    tstart=1;
    N_lambda=c/(frequency*delta);
    Ez(xsource)=sin(((2*pi*(c/(delta*N_lambda))*(n-tstart)*deltat)));
    
    
    
    % Mur's Absorbing boundary condition
    if n>1
        Ez(xdim)=dum2+((S-1)/(S+1))*(Ez(xdim-1)-Ez(xdim));
        Ez(2)=dum1+((S-1)/(S+1))*(Ez(2)-Ez(1));
        Ez(1)=Ez(2);
    end
    dum1=Ez(3);
    dum2=Ez(xdim-1);
    
    
    plot((1:1:xdim)*delta,Ez,'color','k');
    titlestring=['\fontsize{20}Plot of Ez vs x for 1D FDTD Mur Absorbing Boundary Condition at time = ',num2str(round(n*deltat/10e-15)),' fs'];
    title(titlestring,'color','k');
    xlabel('x in m','FontSize',20);
    ylabel('Ez in V/m','FontSize',20);
    set(gca,'FontSize',20);
    axis([0 xdim*delta -3 3]);
    getframe;
end

Answer 1

[Abhimenyu]

Hi Mariam,

From the information shared, I could infer that you are trying to introduce the glass sheet into your FDTD simulation. To achieve this, the permittivity ("epsilon") of the region where the glass sheet is located needs to be modified.

The relative permittivity ("ε_r") of the glass sheet is given as 2.25, and the absolute permittivity (ε) of the glass sheet can be calculated using the formula "ε = ε_r * ε_0", where "ε_0" is the permittivity of free space. For lambda= 650nm, the glass sheet thicknesses as "d = λ, λ/2, and λ/4" are to be used for different simulations. The computation window includes the entire air-glass-air system, with a size chosen as "5λ (air) + d (glass) + 5λ (air)".

Let's assume the center of the glass sheet is at the center of the computational domain for simplicity. The starting and ending indices of the glass sheet in the grid will be first calculated and, then the "epsilon" array is to be modified to reflect the presence of the glass sheet.

Please follow the below-mentioned MATLAB code to include the glass sheet for a thickness of "d = λ" as an example:

% The code before this remains unchanged 

% Initialization of permittivity and permeability vectors
epsilon=epsilon0*ones(1,xdim);
mu=mu0*ones(1,xdim);

% Glass sheet properties
lambda=650e-9; % Wavelength of the incident light in meters
d_glass=lambda; % Thickness of the glass sheet in meters
epsilon_r_glass=2.25; % Relative permittivity of the glass
start_glass=xdim/2 - (d_glass/(2*delta)); % Starting index of the glass sheet
end_glass=xdim/2 + (d_glass/(2*delta)); % Ending index of the glass sheet
epsilon(start_glass:end_glass)=epsilon_r_glass*epsilon0; % Setting permittivity for the glass region

% Initializing electric and magnetic field vectors
sigma=4e-4*ones(1,xdim);
sigma_star=4e-4*ones(1,xdim);

% The user can give a frequency of his choice for sinusoidal waves in Hz
frequency=1.5e+13;

% Multiplication factor vectors for H and E field updates
A=((mu-0.5*deltat*sigma_star)./(mu+0.5*deltat*sigma_star));
B=(deltat/delta)./(mu+0.5*deltat*sigma_star);
C=((epsilon-0.5*deltat*sigma)./(epsilon+0.5*deltat*sigma));
D=(deltat/delta)./(epsilon+0.5*deltat*sigma);

% Placeholder for previous boundary values for Mur's condition
dum1=0;
dum2=0;

% Update loop begins
for n=1:1:time_tot
    
    % Update loop for fields, source, and boundary conditions here
    % Similar to the given code, but with the modified epsilon
    
    % Plotting and visualization code remains the same
    
end

Please repeat the simulation with "d_glass" set to "λ/2" and "λ/4" to observe the effects of different glass thicknesses on the wave propagation.

I hope this helps to resolve your query!