"""Module with parameters and functions for modeling of Random transfer matrix of the fiber""" # 4CCF: 4 coupled-core fiber import numpy as np import random import scipy.special as spec from scipy import signal import scipy.optimize as opt import matplotlib as mpl import matplotlib.pyplot as plt import matplotlib.animation as animation from IPython.display import display, HTML import sympy as sp from scipy.linalg import expm, qr, schur mpl.rcParams['agg.path.chunksize'] = 10000 from matplotlib.ticker import FormatStrFormatter SpeedofLight = 2.99792458e8 # speed of light [m/s] """BASIC FUNCTIONS""" def GaussianPulse(t_array, t0): Pulse = np.exp((-t_array ** 2)/(2 * t0 ** 2)) return Pulse def MatrixesMultiplication(matrixes_array): """Performs multiplication of all matrixes_array.shape[0] matrixes""" mult_result = matrixes_array[0] for i in range(1,matrixes_array.shape[0]): mult_result = np.matmul(mult_result, matrixes_array[i]) return mult_result def Sets_of_Matrixes(CouplingMatrixes, RandomPhaseMatrixes): """Performs multiplication of matrixes CouplingMatrixes and RandomPhaseMatrixes at different frequencies and flips the order of them""" Sets_of_Matrixes = np.empty((RandomPhaseMatrixes.shape[0], CouplingMatrixes.shape[0], CouplingMatrixes.shape[1], CouplingMatrixes.shape[1]), dtype=complex) for i in range(RandomPhaseMatrixes.shape[0]): total = np.empty((CouplingMatrixes.shape[0], CouplingMatrixes.shape[1], CouplingMatrixes.shape[1]), dtype=complex) for j in range(CouplingMatrixes.shape[0]): total[j] = np.matmul(CouplingMatrixes[j], RandomPhaseMatrixes[i]) Sets_of_Matrixes[i] = total Sets_of_Matrixes_reverse = np.flip(Sets_of_Matrixes, 0) return Sets_of_Matrixes_reverse, Sets_of_Matrixes """COUPLING COEFFICIENTS, COUPLING MATRIXES""" def ParametersCalculation_UW(NormalizedFrequency): var1 = [] for i in range(NormalizedFrequency.shape[0]): def f(variables): (U, W) = variables first_eq = U**2 + W**2 - NormalizedFrequency[i]**2 second_eq = U*spec.k0(W)*spec.j1(U)-W*spec.k1(W)*spec.j0(U) return [first_eq, second_eq] solution = opt.fsolve(f, (0.1, 1)) var1.append(solution) Solution = np.array(var1) return Solution def CouplingCoefficients(core_radii, core_pitch, frequencies, index_diff): SpeedOfLight = 3e8 wave_vector = 2*np.pi*frequencies/SpeedOfLight lambda_sq = (SpeedOfLight*1e6/frequencies)**2 # Sellmeier equation for n [Saleh, Fundamentals of Fotonics, p.180] n1_sq = 1+(0.6962*lambda_sq)/(lambda_sq - 0.0684**2) + (0.4079*lambda_sq)/(lambda_sq - 0.1162**2)+ (0.8975*lambda_sq)/(lambda_sq - 9.8962**2) n2_sq = (np.sqrt(n1_sq)*(1-index_diff))**2 V = 2*np.pi*frequencies*core_radii*np.sqrt(n1_sq-n2_sq)/SpeedOfLight U = ParametersCalculation_UW(V)[:,0] W = ParametersCalculation_UW(V)[:,1] W_ = W*core_pitch/core_radii # argument of Ko PropagationConstant = np.sqrt(n1_sq*wave_vector**2-(U/core_radii)**2) coupling_coef_function = np.sqrt((n1_sq-n2_sq)/n1_sq)*U**2*spec.k0(W_)/(core_radii*V**3*spec.k1(W)**2) #d_freq = (frequencies[1]-frequencies[0]) d_freq = 2*np.pi*(frequencies[1]-frequencies[0]) dC_dw = np.gradient(coupling_coef_function, d_freq) dBeta_dw = np.gradient(PropagationConstant, d_freq) birefringence = 2*coupling_coef_function*index_diff*W*spec.kn(2, W_)*(spec.iv(1, W)*spec.k0(W)-spec.iv(2, W)*spec.k1(W))/spec.k0(W_) return coupling_coef_function, dC_dw, PropagationConstant, birefringence def RefrIndex(frequencies, index_diff): SpeedOfLight = 3e8 wave_vector = 2*np.pi*frequencies/SpeedOfLight lambda_sq = (SpeedOfLight*1e6/frequencies)**2 # Sellmeier equation for n [Saleh, Fundamentals of Fotonics, p.180] n1_sq = 1+(0.6962*lambda_sq)/(lambda_sq - 0.0684**2) + (0.4079*lambda_sq)/(lambda_sq - 0.1162**2)+ (0.8975*lambda_sq)/(lambda_sq - 9.8962**2) n2_sq = (np.sqrt(n1_sq)*(1-index_diff))**2 return np.sqrt(n1_sq), np.sqrt(n2_sq) def CouplingMatrixIdeal_scalar(PropagationConstant, coupl_coeff, num_cores): matrix_CoupCoeff = np.full((num_cores, num_cores), coupl_coeff) PropConstans = np.full(num_cores, PropagationConstant) matrix_PropConstans = np.diag(PropConstans-coupl_coeff) CouplingMatrix = matrix_CoupCoeff+matrix_PropConstans return CouplingMatrix def CouplingMatrixScalar_4ccf(PropagationConstant, coupl_coeff1, coupl_coeff2): CouplingMatrix = np.array([[PropagationConstant, coupl_coeff1, coupl_coeff1, coupl_coeff2], [coupl_coeff1, PropagationConstant, coupl_coeff2, coupl_coeff1], [coupl_coeff1, coupl_coeff2, PropagationConstant, coupl_coeff1], [coupl_coeff2, coupl_coeff1, coupl_coeff1, PropagationConstant]]) return CouplingMatrix def CouplingMatrix8x8(alpha1, alpha2, b1, b2, c1, c2): Rotation_matrix1 = np.array([[np.cos(alpha1), np.sin(alpha1)], [-np.sin(alpha1), np.cos(alpha1)]]) Rotation_matrix2 = np.array([[np.cos(alpha2), np.sin(alpha2)], [-np.sin(alpha2), np.cos(alpha2)]]) Birefringence_matrix1 = np.array([[-b1, 0], [0, b1]]) Birefringence_matrix2 = np.array([[-b2, 0], [0, b2]]) Rotation_matrix1_inverse = np.linalg.inv(Rotation_matrix1) Rotation_matrix2_inverse = np.linalg.inv(Rotation_matrix2) core1 = Birefringence_matrix1 + Rotation_matrix1_inverse@Birefringence_matrix1@Rotation_matrix1 + Rotation_matrix2@Birefringence_matrix2@Rotation_matrix2_inverse core2 = Rotation_matrix1@Birefringence_matrix1@Rotation_matrix1_inverse + Birefringence_matrix1 + Rotation_matrix2_inverse@Birefringence_matrix2@Rotation_matrix2 core3 = Rotation_matrix1_inverse@Birefringence_matrix1@Rotation_matrix1 + Birefringence_matrix1 + Rotation_matrix2_inverse@Birefringence_matrix2@Rotation_matrix2 core4 = Rotation_matrix1@Birefringence_matrix1@Rotation_matrix1_inverse + Birefringence_matrix1 + Rotation_matrix2@Birefringence_matrix2@Rotation_matrix2_inverse A_Coupling = np.array([[core1[0][0], core1[0][1], c1, 0, c1, 0, c2, 0], [core1[1][0], core1[1][1], 0 , c1, 0, c1, 0, c2], [c1, 0, core2[0][0], core2[0][1], c2, 0, c1, 0], [0, c1, core2[1][0], core2[1][1], 0, c2, 0, c1], [c1, 0, c2, 0, core3[0][0], core3[0][1], c1, 0], [0, c1, 0, c2, core3[1][0], core3[1][1], 0, c1], [c2, 0, c1, 0, c1, 0, core4[1][0], core4[1][1]], [0, c2, 0, c1, 0, c1, core4[1][0], core4[1][1]]]) return A_Coupling """GROUP DELAYS OF THE SUPERMODES""" def GDs_supermodes_4ccf(num_cores, core_radii, core_pitch, frequencies, index_diff, alpha1, alpha2): coupl_coeff1, dC_dw1, PropagationConstant1, birefringence1 = CouplingCoefficients(core_radii, core_pitch, frequencies, index_diff) coupl_coeff2, dC_dw2, PropagationConstant2, birefringence2 = CouplingCoefficients(core_radii, np.sqrt(2)*core_pitch, frequencies, index_diff) Coup_matrix = np.empty((frequencies.shape[0], num_cores, num_cores)) Coup_matrix_NoCoup = np.empty((frequencies.shape[0], 2*num_cores, 2*num_cores)) Coup_matrix_vec = np.empty((frequencies.shape[0], 2*num_cores, 2*num_cores)) for i in range(frequencies.shape[0]): Coup_matrix[i] = CouplingMatrixScalar_4ccf(0, coupl_coeff1[i], coupl_coeff2[i]) Coup_matrix_NoCoup[i] = CouplingMatrix8x8(alpha1, alpha2, birefringence1[i], birefringence2[i], 0, 0) Coup_matrix_vec[i] = CouplingMatrix8x8(alpha1, alpha2, birefringence1[i], birefringence2[i], coupl_coeff1[i], coupl_coeff2[i]) Eigenvalues_scalar, Eigenvectors_scalar = np.linalg.eigh(Coup_matrix) Eigenvalues_NoCoup, Eigenvectors_NoCoup = np.linalg.eigh(Coup_matrix_NoCoup) Eigenvalues_vector, Eigenvectors_vector = np.linalg.eigh(Coup_matrix_vec) dw = 2*np.pi*(frequencies[1]-frequencies[0]) GDs_scalar = np.empty((num_cores, frequencies.shape[0])) GDs_NoCoup = np.empty((2*num_cores, frequencies.shape[0])) GDs_vector = np.empty((2*num_cores, frequencies.shape[0])) for j in range(num_cores): GDs_scalar[j] = np.gradient(Eigenvalues_scalar[:, j], dw) for k in range(2*num_cores): GDs_NoCoup[k] = np.gradient(Eigenvalues_NoCoup[:, k], dw) GDs_vector[k] = np.gradient(Eigenvalues_vector[:, k], dw) return Eigenvalues_scalar, Eigenvectors_scalar, GDs_scalar, Eigenvalues_NoCoup, Eigenvectors_NoCoup, GDs_NoCoup, Eigenvalues_vector, Eigenvectors_vector, GDs_vector """RANDOM MATRIXES""" def RandomPhaseMatrix(num_cores, num_polarisations): """Computes the matrix with random phase shift in every core""" ra = [] for i in range(num_cores): randomPhaseShift = random.uniform(0, 2*np.pi) for j in range(num_polarisations): ra.append(randomPhaseShift) RandomPhaseShift = np.array(ra) RandomPhaseMatrix = np.diag(np.exp(1j*RandomPhaseShift)) return RandomPhaseMatrix, RandomPhaseShift def RandomUnitaryMatrix(dimention): """Creates the random unitary matrix in every core following the QR decomposition: https: // nhigham.com / 2020 / 04 / 22 / what - is -a - random - orthogonal - matrix /""" matr = np.empty((dimention,dimention),dtype=complex) for ll in range(dimention): for lll in range(dimention): # real_part = random.uniform(0,2*np.pi) # imag_part = random.uniform(0,2*np.pi) real_part = np.random.normal() imag_part = np.random.normal() matr[ll,lll] = complex(real_part,imag_part) Q1, R = qr(matr) ComplexSign = np.diag(R)/abs(np.diag(R)) #Q = np.matmul(Q1, np.diag(np.sign(np.diag(R)))) Q = np.matmul(Q1, np.diag(ComplexSign)) return Q """TOTAL RANDOM TRANSFER MATRIX""" def TotalTransferMatrix_4ccf_Scalar(CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, RandomMatrix, a): #a = int(input("Use a constant fiber length - 1, Use a constant DGD element length - 2")) SpeedOfLight = 3e8 f2 = SpeedOfLight/wl_list[0] f1 = SpeedOfLight/wl_list[1] frequencies = np.linspace(f1, f2, freq_points) coupl_coeff1, dC_dw1, PropagationConstant1, birefringence1 = CouplingCoefficients(CoreRadii, CorePitch, frequencies, indDiff) coupl_coeff2, dC_dw2, PropagationConstant2, birefringence2 = CouplingCoefficients(CoreRadii, np.sqrt(2)*CorePitch, frequencies, indDiff) if a==1: distance = fiberlength/num_pieces else: distance = DGD_elementLength M_coup = CouplingMatrixScalar_4ccf(0, coupl_coeff1[0], coupl_coeff2[0]) Coupl_matrix = np.empty((num_pieces, frequencies.shape[0], M_coup.shape[0], M_coup.shape[0]), dtype=complex) for i in range(num_pieces): M_cp = np.empty((freq_points, M_coup.shape[0], M_coup.shape[0]), dtype=complex) for s in range(freq_points): dc_dw11 = dC_dw1[0]*2*np.pi*frequencies[s] dc_dw22 = dC_dw2[0]*2*np.pi*frequencies[s] m_d = np.array([[0, dc_dw11, dc_dw11, dc_dw22], [dc_dw11, 0, dc_dw22, dc_dw11], [dc_dw11, dc_dw22, 0, dc_dw11], [dc_dw22, dc_dw11, dc_dw11, 0]]) m_cp = M_coup+m_d m_cp_z = m_cp*distance m_cp_z_exp = expm(1j*m_cp_z) M_cp[s] = m_cp_z_exp M_cp[s] = np.matmul(m_cp_z_exp, RandomMatrix[i]) # M_cp1 = np.matmul(RandomMatrix[i].conjugate(), m_cp_z_exp) # M_cp[s] = np.matmul(M_cp1, RandomMatrix[i]) Coupl_matrix[i] = M_cp M_total = MatrixesMultiplication(np.flip(Coupl_matrix, 0)) return M_total, frequencies def TotalTransferMatrix_4ccf_Vector(CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, RandomMatrix, rotAngle1, rotAngle2, a): #a = int(input("Use a constant fiber length - 1, Use a constant DGD element length - 2")) SpeedOfLight = 3e8 f2 = SpeedOfLight/wl_list[0] f1 = SpeedOfLight/wl_list[1] frequencies = np.linspace(f1, f2, freq_points) coupl_coeff1, dC_dw1, PropagationConstant1, birefringence1 = CouplingCoefficients(CoreRadii, CorePitch, frequencies, indDiff) coupl_coeff2, dC_dw2, PropagationConstant2, birefringence2 = CouplingCoefficients(CoreRadii, np.sqrt(2)*CorePitch, frequencies, indDiff) if a==1: distance = fiberlength/num_pieces else: distance = DGD_elementLength M_coup = CouplingMatrix8x8(rotAngle1, rotAngle2, birefringence1[0], birefringence2[0], coupl_coeff1[0], coupl_coeff2[0]) Coupl_matrix = np.empty((num_pieces, frequencies.shape[0], M_coup.shape[0], M_coup.shape[0]), dtype=complex) for i in range(num_pieces): M_cp = np.empty((freq_points, M_coup.shape[0], M_coup.shape[0]), dtype=complex) for s in range(freq_points): dc_dw11 = dC_dw1[0]*2*np.pi*frequencies[s] dc_dw22 = dC_dw2[0]*2*np.pi*frequencies[s] # np.array([[core1[0][0], core1[0][1], c1, 0, c1, 0, c2, 0], # [core1[1][0], core1[1][1], 0, c1, 0, c1, 0, c2], # [c1, 0, core2[0][0], core2[0][1], c2, 0, c1, 0], # [0, c1, core2[1][0], core2[1][1], 0, c2, 0, c1], # [c1, 0, c2, 0, core3[0][0], core3[0][1], c1, 0], # [0, c1, 0, c2, core3[1][0], core3[1][1], 0, c1], # [c2, 0, c1, 0, c1, 0, core4[1][0], core4[1][1]], # [0, c2, 0, c1, 0, c1, core4[1][0], core4[1][1]]]) m_d = np.array([[0, 0, dc_dw11, 0, dc_dw11, 0, dc_dw22, 0], [0, 0, 0, dc_dw11, 0, dc_dw11, 0, dc_dw22], [dc_dw11, 0, 0, 0, dc_dw22, 0, dc_dw11, 0], [0, dc_dw11, 0, 0, 0, dc_dw22, 0, dc_dw11], [dc_dw11, 0, dc_dw22, 0, 0, 0, dc_dw11, 0], [0, dc_dw11, 0, dc_dw22, 0, 0, 0, dc_dw11], [dc_dw22, 0, dc_dw11, 0, dc_dw11, 0, 0, 0], [0, dc_dw22, 0, dc_dw11, 0, dc_dw11, 0, 0]]) m_cp = M_coup+m_d m_cp_z = m_cp*distance m_cp_z_exp = expm(1j*m_cp_z) M_cp[s] = m_cp_z_exp M_cp[s] = np.matmul(m_cp_z_exp, RandomMatrix[i]) # M_cp1 = np.matmul(RandomMatrix[i].conjugate(), m_cp_z_exp) # M_cp[s] = np.matmul(M_cp1, RandomMatrix[i]) Coupl_matrix[i] = M_cp M_total = MatrixesMultiplication(np.flip(Coupl_matrix, 0)) return M_total, frequencies def TransferFunction4ccfScalar_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, a): """Calculates the total transfer matrix in the scalar case (no polarization diversity)""" RandomMatrix = np.empty((num_pieces, num_cores, num_cores), dtype=complex) for i in range(num_pieces): RandomMatrix[i] = RandomUnitaryMatrix(num_cores) M_tot, freqqs = TotalTransferMatrix_4ccf_Scalar(CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces, freq_points, wl_list, RandomMatrix, a) return M_tot, freqqs def TransferFunction4ccfVector_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, a): """Calculates the total transfer matrix in the scalar case (no polarization diversity)""" RandomMatrix = np.empty((num_pieces, 2*num_cores, 2*num_cores), dtype=complex) for i in range(num_pieces): RandomMatrix[i] = RandomUnitaryMatrix(2*num_cores) M_tot, freqqs = TotalTransferMatrix_4ccf_Vector(CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces, freq_points, wl_list, RandomMatrix, rotAngle1, rotAngle2, a) return M_tot, freqqs def TransferFunction4ccfScalar_averagePerformance(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, num_scenarious, a): """Averages total transfer matrix over different scenarious""" ab = np.empty((num_scenarious, freq_points, num_cores, num_cores), dtype=complex) for l in range(num_scenarious): ab[l] = TransferFunction4ccfScalar_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, a)[0] M_tots_av = ab.mean(axis = 0) return M_tots_av def TransferFunction4ccfVector_averagePerformance(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, num_scenarious, a): """Averages total transfer matrix over different scenarious""" ab = np.empty((num_scenarious, freq_points, 2*num_cores, 2*num_cores), dtype=complex) for l in range(num_scenarious): ab[l] = TransferFunction4ccfVector_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, a)[0] M_tots_av = ab.mean(axis = 0) return M_tots_av """PROPAGATION""" def PropagationSimple(TrMatrix, Amplitude): Amplitude_out = np.empty((TrMatrix.shape[0], Amplitude.shape[0]), dtype=complex) for i in range(TrMatrix.shape[0]): Amplitude_out[i] = np.matmul(TrMatrix[i], Amplitude) return Amplitude_out def Propagation_GaussianPulse(TrFunc, Amp0, freq_array, T0): T_array = np.fft.fftshift(np.fft.fftfreq(freq_array.shape[0], freq_array[1] - freq_array[0])) Pulse = GaussianPulse(T_array, T0) Pulse_fft = np.fft.fftshift(np.fft.fft(Pulse)) Amplitude_in = np.empty((TrFunc.shape[0], TrFunc.shape[-1]), dtype=complex) Amplitude_out = np.empty((TrFunc.shape[0], TrFunc.shape[-1]), dtype=complex) IR = np.empty((TrFunc.shape[-1], TrFunc.shape[0])) for i in range(TrFunc.shape[0]): #Amplitude_in[i] = Amp0*Pulse[i] Amplitude_in[i] = Amp0*Pulse_fft[i] for k in range(TrFunc.shape[0]): Amplitude_out[k] = np.matmul(TrFunc[k], Amplitude_in[k]) for l in range(TrFunc.shape[-1]): IR[l] = abs(np.fft.fftshift(np.fft.fft(Amplitude_out[:, l])))**2 return Amplitude_out, np.fft.fftshift(IR), T_array def Propagation_GaussianPulse_4ccfVector(num_cores, CoreRadii, CorePitch, indDiff, num_pieces, wl_list, freq_points, fiberlength, DGD_elementLength, rotAngle1, rotAngle2, Amp0, T0,a): TrFunc, freq_array = TransferFunction4ccfVector_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, a) T_array = np.fft.fftshift(np.fft.fftfreq(freq_array.shape[0], freq_array[1] - freq_array[0])) Pulse = GaussianPulse(T_array, T0) Pulse_fft = np.fft.fftshift(np.fft.fft(Pulse)) Amplitude_in = np.empty((TrFunc.shape[0], TrFunc.shape[-1]), dtype=complex) Amplitude_out = np.empty((TrFunc.shape[0], TrFunc.shape[-1]), dtype=complex) IR = np.empty((TrFunc.shape[-1], TrFunc.shape[0])) for i in range(TrFunc.shape[0]): Amplitude_in[i] = Amp0*Pulse_fft[i] for k in range(TrFunc.shape[0]): Amplitude_out[k] = np.matmul(TrFunc[k], Amplitude_in[k]) for l in range(TrFunc.shape[-1]): IR[l] = abs(np.fft.fftshift(np.fft.fft(Amplitude_out[:, l])))**2 return Amplitude_out, np.fft.fftshift(IR), T_array """IMPULSE RESPONSE""" def Impulse_response(TransferFunction, frequencies): """Calculates impulse response of every element of the transfer function""" ImpulseResponse = np.empty((TransferFunction.shape[1],TransferFunction.shape[1],TransferFunction.shape[0])) for y in range(TransferFunction[1].shape[0]): for x in range(TransferFunction[1].shape[0]): ImpulseResponse[y,x] = abs(np.fft.fftshift(np.fft.fft((TransferFunction[:, y, x]))))**2 time_array = np.fft.fftshift(np.fft.fftfreq(ImpulseResponse.shape[-1], (frequencies[1]-frequencies[0]))) return ImpulseResponse, time_array def ImpulseResponse_av_GaussianPulse_scalar(num_cores, CoreRadii, CorePitch, indDiff, fiberlength,num_pieces,freq_points,wl_list, num_scenarious, Amp0, T0): Amplitude_in = np.empty((freq_points, num_cores), dtype=complex) Amplitude_out = np.empty((num_scenarious, freq_points, num_cores), dtype=complex) SpeedOfLight = 3e8 freq2 = SpeedOfLight/wl_list[0] freq1 = SpeedOfLight/wl_list[1] freqs = np.linspace(freq1, freq2, freq_points) #freqs = TransferFunctionScalar_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength,num_pieces,freq_points,wl_list)[1] #T_array = np.fft.fftshift(np.fft.fftfreq(freqs.shape[0], 2*np.pi*(freqs[1] - freqs[0]))) T_array = np.fft.fftshift(np.fft.fftfreq(freqs.shape[0], freqs[1] - freqs[0])) Pulse = GaussianPulse(T_array, T0) Pulse_fft = np.fft.fftshift(np.fft.fft(Pulse)) for i in range(freq_points): #Amplitude_in[i] = Amp0*Pulse[i] Amplitude_in[i] = Amp0*Pulse_fft[i] for j in range(num_scenarious): TrFunc = TransferFunctionScalar_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength,num_pieces,freq_points,wl_list)[0] for k in range(freq_points): Amplitude_out[j][k] = np.matmul(TrFunc[k], Amplitude_in[k]) IR = np.empty((num_scenarious, num_cores, freq_points)) for m in range(num_scenarious): for l in range(num_cores): IR[m][l] = abs(np.fft.fftshift(np.fft.fft(Amplitude_out[m, :, l])))**2 IR_total = np.empty((num_scenarious, freq_points)) for t in range(num_scenarious): IR_total[t] = np.fft.fftshift(np.sum(IR[t], axis=0)) IR_total_av = IR_total.mean(axis=0) return Amplitude_out, IR, IR_total, IR_total_av, T_array def ImpulseResponse_av_GaussianPulse_4CCFscalar(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength,num_pieces,freq_points,wl_list, num_scenarious, Amp0, T0, a): Amplitude_in = np.empty((freq_points, num_cores), dtype=complex) Amplitude_out = np.empty((num_scenarious, freq_points, num_cores), dtype=complex) SpeedOfLight = 3e8 freq2 = SpeedOfLight/wl_list[0] freq1 = SpeedOfLight/wl_list[1] freqs = np.linspace(freq1, freq2, freq_points) T_array = np.fft.fftshift(np.fft.fftfreq(freqs.shape[0], freqs[1] - freqs[0])) Pulse = GaussianPulse(T_array, T0) Pulse_fft = np.fft.fftshift(np.fft.fft(Pulse)) for i in range(freq_points): Amplitude_in[i] = Amp0*Pulse_fft[i] for j in range(num_scenarious): TrFunc = TransferFunction4ccfScalar_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, a)[0] for k in range(freq_points): Amplitude_out[j][k] = np.matmul(TrFunc[k], Amplitude_in[k]) IR = np.empty((num_scenarious, num_cores, freq_points)) for m in range(num_scenarious): for l in range(num_cores): IR[m][l] = abs(np.fft.fftshift(np.fft.fft(Amplitude_out[m, :, l])))**2 IR_total = np.empty((num_scenarious, freq_points)) for t in range(num_scenarious): IR_total[t] = np.fft.fftshift(np.sum(IR[t], axis=0)) IR_total_av = IR_total.mean(axis=0) return Amplitude_out, IR, IR_total, IR_total_av, T_array def ImpulseResponse_av_GaussianPulse_4CCFvector(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength,num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, num_scenarious, Amp0, T0, a): Amplitude_in = np.empty((freq_points, num_cores*2), dtype=complex) Amplitude_out = np.empty((num_scenarious, freq_points, num_cores*2), dtype=complex) freqs = TransferFunction4ccfVector_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, a)[1] T_array = np.fft.fftshift(np.fft.fftfreq(freqs.shape[0], freqs[1] - freqs[0])) Pulse = GaussianPulse(T_array, T0) Pulse_fft = np.fft.fftshift(np.fft.fft(Pulse)) for i in range(freq_points): Amplitude_in[i] = Amp0*Pulse_fft[i] for j in range(num_scenarious): TrFunc = TransferFunction4ccfVector_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, a)[0] for k in range(freq_points): Amplitude_out[j][k] = np.matmul(TrFunc[k], Amplitude_in[k]) IR = np.empty((num_scenarious, num_cores*2, freq_points)) for m in range(num_scenarious): for l in range(6): IR[m][l] = abs(np.fft.fftshift(np.fft.fft(Amplitude_out[m, :, l])))**2 IR_total = np.empty((num_scenarious, freq_points)) for t in range(num_scenarious): IR_total[t] = np.fft.fftshift(np.sum(IR[t],axis=0)) IR_total_av = IR_total.mean(axis = 0) return Amplitude_out, IR, IR_total, IR_total_av, T_array """DELAY OPERATOR, DELAYS""" def DelayOperator(totalTransferMatrix,frequencies): """Calculates the Delay operator using gradient function""" mm_ = [] dw = 2*np.pi*(frequencies[1]-frequencies[0]) M_diff = np.gradient(totalTransferMatrix, dw, axis=0, edge_order=2) for s in range(M_diff.shape[0]): M_ConjTransp = (totalTransferMatrix[s].conjugate()).transpose() delayOperator = 1j*np.matmul(M_ConjTransp, M_diff[s]) #delayOperator = np.real_if_close(delayOperator1, tol=100) #delayOperator = np.around(delayOperator11, 14) mm_.append(delayOperator) DelayOperator = np.array(mm_) return DelayOperator, M_diff def Delays(DelayOperator): Eigenvalues1, Eigenvectors = np.linalg.eigh(DelayOperator[1:-1,:,:]) Eigenvalues = np.real_if_close(Eigenvalues1, tol=10000) #1000 2cf return Eigenvalues def Delays_Compute(M_tot, freqs): DelOp = DelayOperator(M_tot, freqs)[0] Eigenvalues1, Eigenvectors = np.linalg.eigh(DelOp[1:-1, :, :]) Eigenvalues = np.real_if_close(Eigenvalues1, tol=10000) #1000 2cf return Eigenvalues def DGD_Compute(M_tot, freqs): DelOp = DelayOperator(M_tot, freqs)[0] Eigenvalues1, Eigenvectors = np.linalg.eigh(DelOp[1:-1, :, :]) Eigenvalues = np.real_if_close(Eigenvalues1, tol=10000) #DGD = Eigenvalues.max()-Eigenvalues.min() DGD = np.sqrt(np.sum(Eigenvalues**2))/M_tot.shape[1] return DGD def DGD_theory(num_pieces, piece_length, GDs_per_m): #DGD_sq = num_pieces*(DGD_per_m*piece_length)**2 eigenvalues = GDs_per_m*piece_length DGD_sq = num_pieces*np.sum(eigenvalues**2)/GDs_per_m.shape[0] return np.sqrt(DGD_sq) def DGD_caculationScalar_4ccf(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, num_scenarious, a): freqs = TransferFunction4ccfScalar_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, a)[1] Eigenvalues = np.empty((num_scenarious, num_cores, freq_points-2)) for l in range(num_scenarious): M_tot = TransferFunction4ccfScalar_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, a)[0] Eigenvalues_ = Delays_Compute(M_tot, freqs) for ll in range(num_cores): Eigenvalues[l, ll] = Eigenvalues_[:, ll] EigenvaluesSq_av = np.empty((num_cores, 1)) for i in range(num_cores): EigenvaluesSq_av[i] = np.mean(Eigenvalues[:, i, int(freq_points/2)]**2) DGD1 = np.sqrt(np.sum(EigenvaluesSq_av)/M_tot.shape[1]) Eigenvalues_dif_Sq = np.empty((num_scenarious, 1)) for j in range(num_scenarious): Eigenvalues_dif_Sq[j] = (Eigenvalues[j, :, int(freq_points/2)].max() - Eigenvalues[j, :, int(freq_points/2)].min())**2 DGD2 = np.sqrt(np.mean(Eigenvalues_dif_Sq)) return DGD2, DGD1 def DGD_caculationVector_4ccf(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, num_scenarious, a): freqs = TransferFunction4ccfVector_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, a)[1] Eigenvalues = np.empty((num_scenarious, num_cores*2, freq_points-2)) for l in range(num_scenarious): M_tot = TransferFunction4ccfVector_Compute(num_cores, CoreRadii, CorePitch, indDiff, fiberlength, DGD_elementLength, num_pieces,freq_points,wl_list, rotAngle1, rotAngle2, a)[0] Eigenvalues_ = Delays_Compute(M_tot, freqs) for ll in range(num_cores*2): Eigenvalues[l, ll] = Eigenvalues_[:, ll] EigenvaluesSq_av = np.empty((num_cores*2, 1)) for i in range(num_cores*2): EigenvaluesSq_av[i] = np.mean(Eigenvalues[:, i, int(freq_points / 2)]**2) DGD1 = np.sqrt(np.sum(EigenvaluesSq_av)/M_tot.shape[1]) Eigenvalues_dif_Sq = np.empty((num_scenarious, 1)) for j in range(num_scenarious): Eigenvalues_dif_Sq[j] = (Eigenvalues[j, :, int(freq_points/2)].max() - Eigenvalues[j, :, int(freq_points/2)].min())** 2 DGD2 = np.sqrt(np.mean(Eigenvalues_dif_Sq)) return DGD2, DGD1 """PLOTTING FUNCTIONS""" def Supermodes_4ccf_scalar(size_a, size_b): """For plotting of all 4 eigenvectors of 4CCF""" fig, ax = plt.subplots(nrows=2, ncols=8, figsize=(size_a, size_b)) theta = np.linspace(0, 2*np.pi, 100) r = np.sqrt(1.2) x1 = r*np.cos(theta) x2 = r*np.sin(theta) limit1 = 1.2 fs = 12 coord_x = -0.38 coord_y = -0.38 ax[0,0].plot(x1, x2) ax[0,0].text(coord_x, coord_y, "-", fontsize=16) ax[0,0].set_ylim(-limit1, limit1) ax[0,0].set_xlim(-limit1, limit1) ax[0,0].tick_params(axis='y', labelsize=5) ax[0,0].tick_params(axis='x', labelsize=5) ax[0,0].axis('off') ax[0,1].plot(x1, x2) ax[0,1].text(coord_x, coord_y, "0", fontsize=fs) ax[0,1].set_ylim(-limit1, limit1) ax[0,1].set_xlim(-limit1, limit1) ax[0,1].tick_params(left = False, right = False , labelleft = False , bottom = True, labelbottom = True) ax[0,1].tick_params(axis='y', labelsize=5) ax[0,1].tick_params(axis='x', labelsize=5) ax[0,1].axis('off') ax[0,2].plot(x1,x2, color="red") ax[0,2].text(coord_x,coord_y, "0", fontsize=fs) ax[0,2].set_ylim(-limit1, limit1) ax[0,2].set_xlim(-limit1, limit1) ax[0,2].tick_params(left = False, right = False , labelleft = False , bottom = True, labelbottom = True) ax[0,2].tick_params(axis='y', labelsize=5) ax[0,2].tick_params(axis='x', labelsize=5) ax[0,2].axis('off') ax[0,3].plot(x1,x2, color="red") ax[0,3].text(coord_x,coord_y, "-", fontsize=16) ax[0,3].set_ylim(-limit1, limit1) ax[0,3].set_xlim(-limit1, limit1) ax[0,3].tick_params(left = False, right = False , labelleft = False , bottom = True, labelbottom = True) ax[0,3].tick_params(axis='y', labelsize=5) ax[0,3].tick_params(axis='x', labelsize=5) ax[0,3].axis('off') ax[0,4].plot(x1,x2, color="green") ax[0,4].text(coord_x,coord_y, "+", fontsize=fs) ax[0,4].set_ylim(-limit1, limit1) ax[0,4].set_xlim(-limit1, limit1) ax[0,4].tick_params(left = False, right = False , labelleft = False , bottom = True, labelbottom = True) ax[0,4].tick_params(axis='y', labelsize=5) ax[0,4].tick_params(axis='x', labelsize=5) ax[0,4].axis('off') ax[0,5].plot(x1,x2, color="green") ax[0,5].text(coord_x,coord_y, "-", fontsize=16) ax[0,5].set_ylim(-limit1, limit1) ax[0,5].set_xlim(-limit1, limit1) ax[0,5].tick_params(left = False, right = False , labelleft = False , bottom = True, labelbottom = True) ax[0,5].tick_params(axis='y', labelsize=5) ax[0,5].tick_params(axis='x', labelsize=5) ax[0,5].axis('off') ax[0,6].plot(x1,x2, color="orange") ax[0,6].text(coord_x,coord_y, "+", fontsize=fs) ax[0,6].set_ylim(-limit1, limit1) ax[0,6].set_xlim(-limit1, limit1) ax[0,6].tick_params(left = False, right = False , labelleft = False , bottom = True, labelbottom = True) ax[0,6].tick_params(axis='y', labelsize=5) ax[0,6].tick_params(axis='x', labelsize=5) ax[0,6].axis('off') ax[0,7].plot(x1,x2, color="orange") ax[0,7].text(coord_x,coord_y, "+", fontsize=fs) ax[0,7].set_ylim(-limit1, limit1) ax[0,7].set_xlim(-limit1, limit1) ax[0,7].tick_params(left = False, right = False , labelleft = False , bottom = True, labelbottom = True) ax[0,7].tick_params(axis='y', labelsize=5) ax[0,7].tick_params(axis='x', labelsize=5) ax[0,7].axis('off') ax[1,0].plot(x1,x2) ax[1,0].text(coord_x,coord_y, "o", fontsize=fs) ax[1,0].set_ylim(-limit1, limit1) ax[1,0].set_xlim(-limit1, limit1) ax[1,0].tick_params(axis='y', labelsize=5) ax[1,0].tick_params(axis='x', labelsize=5) ax[1,0].axis('off') ax[1,1].plot(x1,x2) ax[1,1].text(coord_x,coord_y, "+", fontsize=fs) ax[1,1].set_ylim(-limit1, limit1) ax[1,1].set_xlim(-limit1, limit1) ax[1,1].tick_params(axis='y', labelsize=5) ax[1,1].tick_params(axis='x', labelsize=5) ax[1,1].axis('off') ax[1,2].plot(x1,x2, color="red") ax[1,2].text(coord_x,coord_y, "+", fontsize=fs) ax[1,2].set_ylim(-limit1, limit1) ax[1,2].set_xlim(-limit1, limit1) ax[1,2].tick_params(axis='y', labelsize=5) ax[1,2].tick_params(axis='x', labelsize=5) ax[1,2].axis('off') ax[1,3].plot(x1,x2, color="red") ax[1,3].text(coord_x,coord_y, "o", fontsize=fs) ax[1,3].set_ylim(-limit1, limit1) ax[1,3].set_xlim(-limit1, limit1) ax[1,3].tick_params(axis='y', labelsize=5) ax[1,3].tick_params(axis='x', labelsize=5) ax[1,3].axis('off') ax[1,4].plot(x1,x2, color="green") ax[1,4].text(coord_x,coord_y, "-", fontsize=16) ax[1,4].set_ylim(-limit1, limit1) ax[1,4].set_xlim(-limit1, limit1) ax[1,4].tick_params(axis='y', labelsize=5) ax[1,4].tick_params(axis='x', labelsize=5) ax[1,4].axis('off') ax[1,5].plot(x1,x2, color="green") ax[1,5].text(coord_x,coord_y, "+", fontsize=fs) ax[1,5].set_ylim(-limit1, limit1) ax[1,5].set_xlim(-limit1, limit1) ax[1,5].tick_params(axis='y', labelsize=5) ax[1,5].tick_params(axis='x', labelsize=5) ax[1,5].axis('off') ax[1,6].plot(x1, x2, color="orange") ax[1,6].text(coord_x, coord_y, "+", fontsize=fs) ax[1,6].set_ylim(-limit1, limit1) ax[1,6].set_xlim(-limit1, limit1) ax[1,6].tick_params(axis='y', labelsize=5) ax[1,6].tick_params(axis='x', labelsize=5) ax[1,6].axis('off') ax[1,7].plot(x1, x2, color="orange") ax[1,7].text(coord_x, coord_y, "+", fontsize=fs) ax[1,7].set_ylim(-limit1, limit1) ax[1,7].set_xlim(-limit1, limit1) ax[1,7].tick_params(axis='y', labelsize=5) ax[1,7].tick_params(axis='x', labelsize=5) ax[1,7].axis('off') return def Supermodes_4ccf_vector(EigVect, size_a, size_b): """For plotting of all 8 eigenvectors of 4CCF""" firstrow = np.empty((8, 4)) secondrow = np.empty((8, 4)) for j in range(0, 8): firstrow[j] = EigVect[j][:4] secondrow[j] = EigVect[j][4:8] fig, ax = plt.subplots(nrows=2, ncols=16, figsize=(size_a, size_b)) theta = np.linspace(0, 2 * np.pi, 100) r = np.sqrt(1.2) x1 = r * np.cos(theta) x2 = r * np.sin(theta) limit1 = 1.2 colors = ['r', 'g', 'b', 'y', 'lime', 'orange', "black", "darkblue"] for i in range(0, 16, 2): k = int(0.5 * i) ax[0, i].arrow(0, 0, firstrow[k][0], firstrow[k][1], head_width=0.15, head_length=0.15, length_includes_head=True, linewidth=1.5) ax[0, i + 1].arrow(0, 0, firstrow[k][2], firstrow[k][3], head_width=0.15, head_length=0.15, length_includes_head=True, linewidth=1.5) ax[1, i].arrow(0, 0, secondrow[k][0], secondrow[k][1], head_width=0.15, head_length=0.15, length_includes_head=True, linewidth=1.5) ax[1, i + 1].arrow(0, 0, secondrow[k][2], secondrow[k][3], head_width=0.15, head_length=0.15, length_includes_head=True, linewidth=1.5) ax[0, i].plot(x1, x2, color=colors[k]) ax[0, i].set_ylim(-limit1, limit1) ax[0, i].set_xlim(-limit1, limit1) ax[0, i].tick_params(axis='y', labelsize=5) ax[0, i].tick_params(axis='x', labelsize=5) ax[0, i].axis('off') ax[0, i + 1].plot(x1, x2, color=colors[k]) ax[0, i + 1].set_ylim(-limit1, limit1) ax[0, i + 1].set_xlim(-limit1, limit1) ax[0, i + 1].tick_params(axis='y', labelsize=5) ax[0, i + 1].tick_params(axis='x', labelsize=5) ax[0, i + 1].axis('off') ax[1, i].plot(x1, x2, color=colors[k]) ax[1, i].set_ylim(-limit1, limit1) ax[1, i].set_xlim(-limit1, limit1) ax[1, i].tick_params(axis='y', labelsize=5) ax[1, i].tick_params(axis='x', labelsize=5) ax[1, i].axis('off') ax[1, i + 1].plot(x1, x2, color=colors[k]) ax[1, i + 1].set_ylim(-limit1, limit1) ax[1, i + 1].set_xlim(-limit1, limit1) ax[1, i + 1].tick_params(axis='y', labelsize=5) ax[1, i + 1].tick_params(axis='x', labelsize=5) ax[1, i + 1].axis('off') return