diff --git a/FrequencyConversion/FC_solve_EoM.py b/FrequencyConversion/FC_solve_EoM.py new file mode 100644 index 0000000..ed63112 --- /dev/null +++ b/FrequencyConversion/FC_solve_EoM.py @@ -0,0 +1,237 @@ +"""Functions to calculate the numerical solution to the FC equations of motion using linear algebra""" + +import numpy as np +from scipy.linalg import expm +from scipy.special import hermite + +# pylint: disable=invalid-name +# pylint: disable=too-many-arguments +# pylint: disable=too-many-locals +# pylint: disable=consider-using-enumerate + + +def dag(M): + """ + Gives the Hermitian conjugate of input matrix. + + Args: + M (array): input matrix + + Returns: + array: the Hermitian conjugate matrix + """ + return M.conj().T + + +def pulse_shape(dw, sig, n=0): + """ + Gives a Hermite-Gaussian pulse shape. + + Args: + dw (array): frequency vector or grid + sig (float): width of the pulse + n (int): index of Hermite polynomial used. Default is 0 (first Hermite polynomial) + + Returns: + array: pulse shape + """ + gaussian = np.exp(-(dw**2) / (2 * sig**2)) / (np.pi * sig**2)**(1 / 4) + return gaussian * hermite(n)(dw / sig) + + +def get_blocks(M): + """ + Separates a matrix into 4 blocks. + + Args: + M (array): matrix to split + + Returns: + (array, :top left block + array, :top right block + array, :bottom left block + array) :bottom right block + """ + side = M.shape[0] + Ua = M[0 : side // 2, 0 : side // 2] + Va = M[0 : side // 2, side // 2 : side] + Vc = M[side // 2 : side, 0 : side // 2] + Uc = M[side // 2 : side, side // 2 : side] + return Ua, Va, Vc, Uc + + +def get_prop(dw, vp, va, vc, pump, L, D=1): + """ + Generates the propagator given the crystal and pump parameters, and a frequency vector. + + Args: + dw (array): frequency vector + vp (float): group velocity for pump + va (float): group velocity for a beam + vc (float): group velocity for c beam + pump (array): pump pulse shape as a matrix (already scaled) + L (float): length of propagation in the crystal + D (float): number depicting the strength of the nonlinear interaction + + Returns: + (array, :logm of the propagator + array) :propagator + """ + + # Upper left and bottom right (diagonal) blocks + ka = 1 / va - 1 / vp + kc = 1 / vc - 1 / vp + + Ka = np.diag(1j * dw * ka) + Kc = np.diag(1j * dw * kc) + + # Bottom left and upper right (non-diagonal) blocks + offdiag = np.conj(D) * pump * (dw[-1] - dw[0]) / (dw.size - 1) + + J = np.block([[Ka, -1j * offdiag], [-1j * offdiag.conj().T, Kc]]) + prop = expm(J * L) + return J, prop + + +def check_prop(prop, print_checks=False): + """ + Verifies that the obtained propagator fulfills the 6 conditions imposed in the paper + by Andreas Christ from 2013 (equations A3-4 and A6-7). + Returns the maximum error for each condition. + + Args : + prop (matrix): the propagator to perform the checks on + print_checks (bool): if True, prints if the condition is True or False. False by default. + + Returns : + array: maximum value of the absolute value of the condition mismatch, + for each condition (always an array of 6 floats) + """ + Ua, Va, Vc, Uc = get_blocks(prop) + Vc = -Vc + N = Ua.shape[0] + cond1a = Ua @ dag(Ua) + Va @ dag(Va) + cond1b = Uc @ dag(Uc) + Vc @ dag(Vc) + cond3a = dag(Ua) @ Ua + dag(Vc) @ Vc + cond3b = dag(Uc) @ Uc + dag(Va) @ Va + cond2 = Ua @ dag(Vc) - Va @ dag(Uc) + cond4 = dag(Ua) @ Va - dag(Vc) @ Uc + I = np.identity(N) + O = np.zeros((N, N)) + conds = np.array([cond1a, cond1b, cond2, cond3a, cond3b, cond4]) + comps = np.array([I, I, O, I, I, O]) + labels = ["1a", "1b", "2", "3a", "3b", "4"] + to_return = np.array([]) + for i in range(len(labels)): + cond = conds[i] + comp = comps[i] + txt = "Condition {name} is {TF}" + if print_checks: + print(txt.format(name=labels[i], TF=np.allclose(cond, comp))) + to_return = np.append(to_return, np.max(np.abs(cond - comp))) + return to_return + + +def remove_phases(prop, dw, L, vp, va, vc): + """ + Removes the free space propagation phases from the propagator. + + Args: + prop (array): propagator with phases + dw (array): frequency vector + L (float): length of propagation in crystal + vp (float): group velocity for pump + va (float): group velocity for a beam + vc (float): group velocity for c beam + + Returns: + array :the phase-free propagator as a matrix + """ + kc = 1 / vp - 1 / vc + ka = 1 / va - 1 / vp + + # separate into blocks + Ua, Va, Vc, Uc = get_blocks(prop) + # remove phases + Ua = np.diag(np.exp(-1j * dw * ka * L / 2)) @ Ua @ np.diag(np.exp(-1j * dw * ka * L / 2)) + Va = np.diag(np.exp(-1j * dw * ka * L / 2)) @ Va @ np.diag(np.exp(1j * dw * kc * L / 2)) + Vc = np.diag(np.exp(1j * dw * kc * L / 2)) @ Vc @ np.diag(np.exp(-1j * dw * ka * L / 2)) + Uc = np.diag(np.exp(1j * dw * kc * L / 2)) @ Uc @ np.diag(np.exp(1j * dw * kc * L / 2)) + + return np.block([[Ua, Va], [Vc, Uc]]) + + +def get_JCA(prop, dw): + """ + Gives the Joint Spectral Amplitude based on the propagator. + + Args: + prop (array): propagator with phases + dw (array): frequency vector + + Returns: + array :JCA as a matrix + """ + Ua, _, Vc, _ = get_blocks(prop) + M = Ua @ np.conj(Vc).T + U, s, V = np.linalg.svd(M) + S = np.diag(s) + R = np.arcsin(2 * S) / 2 + delta_omega = (dw[-1] - dw[0]) / (dw.size - 1) + J = U @ R @ V / delta_omega + return J + + +def get_schmidt(prop): + """ + Gives average number of photons, eigenfunctions (Schmidt modes) + and eigenvalues (Schmidt coefficients squared) from propagator. + This is done using an eigendecomposition. + For the method with straight-up svds, see get_schmidt2. + + Args: + prop(array): propagator + dw (array): frequency vector + L (float): length of propagation crystal + vp (float): group velocity for pump + va (float): group velocity for a beam + vc (float): group velocity for c beam + + Returns: + (float, :avg number of photons + array, :Schmidt modes + array) :Schmidt coefficients squared + """ + _, Va, _, _ = get_blocks(prop) + Numa = Va @ np.conj(Va).T # Numbers matrix for a-beam + Na = np.real(np.trace(Numa)) # average number of photons in a-beam + val, u = np.linalg.eigh(Numa) # schmidt coefficients squared and modes from eigendecomposition + val = np.flip(val) + u = np.flip(u, axis=1) + return Na, val, u + + +def get_schmidt2(prop): + """ + Alternative method to get Schmidt modes and coefficients + by doing svds of the propagator blocks directly. + + Args : + prop (array): the propagator + + Returns : + (array, :the U matrix from the svd of Va + array, :the array of singular values from the svd of Va + array, :the V matrix from the svd of Va + array, :the U matrix from the svd of Vc + array, :the array of singular values from the svd of Vc + array) :the V matrix from the svd of Vc + """ + _, Va, Vc, _ = get_blocks(prop) + varphiV, sinrk1, phiV = np.linalg.svd(Va) + xiV, sinrk2, psiV = np.linalg.svd(Vc) + varphiV = np.flip(varphiV, axis=0) + phiV = np.flip(phiV, axis=0) + xiV = np.flip(xiV, axis=0) + psiV = np.flip(psiV, axis=0) + return np.conj(varphiV), sinrk1, phiV.T, np.conj(xiV), sinrk2, psiV.T diff --git a/FrequencyConversion/example_solve_EoM.ipynb b/FrequencyConversion/example_solve_EoM.ipynb new file mode 100644 index 0000000..7cab520 --- /dev/null +++ b/FrequencyConversion/example_solve_EoM.ipynb @@ -0,0 +1,9303 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Short example of usage for generating and analyzing the propagator for frequency conversion, for an unpoled waveguide." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from FC_solve_EoM import get_prop, get_JCA, get_schmidt, get_schmidt2, pulse_shape, remove_phases\n", + "%config InlineBackend.figure_formats=['svg']\n", + "%load_ext autoreload\n", + "%autoreload 2\n", + "np.set_printoptions(linewidth=170)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We begin by defining the parameters associated with frequency conversion (of the pump and crystal)." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "vc, vp, va = 1 / 1.5, 1 / 3, 1 / 4.5 # group velocities\n", + "\n", + "dw_range = 10 # Range of frequency values\n", + "N = 301 # Number of frequency points\n", + "dw = np.linspace(-dw_range, dw_range, N) # Frequency vector\n", + "L = 1.35 # Length of propagation in FC medium\n", + "sig = 1.1 # Width of the (Hermite-Gaussian) pump shape\n", + "amp = 0.8 # Pump amplitude\n", + "pump = amp * pulse_shape(dw[:, np.newaxis] - dw, sig, n = 0) # Gaussian-shaped pump" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We then generate the propagator by calling the appropriate function. We then remove its free propagation phases, so we change back to the original electric field operators." + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "_, prop = get_prop(dw, vp, va, vc, pump, L) # Matrix to exponentiate and propagator\n", + "prop = remove_phases(prop, dw, L, vp, va, vc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can visualize the input pulses by defining them and comparing them to the pump pulse shape." + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "# Initial pulses\n", + "initial_Ea = np.zeros(dw.size)\n", + "initial_Ec = 0.05 * pulse_shape(dw, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:24.451744\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(dw, initial_Ea, label=\"Initial $E_a$\")\n", + "plt.plot(dw, initial_Ec, label=\"Initial $E_c$\")\n", + "plt.plot(dw, pulse_shape(dw, sig), label=\"Normalized pump pulse\")\n", + "plt.ylabel(\"Amplitude\")\n", + "plt.title(\"Input pulses\")\n", + "plt.legend(loc=\"upper left\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To get what the output pulses look like, we do matrix multiplication with the propagator, with its phases removed." + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "kc = 1 / vp - 1 / vc\n", + "ka = 1 / va - 1 / vp\n", + "sol = prop @ np.concatenate((initial_Ea, initial_Ec))\n", + "Ea, Ec = np.split(sol, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'a beam / down-converted beam')" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:24.720872\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(dw, np.real(Ea), label=\"Re\")\n", + "plt.plot(dw, np.imag(Ea), label=\"Im\")\n", + "plt.plot(dw, np.abs(Ea), label=\"Abs\")\n", + "plt.legend()\n", + "plt.xlabel(\"$\\delta\\omega$\")\n", + "plt.ylabel(\"Amplitude\")\n", + "plt.title(\"a beam / down-converted beam\")" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'c beam / non-converted beam')" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:24.914946\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(dw, np.real(Ec), label=\"Re\")\n", + "plt.plot(dw, np.imag(Ec), label=\"Im\")\n", + "plt.plot(dw, np.abs(Ec), label=\"Abs\")\n", + "plt.legend()\n", + "plt.xlabel(\"$\\delta\\omega$\")\n", + "plt.ylabel(\"Amplitude\")\n", + "plt.title(\"c beam / non-converted beam\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To get the Joint Conversion Amplitude, we call the appropriate function and show the matrix." + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "J = get_JCA(prop, dw)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, '$\\\\delta\\\\omega_{out}$')" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:25.317534\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(np.abs(J), extent=[dw.min(), dw.max(), dw.min(), dw.max()])\n", + "plt.title(\"JCA\")\n", + "plt.xlabel(\"$\\delta\\omega_{in}$\")\n", + "plt.ylabel(\"$\\delta\\omega_{out}$\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To extract information from the propagator, we can do a Schmidt decomposition. This yields the same information as an eigendecomposition of the numbers matrix, which is what we do by calling `get_schmidt`. From the function, we get the average number of photons generated, the vector of eigenvalues (Schmidt coefficients squared) and the eigenvectors (Schmidt modes)." + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Average number of signal photons is 9.751731e-01\n" + ] + } + ], + "source": [ + "Na, val, u = get_schmidt(prop)\n", + "txt = \"Average number of signal photons is {val:e}\"\n", + "print(txt.format(val=Na))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can then select a mode to plot. Here is the first mode." + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'Schmidt mode 1')" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:25.618635\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "mark = 1 # number of the Schmidt mode to plot\n", + "\n", + "mode = np.exp(-1j * np.angle(u[(N - 1) // 2, mark-1])) * u[:, mark-1] # applying a phase so that the imaginary part is 0 at 0\n", + "plt.plot(dw, np.real(mode), label=\"Re\")\n", + "plt.plot(dw, np.imag(mode), label=\"Im\")\n", + "plt.plot(dw, np.abs(mode), label=\"Abs\")\n", + "plt.legend()\n", + "plt.ylabel(\"Amplitude\")\n", + "plt.title(\"Schmidt mode \" + str(mark))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also make a bar plot of the conversion efficiencies for each mode. Here we do the first 10 modes." + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.0, 1.1)" + ] + }, + "execution_count": 35, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:25.818288\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "effs = val[:10] # efficiencies for the first 10 modes\n", + "plt.bar(x=np.arange(1, 11), height=effs)\n", + "plt.xlabel(\"Schmidt mode number\")\n", + "plt.ylabel(\"$sin^2(r_k)$\")\n", + "plt.title(\"FC Efficiency\")\n", + "plt.ylim(0, 1.1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To check that our results are good, we can also make use of the `get_schmidt2` function, which does the same Schmidt decomposition but by doing the SVDs directly." + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "u, val, *_ = get_schmidt2(prop)" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'Schmidt mode 1')" + ] + }, + "execution_count": 37, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:26.098871\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "mark = 1 # number of the Schmidt mode to plot\n", + "\n", + "mode = np.exp(-1j * np.angle(u[(N - 1) // 2, mark-1])) * u[:, mark-1] # applying a phase so that the imaginary part is 0 at 0\n", + "plt.plot(dw, np.real(mode), label=\"Re\")\n", + "plt.plot(dw, np.imag(mode), label=\"Im\")\n", + "plt.plot(dw, np.abs(mode), label=\"Abs\")\n", + "plt.legend()\n", + "plt.ylabel(\"Amplitude\")\n", + "plt.title(\"Schmidt mode \" + str(mark))" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.0, 1.1)" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-10-02T12:33:26.270307\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.7.2, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "effs = val[:10]**2 # efficiencies for the first 10 modes. We square them because SVDs yield the Schmidt coefficients not squared.\n", + "plt.bar(x=np.arange(1, 11), height=effs)\n", + "plt.xlabel(\"Schmidt mode number\")\n", + "plt.ylabel(\"$sin^2(r_k)$\")\n", + "plt.title(\"FC Efficiency\")\n", + "plt.ylim(0, 1.1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As we can see, the SVD results match the eigendecomposition results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "base", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/FrequencyConversion/test_fermionic_bloch_messiah.py b/FrequencyConversion/test_fermionic_bloch_messiah.py new file mode 100644 index 0000000..6664743 --- /dev/null +++ b/FrequencyConversion/test_fermionic_bloch_messiah.py @@ -0,0 +1,47 @@ +"""Tests for the implementation of the fermionic bloch messiah decomposition""" + +import pytest +import numpy as np +from FC_solve_EoM import get_prop, get_JCA, get_schmidt, get_schmidt2, pulse_shape, get_blocks +from fermionic_bloch_messiah import fermionic_bloch_messiah_1 + +# pylint: disable=invalid-name +# pylint: disable=too-many-arguments +# pylint: disable=too-many-locals +# pylint: disable=consider-using-enumerate + +@pytest.mark.parametrize("N", [101, 301, 1001]) +@pytest.mark.parametrize("amp", np.arange(0.01, 0.31, 0.03)) +def test_fbm_decomp(amp, N): + """Tests that the product of decomposition matrices gives back the propagator""" + vc, vp, va = 1 / 1.5, 1 / 3, 1 / 4.5 # group velocities + dw_range = 10 # Range of frequency values + dw = np.linspace(-dw_range, dw_range, N) # frequency difference vector + L = 1.35 # length of propagation in FC medium + sig = 1.1 # width of the (Hermite-Gaussian) pump shape + pump = amp * pulse_shape(dw[:, np.newaxis] - dw, sig) # pump + _, prop = get_prop(dw, vp, va, vc, pump, L) # propagator + + A, D, Bdag = fermionic_bloch_messiah_1(prop) + + assert np.allclose(A @ D @ Bdag, prop, atol=1e-7) + # We set absolute tolerance to 1e-7 to pass the tests, but even with that on test is failed + + +@pytest.mark.parametrize("N", [101, 301, 1001]) +@pytest.mark.parametrize("amp", np.linspace(0.01, 0.3, 10)) +def test_fbm_sv(amp, N): + """Tests that the singular values match the expected cos^2+sin^2 relation""" + vc, vp, va = 1 / 1.5, 1 / 3, 1 / 4.5 # group velocities + dw_range = 10 # Range of frequency values + dw = np.linspace(-dw_range, dw_range, N) # frequency difference vector + L = 1.35 # length of propagation in FC medium + sig = 1.1 # width of the (Hermite-Gaussian) pump shape + pump = amp * pulse_shape(dw[:, np.newaxis] - dw, sig) # pump + _, prop = get_prop(dw, vp, va, vc, pump, L) # propagator + + A, D, Bdag = fermionic_bloch_messiah_1(prop) + + Du, Dv, *_ = get_blocks(D) + + assert np.allclose(Du @ Du + Dv @ Dv, np.identity(N)) diff --git a/FrequencyConversion/test_solve_EoM.py b/FrequencyConversion/test_solve_EoM.py new file mode 100644 index 0000000..32cc57c --- /dev/null +++ b/FrequencyConversion/test_solve_EoM.py @@ -0,0 +1,103 @@ +"""Tests for the solve_EoM functions""" + +import pytest +import numpy as np +from scipy.optimize import curve_fit +from FC_solve_EoM import get_prop, get_JCA, get_schmidt, get_schmidt2, pulse_shape, get_blocks + +# pylint: disable=invalid-name +# pylint: disable=too-many-arguments +# pylint: disable=too-many-locals +# pylint: disable=consider-using-enumerate + + +@pytest.mark.parametrize("N", [101, 301, 1001]) +@pytest.mark.parametrize("amp", np.arange(0.01, 0.3, 0.02)) +def test_prop_blocks_1(N, amp): + """Tests that the propagator blocks follow the right relations""" + vc, vp, va = 1 / 1.5, 1 / 3, 1 / 4.5 # group velocities + dw_range = 10 # Range of frequency values + dw = np.linspace(-dw_range, dw_range, N) # frequency difference vector + L = 1.35 # length of propagation in FC medium + sig = 1.1 # width of the (Hermite-Gaussian) pump shape + pump = amp * pulse_shape(dw[:, np.newaxis] - dw, sig) # pump + _, prop = get_prop(dw, vp, va, vc, pump, L) # matrix to exponentiate and propagator + + Ua, Va, Vc, Uc = get_blocks(prop) + Vc = -Vc + + cond1a = np.allclose(Ua @ Ua.T.conj() + Va @ Va.T.conj(), np.identity(N)) + cond1b = np.allclose(Uc @ Uc.T.conj() + Vc @ Vc.T.conj(), np.identity(N)) + assert cond1a and cond1b + + +@pytest.mark.parametrize("N", [101, 301, 1001]) +@pytest.mark.parametrize("amp", np.linspace(0.01, 0.3, 10)) +def test_prop_blocks_2(N, amp): + """Tests that the propagator blocks follow the right relations""" + vc, vp, va = 1 / 1.5, 1 / 3, 1 / 4.5 # group velocities + dw_range = 10 # Range of frequency values + dw = np.linspace(-dw_range, dw_range, N) # frequency difference vector + L = 1.35 # length of propagation in FC medium + sig = 1.1 # width of the (Hermite-Gaussian) pump shape + pump = amp * pulse_shape(dw[:, np.newaxis] - dw, sig) # pump + _, prop = get_prop(dw, vp, va, vc, pump, L) # matrix to exponentiate and propagator + + Ua, Va, Vc, Uc = get_blocks(prop) + Vc = -Vc + + assert np.allclose(Ua @ Vc.T.conj(), Va @ Uc.T.conj()) + + +@pytest.mark.parametrize("N", [101, 301]) +@pytest.mark.parametrize("amp", np.linspace(0.01, 0.3, 10)) +def test_noint(N, amp): + """Tests that when the nonlinear interaction is off, the propagator is diagonal""" + vc, vp, va = 1 / 1.5, 1 / 3, 1 / 4.5 # group velocities + dw_range = 10 # Range of frequency values + dw = np.linspace(-dw_range, dw_range, N) # frequency difference vector + L = 1.35 # length of propagation in FC medium + sig = 1.1 # width of the (Hermite-Gaussian) pump shape + pump = amp * pulse_shape(dw[:, np.newaxis] - dw, sig) # pump + _, prop = get_prop(dw, vp, va, vc, pump, L, D=0) # propagator + + assert np.allclose(prop - np.diag(np.diagonal(prop)), np.zeros(prop.shape)) + + +@pytest.mark.parametrize("N", [101, 301]) +@pytest.mark.parametrize("L", np.arange(0.5, 3, 0.5)) +def test_JCA_low_gain(N, L): + """Tests that the JCA at low gain matches a sinc times a gaussian""" + vc, vp, va = 1 / 1.5, 1 / 3, 1 / 4.5 # group velocities + dw_range = 10 # Range of frequency values + dw = np.linspace(-dw_range, dw_range, N) # frequency difference vector + sig = 1.1 # width of the (Hermite-Gaussian) pump shape + amp = 1e-2 # pump amplitude for low gain + pump = amp * pulse_shape(dw[:, np.newaxis] - dw, sig) # pump + _, prop = get_prop(dw, vp, va, vc, pump, L, D=0) # propagator + + def phase_matching(dw_a, dw_c, L): + dk = dw_a * (1 / va - 1 / vp) + dw_c * (1 / vp - 1 / vc) + return L / 2 * np.sinc(dk * L / 2) + + def JSA_model(X, Y, sig, L, amp): + return np.abs(phase_matching(X, Y, L) * amp * pulse_shape(X - Y, sig)) + + def _JSA_model(M, sig, L, amp): + x, y = M + return JSA_model(x, y, sig, L, amp) + + J = get_JCA(prop) + + X, Y = np.meshgrid(dw, dw) + + xdata = np.vstack((np.ravel(X), np.ravel(Y))) + ydata = np.ravel(np.abs(J)) + p0 = (sig, amp, 6.28) + + popt, pcov = curve_fit(_JSA_model, xdata, ydata, p0) + + fit = JSA_model(X, Y, popt[0], popt[1], popt[2]) + + assert np.allclose(J, fit) +