pyMMF vs FDTD simulations (Comsol) for GRIN#
For a parabolic graded index fiber, we compare here the results of pyMMF with the solvers radial, eig, and WKB with FDTD simulations using COMSOL.
[28]:
import matplotlib.pyplot as plt
import numpy as np
import sys
import pyMMF
SMALL_SIZE = 14
MEDIUM_SIZE = 16
BIGGER_SIZE = 18
plt.rcParams.update({'font.size': MEDIUM_SIZE})
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
1. Parameters of the fiber#
Comsol simulation were done using those same parameters
[29]:
NA = 0.2
radius = 25 # in microns
areaSize = 2.4*radius # calculate the field on an area larger than the diameter of the fiber
npoints = 2**7 # resolution of the window
n1 = 1.45
wl = 1.55 # wavelength in microns
curvature = None
k0 = 2.*np.pi/wl
r_max = 3.2*radius
npoints_search = 2**8
dh = 2*radius/npoints_search
2. Compute transverse modes using pyMMF’s axisymmetric solver#
[30]:
profile = pyMMF.IndexProfile(npoints = npoints, areaSize = areaSize)
profile.initParabolicGRIN(n1=n1,a=radius,NA=NA)
[31]:
plt.figure()
plt.imshow(profile.n.reshape([npoints]*2))
plt.title(r'Index porfile')
plt.axis('off')
plt.colorbar()
[31]:
<matplotlib.colorbar.Colorbar at 0x788a975ac230>
[32]:
solver = pyMMF.propagationModeSolver()
solver.setIndexProfile(profile)
solver.setWL(wl)
solver_options = {
'r_max': r_max, # max radius to calculate (and first try for large radial boundary condition)
'dh' : dh, # radial resolution during the computation
'N_beta_coarse' : 1_000, # number of steps of the initial coarse scan
}
modes = solver.solve(
solver='radial',
curvature = None,
options = solver_options
)
2024-09-13 16:11:30,654 - pyMMF.core [DEBUG ] Debug mode ON.
2024-09-13 16:11:30,655 - pyMMF.solv [INFO ] Searching for modes with beta_min=5.821637357564584, beta_max=5.877818513168
2024-09-13 16:11:30,666 - pyMMF.solv [INFO ] Found 5 radial mode(s) for m=0
2024-09-13 16:11:30,666 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,668 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.05066217542815342 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,669 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,669 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,670 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.05066217542815342 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,671 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,671 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,672 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.05066217542815342 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,672 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,672 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,673 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.05066217542815342 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,674 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,675 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.10a
2024-09-13 16:11:30,675 - pyMMF.solv [ERROR ] Field limit 0.012755129839120789 at the founded beta=0.050662175428153364 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,676 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,676 - pyMMF.solv [WARNING] Retrying by changing r_max to 1.89a
2024-09-13 16:11:30,678 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,679 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.03961031292887776 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,680 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,680 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,681 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.03961031292887776 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,681 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,682 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,683 - pyMMF.solv [ERROR ] Field limit 0.056019302460125854 at the founded beta=0.03961031292887755 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,683 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,684 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,686 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 3
2024-09-13 16:11:30,688 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.028541272457537663 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,688 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,688 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,690 - pyMMF.solv [ERROR ] Field limit -0.05253356524115814 at the founded beta=0.028541272457537507 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,691 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,691 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,693 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 4
2024-09-13 16:11:30,696 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 5
2024-09-13 16:11:30,709 - pyMMF.solv [INFO ] Found 5 radial mode(s) for m=1
2024-09-13 16:11:30,710 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,711 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.04513477880645375 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,711 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,712 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,713 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.04513477880645375 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,713 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,714 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,715 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.04513477880645375 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,716 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,716 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,717 - pyMMF.solv [ERROR ] Field limit 0.14218744703877784 at the founded beta=0.045134778806453925 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,717 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,717 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.10a
2024-09-13 16:11:30,718 - pyMMF.solv [ERROR ] Field limit 0.0021314917971200368 at the founded beta=0.045134778806453925 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,719 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,719 - pyMMF.solv [WARNING] Retrying by changing r_max to 1.89a
2024-09-13 16:11:30,722 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,723 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.034070167189183195 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,723 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,724 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,725 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.034070167189183195 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,726 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,727 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,728 - pyMMF.solv [ERROR ] Field limit 0.04819796399196991 at the founded beta=0.03407016718918313 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,729 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,730 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,735 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 3
2024-09-13 16:11:30,736 - pyMMF.solv [ERROR ] Field limit -0.19714650798060052 at the founded beta=0.02298945455284462 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,737 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,738 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,739 - pyMMF.solv [ERROR ] Field limit -0.003288172064295251 at the founded beta=0.02298945455284459 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,739 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,739 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,743 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 4
2024-09-13 16:11:30,745 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 5
2024-09-13 16:11:30,757 - pyMMF.solv [INFO ] Found 4 radial mode(s) for m=2
2024-09-13 16:11:30,758 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,763 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.039606896670381544 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,763 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,764 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,766 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.039606896670381544 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,766 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,767 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,768 - pyMMF.solv [ERROR ] Field limit -0.1675003460169947 at the founded beta=0.039606896670381475 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,768 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,768 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,769 - pyMMF.solv [ERROR ] Field limit -0.0021978546918678364 at the founded beta=0.039606896670381524 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,770 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,771 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.10a
2024-09-13 16:11:30,773 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,774 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.02853435051797808 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,774 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,775 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,776 - pyMMF.solv [ERROR ] Field limit -0.12817439525521637 at the founded beta=0.028534350517978212 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,776 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,776 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,777 - pyMMF.solv [ERROR ] Field limit -0.0020810532272158865 at the founded beta=0.028534350517978292 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,777 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,778 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,781 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 3
2024-09-13 16:11:30,782 - pyMMF.solv [ERROR ] Field limit -0.005271842278740631 at the founded beta=0.017450798239859608 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,782 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,783 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,785 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 4
2024-09-13 16:11:30,796 - pyMMF.solv [INFO ] Found 4 radial mode(s) for m=3
2024-09-13 16:11:30,797 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,798 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.034071439795895594 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,798 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,799 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,799 - pyMMF.solv [ERROR ] Field limit 1.0 at the founded beta=0.034071439795895594 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,800 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,800 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,801 - pyMMF.solv [ERROR ] Field limit -0.06289455472788127 at the founded beta=0.03407143979589552 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,801 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,801 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,803 - pyMMF.solv [ERROR ] Field limit -0.0011262124056772855 at the founded beta=0.034071439795895476 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,803 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,803 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.10a
2024-09-13 16:11:30,807 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,809 - pyMMF.solv [ERROR ] Field limit 0.5055215288595848 at the founded beta=0.02298910620795198 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,809 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,810 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,811 - pyMMF.solv [ERROR ] Field limit 0.00833954404674791 at the founded beta=0.022989106207951988 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,812 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,812 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,815 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 3
2024-09-13 16:11:30,818 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 4
2024-09-13 16:11:30,834 - pyMMF.solv [INFO ] Found 3 radial mode(s) for m=4
2024-09-13 16:11:30,835 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,836 - pyMMF.solv [ERROR ] Field limit -1.0 at the founded beta=0.028530789290099768 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,836 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,836 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,837 - pyMMF.solv [ERROR ] Field limit 0.3411508905024992 at the founded beta=0.028530789290099973 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,837 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,837 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,838 - pyMMF.solv [ERROR ] Field limit 0.005449900686374114 at the founded beta=0.028530789290099997 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,838 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,839 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.33a
2024-09-13 16:11:30,841 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,842 - pyMMF.solv [ERROR ] Field limit 0.003117664564555346 at the founded beta=0.01744070049566171 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,842 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,842 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,845 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 3
2024-09-13 16:11:30,857 - pyMMF.solv [INFO ] Found 3 radial mode(s) for m=5
2024-09-13 16:11:30,857 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,859 - pyMMF.solv [ERROR ] Field limit 0.8563321063878859 at the founded beta=0.022985043200470234 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,859 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,860 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,861 - pyMMF.solv [ERROR ] Field limit 0.013825790229949975 at the founded beta=0.02298504320047025 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,862 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,862 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.59a
2024-09-13 16:11:30,865 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,867 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 3
2024-09-13 16:11:30,890 - pyMMF.solv [INFO ] Found 2 radial mode(s) for m=6
2024-09-13 16:11:30,891 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,894 - pyMMF.solv [ERROR ] Field limit -0.02120303092848433 at the founded beta=0.017434568053403994 is greater than field_limit_tol=0.001
2024-09-13 16:11:30,896 - pyMMF.solv [WARNING] Boundary condition could not be met.
2024-09-13 16:11:30,896 - pyMMF.solv [WARNING] Retrying by changing r_max to 2.88a
2024-09-13 16:11:30,901 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,925 - pyMMF.solv [INFO ] Found 2 radial mode(s) for m=7
2024-09-13 16:11:30,926 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,928 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 2
2024-09-13 16:11:30,939 - pyMMF.solv [INFO ] Found 1 radial mode(s) for m=8
2024-09-13 16:11:30,940 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,951 - pyMMF.solv [INFO ] Found 1 radial mode(s) for m=9
2024-09-13 16:11:30,952 - pyMMF.solv [INFO ] Searching propagation constant for |l| = 1
2024-09-13 16:11:30,967 - pyMMF.solv [INFO ] Found 0 radial mode(s) for m=10
2024-09-13 16:11:30,967 - pyMMF.solv [INFO ] Solver found 55 modes is 0.31 seconds.
2024-09-13 16:11:30,967 - pyMMF.core [DEBUG ] Mode data stored in memory.
Get the mode matrix and rearrange to mode order#
(to fit with comsol data later on)
[33]:
M0_rad = modes.getModeMatrix(npola = 2)
Nmodes = modes.number
new_ind = [i//2 if i%2 == 0 else (i-1)//2+Nmodes for i in range(2*Nmodes)]
M0_rad = M0_rad[:,new_ind]
betas_as = np.sort(np.concatenate([modes.betas]*2))[::-1]
M = modes.m
L = modes.l
3. Compute transverse modes using pyMMF’s 2D finite difference eigenvalue solver#
See wavefrontshapin.net tutorial on solving the Helmholtz discretized equation and on pyMMF implementation.
[34]:
modes = solver.solve(
solver='eig',
options = {'nmodesMax' : Nmodes}
)
2024-09-13 16:11:31,079 - pyMMF.solv [INFO ] Solving the spatial eigenvalue problem for mode finding.
2024-09-13 16:11:31,081 - pyMMF.solv [INFO ] Use close boundary condition.
2024-09-13 16:11:37,034 - pyMMF.solv [INFO ] Solver found 55 modes is 5.95 seconds.
2024-09-13 16:11:37,034 - pyMMF.solv [WARNING] The solver reached the maximum number of modes set.
2024-09-13 16:11:37,035 - pyMMF.solv [WARNING] Some propagating modes may be missing.
2024-09-13 16:11:37,035 - pyMMF.core [DEBUG ] Mode data stored in memory.
[35]:
M0_eig = modes.getModeMatrix(npola = 2)
M0_eig = M0_eig[:,new_ind]
# dupliate the propagation constants for the two polarizations
betas_eig = np.abs(np.sort(np.concatenate([modes.betas]*2))[::-1])
4. WKB approximation#
[36]:
b = radius*n1/NA
f_parabolic = lambda r: np.sqrt(n1**2*(1.-(r/b)**2))
[37]:
r_vec = np.linspace(0, 0.7*areaSize, 150)
real_profile = [profile.radialFunc(r) for r in r_vec]
infinite_GRIN_profile = [f_parabolic(r) for r in r_vec]
plt.figure()
plt.plot(r_vec,real_profile)
plt.plot(r_vec,infinite_GRIN_profile, '--')
[37]:
[<matplotlib.lines.Line2D at 0x788a9759bdd0>]
propagation constants under WKB approximation#
\(\beta_{l,m} = \sqrt{k_o^2 n_1^2-2\alpha \left( |l|+2m+1\right)}\)
\(\alpha = k_o n_1/b\)
and
\(b = \frac{radius \times n_1}{NA}\)
Mode profile under WKB approximation#
$:nbsphinx-math:psi_{l,m}(r, \phi) = A e^{- \frac{\alpha r^2}{2}} (\alpha `r\ :sup:`2){|m|/2} L_l{|m|}(:nbsphinx-math:alpha r2)e^{im:nbsphinx-math:phi} $
\(L_l^{|m|}\) Laguerre polynomial
[38]:
modes = solver.solve(
solver='WKB',
options = {'degenerate_mode' : "sin"}
)
M0_wkb = modes.getModeMatrix(npola = 2)
betas_wkb = np.concatenate([modes.betas]*2)
new_ind = np.argsort(betas_wkb)[::-1]
betas_wkb = betas_wkb[new_ind]
M0_wkb = M0_wkb[:,new_ind]
2024-09-13 16:11:37,186 - pyMMF.solv [INFO ] Finding analytical LP mode profiles associated to the propagation constants.
2024-09-13 16:11:37,186 - pyMMF.solv [INFO ] modes = 55
2024-09-13 16:11:37,225 - pyMMF.solv [INFO ] Found 55 LP mode profiles in 0.0 minutes.
2024-09-13 16:11:37,225 - pyMMF.core [DEBUG ] Mode data stored in memory.
5. Load Comsol data#
We used the same paramters as for pyMMF calculations. Simulation are full vectorial FDTD simulation, so we do have a (small) longitudinal component of the optical field.
[45]:
data_Ex = np.load('data/modes_transverse_comsol_128_Ex.npz')['Ex']
data_Ey = np.load('data/modes_transverse_comsol_128_Ey.npz')['Ey']
data_Ez = np.load('data/modes_transverse_comsol_128_Ez.npz')['Ez']
# matrix of mode profiles
M0_cs = np.concatenate([data_Ex,data_Ey,data_Ez],axis=0)
# propagation constats
betas_cs = np.load('data/modes_transverse_comsol_128_betas.npz')['betas']
Show one mode#
[40]:
ind = 45
Ex = M0_cs[:npoints**2,ind].reshape([npoints]*2)
Ey = M0_cs[npoints**2:2*npoints**2,ind].reshape([npoints]*2)
Ez = M0_cs[2*npoints**2:,ind].reshape([npoints]*2)
max_E = np.max([np.max(np.abs(E)) for E in [Ex,Ey]])
plt.figure(figsize=[12,5])
plt.subplot(131)
plt.imshow(np.abs(Ex), vmax = max_E)
plt.axis('off')
plt.title(r'$|E_x|$')
plt.subplot(132)
plt.imshow(np.abs(Ey), vmax = max_E)
plt.axis('off')
plt.title(r'$|E_y|$')
plt.subplot(133)
plt.imshow(np.abs(Ez), vmax = max_E)
plt.axis('off')
plt.title(r'$|E_z|$')
plt.suptitle(f'Full vectorial FDTD (Comsol)\n Mode {ind}')
plt.figure(figsize=[8,5])
plt.subplot(121)
plt.imshow(np.sqrt(np.abs(Ex)**2+np.abs(Ey)**2), vmax = max_E)
plt.axis('off')
plt.title(r'$|E_t|=\sqrt{|E_x|^2+|E_y|^2}$')
plt.subplot(122)
plt.imshow(np.abs(Ey), vmax = max_E)
plt.axis('off')
plt.title(r'$|E_y|$')
plt.suptitle(f'Full vectorial FDTD (Comsol)\n Mode {ind}')
[40]:
Text(0.5, 0.98, 'Full vectorial FDTD (Comsol)\n Mode 45')
6. Compare propagation constants#
6.1 Whole range#
[16]:
msize = 10
beta_max = k0*n1
beta_min = k0*np.min(profile.n)
plt.figure(figsize = (10,6))
plt.plot(np.abs(betas_cs[:110]),'b+', label = 'Comsol', markersize = msize, markeredgewidth=2)
plt.plot(betas_as,'r*', label = 'pyMMF axisymmetric solver', markersize = msize, markeredgewidth=2)
plt.plot(betas_eig,'gx', label = 'pyMMF 2D eigenvalue solver', markersize = msize, markeredgewidth=2)
plt.plot(betas_wkb,'k.', label = 'WKB', markersize = msize,markeredgewidth=2)
plt.axhline(beta_min, color='g')
plt.axhline(beta_max, color='g')
plt.xlabel('Modes')
plt.ylabel('Propagation constant')
plt.ylim([5.81,5.88])
plt.legend(fontsize=16, loc = 'upper right')
[16]:
<matplotlib.legend.Legend at 0x788a96b835f0>
6.2 Zoom on lower order modes#
[17]:
msize = 14
xlim = [0,10]
ylim = [5.86,5.875]
beta_max = k0*n1
beta_min = k0*np.min(profile.n)
plt.figure(figsize = (10,6))
plt.plot(np.abs(betas_cs[:110]),'b+', label = 'Comsol', markersize = 1.5*msize, markeredgewidth=2)
plt.plot(betas_as,'r*', label = 'pyMMF axisymmetric solver', markersize = msize, markeredgewidth=2)
plt.plot(betas_eig,'gx', label = 'pyMMF 2D eigenvalue solver', markersize = msize, markeredgewidth=2)
plt.plot(betas_wkb,'k.', label = 'WKB', markersize = msize,)
plt.axhline(beta_min, color='g')
plt.axhline(beta_max, color='g')
plt.xlabel('Modes')
plt.ylabel('Propagation constant')
plt.xlim(xlim)
plt.ylim(ylim)
plt.legend(fontsize=16, loc = 'upper right')
plt.title('Low order modes')
[17]:
Text(0.5, 1.0, 'Low order modes')
6.3 Zoom on higher order modes#
[18]:
xlim = [89,110]
ylim = [5.8215,5.824]
beta_max = k0*n1
beta_min = k0*np.min(profile.n)
plt.figure(figsize = (10,6))
plt.plot(np.abs(betas_cs[:110]),'b+', label = 'Comsol', markersize = 1.5*msize, markeredgewidth=2)
plt.plot(betas_as,'r*', label = 'pyMMF axisymmetric solver', markersize = msize, markeredgewidth=2)
plt.plot(betas_eig,'gx', label = 'pyMMF 2D eigenvalue solver', markersize = msize, markeredgewidth=2)
plt.plot(betas_wkb,'k.', label = 'WKB', markersize = msize,)
plt.axhline(beta_min, color='g')
plt.axhline(beta_max, color='g')
plt.xlabel('Modes')
plt.ylabel('Propagation constant')
plt.xlim(xlim)
plt.ylim(ylim)
plt.legend(fontsize=16, loc = 'upper right')
plt.title('High order modes')
[18]:
Text(0.5, 1.0, 'High order modes')
[19]:
i = 4
ex = M0_cs[:npoints**2,i].reshape([npoints]*2)
ey = M0_cs[npoints**2:2*npoints**2,i].reshape([npoints]*2)
ex2 = M0_rad[:npoints**2,i].reshape([npoints]*2)
ey2 = M0_rad[npoints**2:,i].reshape([npoints]*2)
m = np.max([np.max(np.abs(ex)),np.max(np.abs(ey))])
m2 = np.max([np.max(np.abs(ex2)),np.max(np.abs(ey2))])
plt.figure(figsize=(9,9))
plt.subplot(221)
plt.imshow(np.abs(ex), vmax = m)
plt.axis('off')
plt.title(r'$|E_x|$')
plt.subplot(222)
plt.imshow(np.abs(ey), vmax = m)
plt.axis('off')
plt.title(r'$|E_y|$')
plt.subplot(223)
plt.imshow(np.abs(ex2), vmax = m2)
plt.axis('off')
plt.subplot(224)
plt.imshow(np.abs(ey2), vmax = m2)
plt.axis('off')
[19]:
(-0.5, 127.5, 127.5, -0.5)
7. Projection of one basis onto the other one#
The order inside each degenerate group is not exacty the same, but it is normal as they are (quasi) degenerate.
First compute the mask corresponding to groups of degenerate modes#
[20]:
ind_degenerate = modes.getNearDegenerate(tol = 1e-3)
mask_degenerate = np.zeros(shape = [Nmodes*2]*2, dtype = float)
for ind in ind_degenerate:
mask_degenerate[2*min(ind):2*max(ind)+2,2*min(ind):2*max(ind)+2] = 1
# mask_degenerate = np.kron(np.eye(2),mask_degenerate)
plt.close('all')
plt.figure()
plt.imshow(mask_degenerate)
[20]:
<matplotlib.image.AxesImage at 0x788a92101cd0>
7.1 Display conversion matrices#
We project the mode matrices obtained with the different solver onto the modes profiles from FDTD simulations. Ideally, all energy is conserved and we have non-zero values only inside the degenerate block. Note that this matrix is not purely diagonal as modes are degenerate inside each block.
[21]:
# only keep transverse field and the propagating modes
M0_cs_trans = M0_cs[:2*npoints**2,:110]
A_as = M0_cs_trans.transpose().conjugate()@M0_rad
A_eig = M0_cs_trans.transpose().conjugate()@M0_eig
A_wkb = M0_cs_trans.transpose().conjugate()@M0_wkb
We define a loss term as the ratio energy of the conversion matrix that is out of the degenerate blocks.
[22]:
energy_loss_as = 1-np.linalg.norm(A_as*mask_degenerate)/np.linalg.norm(A_as)
energy_loss_eig = 1-np.linalg.norm(A_eig*mask_degenerate)/np.linalg.norm(A_eig)
energy_loss_wkb = 1-np.linalg.norm(A_wkb*mask_degenerate)/np.linalg.norm(A_wkb)
plt.figure(figsize = (10,4))
plt.subplot(131)
plt.imshow(np.abs(A_as), interpolation = 'None')
plt.axis('off')
plt.title(f'radial solver\n loss = {100*energy_loss_as:.2e}%')
plt.subplot(132)
plt.imshow(np.abs(A_eig), interpolation = 'None')
plt.axis('off')
plt.title(f'2D eigenvalue solver\n loss = {100*energy_loss_eig:.2e}%')
plt.subplot(133)
plt.imshow(np.abs(A_wkb), interpolation = 'None')
plt.axis('off')
plt.title(f'WKB\n loss = {100*energy_loss_wkb:.2e}%')
[22]:
Text(0.5, 1.0, 'WKB\n loss = 1.59e-03%')
7.2 Conversion efficiency per mode#
For each input mode found using the different technique, compute the energy that couples to the corresponding group of degenerate modes found using Comsol.
[23]:
Proj_as = np.sqrt(np.sum( np.abs(A_as*mask_degenerate)**2, axis = 1))
Proj_eig = np.sqrt(np.sum( np.abs(A_eig*mask_degenerate)**2, axis = 1))
Proj_wkb = np.sqrt(np.sum( np.abs(A_wkb*mask_degenerate)**2, axis = 1))
plt.figure()
plt.plot(100*Proj_as, label = 'Radial solver')
plt.plot(100*Proj_eig, label = '2D eigenvalue solver')
plt.plot(100*Proj_wkb, label = 'WKB')
plt.ylabel('Conversion efficiency')
plt.xlabel('Modes')
plt.legend()
yticks = np.linspace(95,100,6)
# plt.yticks(yticks, [f'{y:g}%' for y in yticks])
plt.title('Conversion efficiency onto the mode\n basis computed with Comsol')
plt.tight_layout()
7.3 Global conversion losses#
We compute the singular value decomposition of the conversion matrix between the modes found using Comsol and using the pyMMF solvers \(M_\text{cs}^\dagger . M_\text{pyMMF}\). Ideally, if two bases represent the same subspace, the conversion matrix has to be unitary, i.e. all its singular values are equal to one.
[27]:
_,s_as,_ = np.linalg.svd(A_as)
_,s_eig,_ = np.linalg.svd(A_eig)
_,s_wkb,_ = np.linalg.svd(A_wkb)
plt.figure(figsize = (10,6))
plt.plot(np.sort(np.abs(s_as)), label = 'axisymmetric solver')
plt.plot(np.sort(np.abs(s_eig)), label = '2D eigenvalue solver')
plt.plot(np.sort(np.abs(s_wkb)), label = 'WKB')
plt.ylim([0.9, 1.1])
plt.legend()
plt.title('Singular values of conversion matrices')
[27]:
Text(0.5, 1.0, 'Singular values of conversion matrices')
[ ]: