iwannachill

import numpy as np
from math import exp
from matplotlib import pyplot as plt




def plot_e_field(data, timestep, epsilon_high, epsilon_low, cb, label):
    plt.plot(data, color='k', linewidth=1)
    plt.ylabel('E$_x$', fontsize=14)
    plt.xticks(np.arange(0, 199, step=20))
    plt.xlim(0, 199)
    plt.yticks(np.arange(-0.5, 1.2, step=0.5))
    plt.ylim(-0.5, 1.2)
    plt.text(70, 0.5, 'T = {}'.format(timestep), horizontalalignment='center')

    cool1 = ((0.5 / cb - 1) - (epsilon_low-1)) / (epsilon_high - epsilon_low)
    cool2 = ((0.5 / cb - 1) - (epsilon_high-1)) / (epsilon_low - epsilon_high)
    cool2 = cool2 * (epsilon_low / epsilon_high)
    plt.plot(cool1, 'r--', linewidth=0.75)
    plt.plot(cool2, 'b--', linewidth=0.75)

    plt.xlabel('{}'.format(label))


ke = 200
ex = np.zeros(ke)
hy = np.zeros(ke)

t0 = 40
spread = 12

boundary_low = [0, 0]
boundary_high = [0, 0]

# create dielectric profile
cb = np.ones(ke)
cb = 0.5 * cb
cb_start = 100
epsilon_high = 3.4 ** 2 # (GaAs)
epsilon_low = 2.9 ** 2  # (Al0.7Ga0.3As)

# DBR defined here
start_dbr = 100
layer_thickness = 5  # want it thicc

for i in range(start_dbr, ke, 2 * layer_thickness):
    cb[i:i + layer_thickness] = 0.5 / epsilon_high
    cb[i + layer_thickness:i + 2 * layer_thickness] = 0.5 / epsilon_low

nsteps = 540

# dictionary to keep track of desired points for plotting
plotting_points = [
    {'num_steps': 100, 'data_to_plot': None, 'label': ''},
    {'num_steps': 120, 'data_to_plot': None, 'label': ''},
    {'num_steps': 150, 'data_to_plot': None, 'label': ''},
    {'num_steps': 180, 'data_to_plot': None, 'label': ''},
    {'num_steps': 200, 'data_to_plot': None, 'label': ''},
    {'num_steps': 220, 'data_to_plot': None, 'label': ''},
    {'num_steps': 250, 'data_to_plot': None, 'label': ''},
    {'num_steps': 280, 'data_to_plot': None, 'label': ''},
    {'num_steps': 300, 'data_to_plot': None, 'label': ''},
    {'num_steps': 320, 'data_to_plot': None, 'label': ''},
    {'num_steps': 350, 'data_to_plot': None, 'label': ''},
    {'num_steps': 380, 'data_to_plot': None, 'label': ''},
    {'num_steps': 400, 'data_to_plot': None, 'label': ''},
    {'num_steps': 420, 'data_to_plot': None, 'label': 'FDTDcells'}
]

# main FDTD Loop
for time_step in range(1, nsteps + 1):

    # calculate Ex field
    for k in range(1, ke):
        ex[k] = ex[k] + cb[k] * (hy[k - 1] - hy[k])

    # put a gaussian pulse at the low end
    pulse = exp(-0.5 * ((t0 - time_step) / spread) ** 2)
    ex[5] = pulse + ex[5]

    ex[0] = boundary_low.pop(0)
    boundary_low.append(ex[1])

    ex[ke - 1] = boundary_high.pop(0)
    boundary_high.append(ex[ke - 2])

    # calculate the Hy field
    for k in range(ke - 1):
        hy[k] = hy[k] + 0.5 * (ex[k] - ex[k + 1])

    # solve data at certain points for later plotting
    for plotting_point in plotting_points:
        if time_step == plotting_point['num_steps']:
            plotting_point['data_to_plot'] = np.copy(ex)

# plot the outputs in multiple figures
plt.rcParams['font.size'] = 12
num_plots = len(plotting_points)
plots_per_figure = 3
num_figures = (num_plots + plots_per_figure - 1) // plots_per_figure  # Ceiling division

for fig_num in range(num_figures):
    fig = plt.figure(figsize=(8, 7))
    start_index = fig_num * plots_per_figure
    end_index = min(start_index + plots_per_figure, num_plots)

    for subplot_num, plotting_point in enumerate(plotting_points[start_index:end_index]):
        ax = fig.add_subplot(plots_per_figure, 1, subplot_num + 1)
        plot_e_field(plotting_point['data_to_plot'],
                     plotting_point['num_steps'], epsilon_high, epsilon_low, cb,
                     plotting_point['label'])

    plt.subplots_adjust(bottom=0.1, hspace=0.45)
    plt.show()