-
Notifications
You must be signed in to change notification settings - Fork 1
/
meep_pcsel_latest.py
274 lines (225 loc) · 10.8 KB
/
meep_pcsel_latest.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
'''meep code for simulating PCSEL, generated by gpt.
EXTREMELY SENSITIVE: DO NOT SHARE OR DISTRIBUTE.
AUTHORS: Renjie Li, Sixuan Mao, NOEL, CUHKSZ, June 2023. '''
#June 25th 2023 update: this code contains fake/madeup functions generated by GPT (outrageous /:). Need to fix those manually.
import meep as mp
import numpy as np
import matplotlib.pyplot as plt
# PCSEL structure parameterspythoi
a=1
resolution = 35 # pixels per micron
wavelength = 0.98 # operating wavelength (in microns)
fcen=1/wavelength #central frequency
df = 0.2*fcen #The frequency width of the source
pml_thickness = 0.5 # thickness of perfectly matched layer (PML)
hole_radius = 0.1 # radius of the air holes (in microns)
lattice_spacing = 0.3 # spacing between lattice points (in microns)
active_layer_thickness = 0.5 # thickness of the active layer (in microns)
pc_thickness = active_layer_thickness
n_sub_thickness = 0.2 # thickness of n-substrate layer (in microns)
n_clad_thickness = 0.3 # thickness of n-cladding layer (in microns)
p_clad_thickness = 0.3 # thickness of p-cladding layer (in microns)
# Define the materials
n_inp = 3.4 # refractive index of InP
n_air = 1.0 # refractive index of air
n_substrate = 3.0 # refractive index of substrate
n_cladding = 3.2 # refractive index of cladding (n-type)
p_cladding = 3.2 # refractive index of cladding (p-type)
# Define the PCSEL structure
geometry = [
mp.Block(
size=mp.Vector3(mp.inf, mp.inf, n_sub_thickness),
center=mp.Vector3(0, 0, 0),
material=mp.Medium(index=n_substrate)
),
mp.Block(
size=mp.Vector3(mp.inf, mp.inf, n_clad_thickness),
center=mp.Vector3(0, 0, -n_clad_thickness / 2 - n_sub_thickness / 2),
material=mp.Medium(index=n_cladding)
),
mp.Block(
size=mp.Vector3(mp.inf, mp.inf, active_layer_thickness),
center=mp.Vector3(0, 0, -n_clad_thickness - active_layer_thickness / 2 - n_sub_thickness / 2),
material=mp.Medium(index=n_inp)
),
mp.Block(
size=mp.Vector3(mp.inf, mp.inf, pc_thickness),
center=mp.Vector3(0, 0,
-n_clad_thickness - active_layer_thickness - pc_thickness / 2 - n_sub_thickness / 2),
material=mp.Medium(index=n_inp)
),
mp.Block(
size=mp.Vector3(mp.inf, mp.inf, p_clad_thickness),
center=mp.Vector3(0, 0,
-n_clad_thickness - active_layer_thickness - pc_thickness - p_clad_thickness / 2 - n_sub_thickness / 2),
material=mp.Medium(index=p_cladding)
),
]
# Add photonic crystal layer with air holes
n_side=50
num_holes_x = 50
num_holes_y = 50
hole_spacing_x = lattice_spacing
hole_spacing_y = lattice_spacing
pc_x = (num_holes_x+0) * hole_spacing_x
pc_y = (num_holes_y+0) * hole_spacing_y
pc_z=pml_thickness*2 + pc_thickness + n_sub_thickness + n_clad_thickness + active_layer_thickness + p_clad_thickness # size in z direction
pc_x_center = 0
pc_y_center = 0
for i in range(num_holes_x):
for j in range(num_holes_y):
hole_center = mp.Vector3(
pc_x_center - pc_x / 2 + hole_spacing_x * i + hole_spacing_x / 2,
pc_y_center - pc_y / 2 + hole_spacing_y * j + hole_spacing_y / 2,
-n_clad_thickness - active_layer_thickness - active_layer_thickness/2 - n_sub_thickness/2
)
geometry.append(
mp.Cylinder(radius=hole_radius, height=pc_thickness, center=hole_center, material=mp.air)
)
# Set up the simulation parameters
cell_size = mp.Vector3(pc_x+2*pml_thickness, pc_y+2*pml_thickness, pc_z)
boundary_layers = [mp.PML(pml_thickness)] # boundary conditions
#k_points = mp.interpolate(10, [mp.Vector3(), mp.Vector3(0.5, 0, 0)]) # k-points for Bloch boundary conditions
k_points = mp.Vector3(0,0,0)
# Define the source
source_position = mp.Vector3(0, 0)
source = mp.Source(mp.GaussianSource(wavelength, fwidth=0.2 * wavelength), component=mp.Hz, center=source_position)
src_cmpt = mp.Hz
if src_cmpt == mp.Hz:
symmetries = [mp.Mirror(mp.X,phase=-1),
mp.Mirror(mp.Y,phase=-1)]
elif src_cmpt == mp.Ez:
symmetries = [mp.Mirror(mp.X,phase=+1),
mp.Mirror(mp.Y,phase=+1)]
else:
symmetries = []
# Set up the simulation
sim = mp.Simulation(
resolution=resolution,
geometry=geometry,
cell_size=cell_size,
boundary_layers=boundary_layers,
k_point=k_points,
symmetries=symmetries,
sources=[source]
)
#at the very top n_sub layer
resonance_top = mp.FluxRegion(center=mp.Vector3(0.45*n_side*lattice_spacing,0.45*n_side*lattice_spacing,n_sub_thickness/2),
size=mp.Vector3(0.8*n_side*lattice_spacing,0.8*n_side*lattice_spacing),direction=mp.Z)
#resonance_side = mp.FluxRegion(center=mp.Vector3(0.5*n_side*lattice_spacing,0,0),size=mp.Vector3(0,n_side*lattice_spacing,pc_z-2*pml_thickness),direction=mp.X)
pt = mp.Vector3(z=-n_sub_thickness/2-n_clad_thickness-active_layer_thickness-pc_thickness/2) #at the PC layer
#mode_number = 1 # Mode number to monitor
#mode_monitor = sim.add_mode_monitor(1/wavelength, 0.2/wavelength, 41, mp.ModeRegion(volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(1, 1, 0))))
#dft = sim.add_dft_fields([mp.Ez], 1/wavelength, 0.2/wavelength, 41 , center=mp.Vector3(), size=mp.Vector3(1, 1, 1))
# Run the simulation
#sim.run(until_after_sources=200)
# add DFT monitor
resonance_z = sim.add_flux(fcen, df*2, 41, resonance_top)
#a box surrounding the top n_sub layer
nearfield_box = sim.add_near2far(fcen, df, 41,
mp.Near2FarRegion(center=mp.Vector3(z=0.5*n_sub_thickness/2),
size=mp.Vector3(pc_x,pc_y,0)),
mp.Near2FarRegion(center=mp.Vector3(z=0,x=-pc_x/2),
size=mp.Vector3(0,pc_y,0.5*n_sub_thickness),weight=-1),
mp.Near2FarRegion(center=mp.Vector3(z=0,x=pc_x/2),
size=mp.Vector3(0,pc_y,0.5*n_sub_thickness)),
mp.Near2FarRegion(center=mp.Vector3(z=0,y=-pc_y/2),
size=mp.Vector3(pc_x,0,0.5*n_sub_thickness),weight=-1),
mp.Near2FarRegion(center=mp.Vector3(z=0,y=pc_y/2),
size=mp.Vector3(pc_x,0,0.5*n_sub_thickness))
)
# Run thei simulation
harminv = mp.Harminv(mp.Hz, pt, fcen, 2*df, mxbands = 8)
sim.run(harminv, until_after_sources=200/fcen)
sim.run(mp.in_volume(mp.Volume(center=mp.Vector3(x=pc_x/8,y=pc_y/8,z=-n_sub_thickness/2-n_clad_thickness-active_layer_thickness/2),size=mp.Vector3(0.7*pc_x,0.7*pc_y))),
until_after_sources=1/fcen)
#compute mode metrics: Q, V ... using Harmonic inversion method (Harminv)
first_mode = harminv.modes[0] #the 1st order mode
Q = first_mode.Q #extract the Q factor of the 1st mode
freq = first_mode.freq #extract the freq of the 1st mode
Qlist = [m.Q for m in harminv.modes] #store all Q factors
flist = [m.freq for m in harminv.modes] #store all freq
vol = mp.Volume(center=mp.Vector3(x=pc_x/8,y=pc_y/8,z=-n_sub_thickness/2-n_clad_thickness-active_layer_thickness),size=mp.Vector3(0.7*pc_x,0.7*pc_y,active_layer_thickness+pc_thickness))
V = sim.modal_volume_in_box(box=vol) #calculate the modal volume
# transmittance in z direction
z_flux = mp.get_fluxes(resonance_z)
# x_flux = mp.get_fluxes(resonance_x)
flux_freqs = mp.get_flux_freqs(resonance_z)
wl = []
Trans = []
for i in range(len(flux_freqs)):
wl = np.append(wl, 1 / flux_freqs[i])
Trans = np.append(Trans, z_flux[i])
plt.plot(wl,Trans,label='z')
#plt.plot(wl,Ts_x,label='x')
plt.xlabel("wavelength ( 1/4 m)")
plt.ylabel('Transmittance (a.u.)')
plt.savefig('./Transmittance.png')
#plt.show()
# near field Poynting vector in z direction
(x, y, z, w) = sim.get_array_metadata(dft_cell=resonance_z)
Pz = []
i = 0
for _ in resonance_z.freq:
(Ex, Ey, Hx, Hy) = [sim.get_dft_array(resonance_z, c, i) for c in [mp.Ex, mp.Ey, mp.Hx, mp.Hy]]
flux_density = np.real(np.conj(Ex) * Hy - np.conj(Ey) * Hx) # array
flux = np.sum(w * flux_density) # scalar
Pz.append(flux)
i += 1
# Extract electric field data (Ey component)
ey_data = sim.get_array(center=mp.Vector3(), size=cell_size, component=mp.Ey)
ey_data = ey_data[0, :, :].squeeze()
# Define the grid coordinates
x = np.linspace(-pc_x / 2, pc_x / 2, ey_data.shape[0])
y = np.linspace(-pc_y / 2, pc_y / 2, ey_data.shape[1])
# Plot the electric field profile (Ey component)
plt.figure(figsize=(8, 6))
plt.imshow(ey_data.transpose(), interpolation='spline36', cmap='RdBu', alpha=0.9, extent=[x.min(), x.max(), y.min(), y.max()])
plt.xlabel('x (microns)')
plt.ylabel('y (microns)')
plt.title('Electric Field Profile (|Ey|^2)')
plt.colorbar(label='Electric Field Intensity')
plt.axis('off')
#plt.show()
plt.savefig('./Electrical_field.png')
# Print the computed values
print("Q factor:", Q)
#print("FWHM:", FWHM) #ToDo
print("Resonance wavelength:", freq)
print("Emitted power:", max(Pz))
print("Mode volume:", V)
#compute the far field
# half side length of far-field square box
r = 1e3 # 1mm
# resolution of far fields (points/μm)
res_ff = 50/(2*r) # about 200 points per 2r, total 200*200 points
far_field = sim.get_farfields(nearfield_box,res_ff,center=mp.Vector3(z=r),size=mp.Vector3(pc_x+2*r,pc_y+2*r))
# Define a large sphere (radius=1000) to compute the far-field
farfield_box = mp.Volume(center=mp.Vector3(0, 0, 0), size=mp.Vector3(1000, 1000, 0))
# Run the simulation for some time to let it reach steady state
sim.run(until=200)
# Use the add_flux function to compute the E field on the large sphere
E_far = sim.add_flux(1/wvl_min, 0, 1, farfield_box)
# Run the simulation for some more time to let it reach steady state
sim.run(until=200)
# Get the far field data
far_field_data = sim.get_farfield(E_far, farfield_box)
# Calculate the angles and the radiation pattern
angles = np.linspace(0, 2*np.pi, len(far_field_data))
radiation_pattern = np.abs(far_field_data)
# Find the max radiation direction
max_index = np.argmax(radiation_pattern)
max_angle = angles[max_index]
# Find the FWHM in the far field to get the divergence
half_max = max(radiation_pattern) / 2.
indices = np.where(radiation_pattern > half_max)[0]
fwhm_angle = angles[indices[-1]] - angles[indices[0]]
print("Divergence Angle: ", fwhm_angle)
# Plot the far-field pattern
# plt.figure()
# plt.plot(np.degrees(far_field[:, 0]), far_field[:, 1])
# plt.xlabel("Angle (degrees)")
# plt.ylabel("Power (arb. units)")
# plt.title("Far-Field Pattern")
# #plt.show()
# plt.savefig('./Far_field.png')