1D FDTD plane wave propagation simulation

I'm currently trying to simulate a simple case of wave propagation in free space before adding in more complexities, and already I'm stumped. Most examples I see online set the Courant stability condition to 1, i.e. the wave travels 1 spatial step in 1 time step. However, I set mine to 2 with the intention of increasing it further should I require more resolution.

Problem: My source is a unity amplitude sine wave (in exponential form), and when I plot it the amplitude is around 2 with increasing ringing with increasing distance.

I'm using these slides as my reference (https://empossible.net/wp-content/uploads/2020/01/Lecture-Formulation-of-1D-FDTD.pdf). The relevant slides are on p.12 and 20.

c0=299792458; %speed of light

u0 = 1.256637e-6;

e0 = 8.85418782e-12;

n1 = 1.; er1=n1^2; %refractive index and relative permittivity

fmax = 1e6; %frequency

lambda_min = c0/(fmax*n1); %minimum wavelength

z_step = lambda_min/(100); %step size in space

t_step = z_step/(2*c0); %Courant stability condition

t_total = t_step*4000; %total runtime

T = t_total/t_step; %total number of timesteps

Z = 400; %total number of spatial steps

Zarray = linspace(1,Z,Z);

Hx = zeros(Z,1); Ey = zeros(Z,1);

C_Ey = ones(Z,1).*c0*t_step/z_step/er1;

C_Hx = ones(Z,1).*c0*t_step/z_step;

H1=0;H2=0;E1=0;E2=0;

for t_count = 1:T

disp(t_count)

H2=H1; H1=Hx(1);

for z_count = 1:Z-1

Hx(z_count) = Hx(z_count) + C_Hx(z_count)*(Ey(z_count+1)-Ey(z_count));

end

Hx(Z) = Hx(Z) + C_Hx(Z)*(E1-Ey(Z));

E2=E1; E1=Ey(Z);

Ey(1) = Ey(1) + C_Ey(1)*(Hx(1)-H1);

for z_count = 2:Z

Ey(z_count) = Ey(z_count) + C_Ey(z_count)*(Hx(z_count)-Hx(z_count-1));

end

Ey(50) = Ey(50)+exp(1i*2*pi*fmax*t_count*t_step); %soft source

plot(Zarray,real(Ey))

drawnow

end

Answer 1

[Romain]

Hi Jerry,

At line 33, change:

Ey(50) = Ey(50)+exp(1i*2*pi*fmax*t_count*t_step);             %soft source

To:

Ey(1) = Ey(1)+exp(1i*2*pi*fmax*t_count*t_step);             %soft source

And you get a much smoother wave:

Additionally, a small error is present in the definition of lambda_min: replace n2 with n1 or define n2 before.

If you sorely need the wave to generate at abscissa 50, do not perform the above modification but change the for range at lines 22 and 29:

 for z_count = 50:Z-1
        Hx(z_count) = Hx(z_count) + C_Hx(z_count)*(Ey(z_count+1)-Ey(z_count));
    end

for z_count = 51:Z
        Ey(z_count) = Ey(z_count) + C_Ey(z_count)*(Hx(z_count)-Hx(z_count-1));
    end

And you will get:

[Jerry Yeung]

Thanks for the reply. Unfortunately this does not solve my issue. In your figures the wave has a peak magnitude of >20, whereas my source has a magnitude of unity. This is due to reflection from the left boundary, but you move the source anywhere within the grid and the rightward propagating wave still has a magnitude of 2. Furthermore, I require 1) my source to not be next to the boundaries and 2) simulation of the back-reflected (leftward-travelling) field later down the line. I can't simply not simulate anything before abscissa 50.

[Romain]

Hi Jerry,

Is there a reason why you generate your sine wave in exponential form ?

How do you deal with boundaries Z=50 and Z=400 ?

Hx(Z) = Hx(Z) + C_Hx(Z)*(E1-Ey(Z));

But you set a each iteration

E1 = Ey(Z)

So Hx(Z) is never modified. Same goes for Ey(50), where you set (I deliberately changed indexes from 1 to 50):

H1 = Hx(50)

And then

Ey(50) = Ey(50) + C_Ey(50)*(Hx(50)-H1)

Looking at the last slide of the pdf you mentioned:

Should not Ey and Hx be equal to 0 outside of the boundaries ? So E1=0 and H1=0 at any time ?

If you set it to 0, I get the reflecting behavior you asked for