```python
import numpy as np
import pykat
import matplotlib.pyplot as plt
```

# Analytics

## Input - output (2x2 matrices)

**Quadrature picture, aka 2-photon formalism** (The one we use for the analytics)

In the frequency domain, the vacuum noise incident on the cavity can in the quadrature-picture be expressed as:

\begin{align}
\mathbf{a}(\Omega) = \begin{bmatrix}
    \hat{a}_1(\Omega) \\
    \hat{a}_2(\Omega)
\end{bmatrix}
\end{align}

The field reflected off an overcoupled cavity ($r_2 = 1$)  is given by:


\begin{align}
\mathbf{b}(\Omega) = \begin{bmatrix}
    \hat{b}_1(\Omega) \\
    \hat{b}_2(\Omega)
\end{bmatrix} = \mathbf{M}\, \mathbf{a}(\Omega)
\end{align}

where the matric describing the reflection off the cavity is:

\begin{align}
\mathbf{M} &= -\mathbf{R} + \mathbf{T} ( \mathbf{I} - \mathbf{N}\, \mathbf{S}\, \mathbf{N}\, \mathbf{R} )^{-1} \mathbf{N}\, \mathbf{S}\, \mathbf{N}\, \mathbf{T} \\ \\
\mathbf{R} &- \text{Reflection off input mirror} \\
\mathbf{T} &- \text{Transmission through input mirror} \\
\mathbf{S} &- \text{Roundtrip propagation} \\
\mathbf{I} &- \text{Identity matrix} \\\\
\end{align}

and the matrix for the non-linear crystal is:

\begin{align}
\mathbf{N} &\equiv \mathbf{N}(r, \phi) = \begin{bmatrix} \cosh(r) + \sinh(r) \cos(2\phi) & \sinh(r) \sin(2\phi) \\ 
                              \sinh(r) \sin(2\phi) & \cosh(r) - \sinh(r) \cos(2\phi)
                     \end{bmatrix}
\end{align}

The homodyne readout power-fluctuations operator is given by:

\begin{align}
\hat{P}(\Omega) =  A_1 \big(\hat{b}_1(\Omega) + \hat{b}_1^\dagger(\Omega) \big) + A_2 \big(\hat{b}_2(\Omega) + \hat{b}_2^\dagger(\Omega) \big)
\end{align}


where $A_1 = \sqrt{P_0}\cos(\theta)$ and $A_2 = \sqrt{P_0}\sin(\theta)$, $P_0$ is the local oscillator power, and $\theta$ is the local oscillator phase. The power spectral density (PSD) and amplitude spectral density (ASD) of the quantum noise becomes:

\begin{align}
\text{PSD} &= ⟨0|\, \hat{P}(\Omega) \hat{P}^\dagger(\Omega) \,| 0 ⟩ =  2 P_0 \hbar \omega_0 \begin{bmatrix} \cos \theta & \sin\theta \end{bmatrix} \mathbf{M} \,\mathbf{M}^\dagger \begin{bmatrix} \cos \theta \\ \sin\theta \end{bmatrix} \\
\text{ASD} &= \sqrt{\text{PSD}}
\end{align}


**Alternatively** Express the power fluctuation operator as:

\begin{align}
\hat{P}(\Omega) =  \mathbf{A}^\text{T} \left(\mathbf{b} + \left(\mathbf{b}^\dagger \right)^\text{T} \right)
\end{align}

where $\mathbf{A}^\text{T} = \sqrt{P_0}\begin{bmatrix} \cos \theta & \sin \theta \end{bmatrix}$, and then 

\begin{align}
\text{PSD} &= ⟨0|\, \hat{P}(\Omega) \hat{P}^\dagger(\Omega) \,| 0 ⟩ =  2 \hbar \omega_0 \mathbf{A}^\text{T} \mathbf{M} \,\mathbf{M}^\dagger \mathbf{A} \\
\text{ASD} &= \sqrt{\text{PSD}}
\end{align}



**Sideband picture, aka 1-photon formalism** (Used internally by Finesse)

\begin{align}
\mathbf{N} &\equiv \mathbf{N}(r, \phi) = \begin{bmatrix} \cosh(r) & \text{e}^{2i\phi} \sinh(r) \\ 
                              \text{e}^{-2i\phi} \sinh(r) & \cosh(r)
                     \end{bmatrix}
\end{align}

\begin{align}
\begin{bmatrix}
\hat{b}_{\omega_0 + \Omega} \\
\hat{b}^\dagger_{\omega_0 - \Omega}  
\end{bmatrix} = \mathbf{N} \begin{bmatrix}
\hat{a}_{\omega_0 + \Omega} \\
\hat{a}^\dagger_{\omega_0 - \Omega}
\end{bmatrix}
\end{align}

# Functions

Functions for doing the analytics:


```python
# Speed of light
c = 299792458.0
# Planck's reduced constant
hbar = 1.0545718e-34

# Propagation through vacuum
def space(phi, L, W=0):
    return np.exp(-1j*W*L/c)*np.array([[np.cos(phi), -np.sin(phi)],
                                       [np.sin(phi), np.cos(phi)]] )

# Transmission
def trans(t):
    return np.eye(2)*t

# Reflection
def refl(r, neg=False):
    if neg:
        r *= -1
    return np.eye(2)*r

# Nonlinear element
def nle(r, phi):
    '''
    Nonlinear element.
    
    r    - Squeeze exponent. 
    phi  - Squeeze angle.
    '''
    return np.array([[np.cosh(r) + np.sinh(r)*np.cos(2*phi), np.sinh(r)*np.sin(2*phi)], 
                     [np.sinh(r)*np.sin(2*phi), np.cosh(r) - np.sinh(r)*np.cos(2*phi)]])


# Reflection off cavity with internal NLE. 
def r_cav(r1, t1, r2, t2, phi, L, r_sq=0, phi_sq=0, W=0):
    # Must to allow second source through end mirror to be correct for quantum noise. 
    # Thus, currenlty incorrect for quantum noise calculations if r2 != 1.
    
    N = nle(r_sq, phi_sq)
    R1p = refl(r1)
    R1n = refl(r1, neg=True)
    T1 = trans(t1)
    R2 = refl(r2)
    S = space(phi, L, W)
    
    # Roundtrip.
    M = np.eye(2) - N@S@R2@S@N@R1p
    M = np.linalg.inv(M)
    
    # M_rt = eye(2) - N*S*R2*S*N*R1p
    # M_rt = M_rt.inv()
        
    return R1n + T1@M@N@S@R2@S@N@T1


def nle_in_cav(r1, t1, r2, t2, phi, L, r_sq, phi_sq, P_LO, F):
    '''
    Computes the ASD of quantum noise reflected off a cavity with an internal NLE.
    
    r1      - Input mirror amplitude reflection
    t1      - Input mirror amplitude transmission
    r2      - End mirror amplitude reflection
    t2      - End mirror amplitude transmission
    phi     - Cavity tuning [deg]
    L       - Cavity length [m]
    r_sq    - Squeeze exponent.
    phi_sq  - Squeeze angle [deg]
    P_LO    - Power of the local oscillator [W].
    F       - Frequency offset at which to compute quantum noise [Hz].
    '''
    
    hbar = 1.0545718e-34
    c0 = 299792458.0

    lam0 = 1.064e-6
    f0 = c0/lam0
    w0 = 2.0*np.pi*f0
    W = 2.0*np.pi*F
    
    phi *= np.pi/180.0
    phi_sq *= np.pi/180.0
    

    # Fields
    N = 5
    a = 0
    b = 2
    c = 4
    d = 6
    e = 8

    SRC = np.zeros([2*N])
    SRC[a:a+2] = 1
    SRC[e:e+2] = 1

    SRC *= np.sqrt(hbar*w0/2.0)

    # The noise matrix
    N_Mat = np.eye(2*N, dtype=complex)*-1

    # a to a (source)
    N_Mat[a:a+2, a:a+2] = np.eye(2)

    # a to b
    N_Mat[b:b+2, a:a+2] = r_cav(r1, t1, r2, t2, phi, L, r_sq, phi_sq, W)

    # b to c
    N_Mat[c:c+2, b:b+2] = refl(np.sqrt(2))

    # b to d
    N_Mat[d:d+2, b:b+2] = trans(np.sqrt(2))

    # e to c
    N_Mat[c:c+2, e:e+2] = trans(np.sqrt(2))

    # e to d
    N_Mat[d:d+2, e:e+2] = refl(np.sqrt(2), neg=True)

    # e to e (source)
    N_Mat[e:e+2, e:e+2] = np.eye(2)

    N = np.linalg.inv(N_Mat)

    # Coherent field
    C = np.sqrt(P_LO)/np.sqrt(2.0)
    D = -np.sqrt(P_LO)/np.sqrt(2.0)

    # ASD readout for both quadratures 
    Q1 = np.sqrt(np.sum(np.abs(N[c,:]*SRC*C - N[d,:]*SRC*D)**2))
    Q2 = np.sqrt(np.sum(np.abs(N[c+1,:]*SRC*C - N[d+1,:]*SRC*D)**2))
    
    return Q1, Q2

# Photodetctor
def pd(A):
    # A is the complex amplitue
    return np.abs(A)**2

# Amplitude detector 
def ad(A):
    return A[0] - 1j*A[1]

# Homodyne readout (not finished. Do not use.)
def HD_readout(M, idx1, idx2, A, a):
    '''
    A    - Coherent part of the local oscillator. 
    a    - Noise part of the local oscillator
    b    - Noise part to be measured
    
    Returns
    --------
    P_HD   - Coherent power
    p_hd   - Noise power operator
    '''
    
    # From cav
    # T1 = trans(1.0/np.sqrt(2.0))
    # R1 = refl(1.0/np.sqrt(2.0))
    
    # From LO
    T = trans(1.0/sqrt(2.0))
    R = refl(1.0/sqrt(2.0), neg=True)
    
    # Coherent
    C = T*A
    D = R*A
    
    # Readout, coherent
    
    P_C = C[0]*C[0].conjugate() + C[1]*C[1].conjugate()
    P_D = D[0]*D[0].conjugate() + D[1]*D[1].conjugate()
    
    P_HD = P_C - P_D
    
    # Readout, quantum noise
    p_c = C[0]*Dagger(c[0]) + C[0].conjugate()*c[0] + C[1]*Dagger(c[1]) + C[1].conjugate()*c[1]
    p_d = D[0]*Dagger(d[0]) + D[0].conjugate()*d[0] + D[1]*Dagger(d[1]) + D[1].conjugate()*d[1]
    
    p_hd = p_c - p_d
    
    return P_HD, p_hd


# Power Spectral Density. Not Finished. Do not use. 
def PSD(p):
    return (p*Dagger(p) + Dagger(p)*p)/2
```

# Cavity response test

No NLE in the cavity. Just testing the "analytic" equations vs. Finesse.


```python
# Controller (setting the parameters here)
# -----------------------------------------

# Input laser power
P0 = 10
# Input laser phase
phi0 = 0

# Power reflection and transmission coeffiecients 
R1 = 0.81
T1 = 1.0-R1

R2 = 0.9
T2 = 1.0-R2

# NLE-parameters
r_sq = 0
phi_sq = 0

# Tuning of end mirror
phi = 0
# Cavity length
L = 1.0
# Frequency offset
F = 10e6
```

## Finesse model


```python
base = pykat.finesse.kat()
base.parse('''

l laser 1 0 0 n0
s s0 0 n0 nIM1

m IM 0.81 0.19 0 nIM1 nIM2
s sCav 1 nIM2 nEM1
m EM 1 0 0 nEM1 nEM2

pd P nIM1
ad A 0 nIM1

yaxis abs:deg
''')     
```


```python
# Finesse run
kat = base.deepcopy()
kat.laser.P = P0
kat.laser.phase = phi0
kat.laser.f = F
kat.A.f = F

kat.IM.R = R1
kat.IM.T = T1
kat.EM.R = R2
kat.EM.T = T2

kat.EM.phi = phi
kat.sCav.L = L

# Sweeping the cavity tuning. Note that positive phi at the end mirror means shorter cavity.
kat.parse('xaxis EM phi lin -90 90 201')

out = kat.run()
P_f = out['P'].real
A_f = out['A']
```

## Analytics


```python
# Analytic
# ---------------

# Transforming parameters
phi0 *= np.pi/180

r1 = np.sqrt(R1)
t1 = np.sqrt(T1)
r2 = np.sqrt(R2)
t2 = np.sqrt(T2)

W = 2.0*np.pi*F

# Input field
E1 = np.sqrt(P0)*np.cos(phi0)
E2 = np.sqrt(P0)*np.sin(phi0)
# ... as vector
E = np.array([[E1],[E2]])

# Arrays for storing output power and complex amplitude
P_a = np.zeros(len(out.x))
A_a = np.zeros(len(out.x), dtype=complex)

for k,phi in enumerate(out.x):
    phi *= np.pi/180
    # Matrix describing reflection off cavity
    M = r_cav(r1, t1, r2, t2, -phi, L, r_sq=0, phi_sq=0, W=W)

    # Output field
    E2 = M @ E
    
    # Computing complex amplitude and power. Rotating 180 deg to make it match Finesse. 
    A_a[k] = ad(E2)*np.exp(-1j*np.pi)
    P_a[k] = pd(A_a[k])
```

## Plotting


```python
fig = plt.figure(dpi=120, figsize=(9,6))
ax = fig.add_subplot(221)
ax.plot(out.x, P_f, label='Finesse')
ax.plot(out.x, P_a,'--', label='Analytic')

ax.legend()
ax.grid()
ax.set_yscale('linear')
ax.set_xlabel(r'Cavity tuning [deg]')
ax.set_ylabel(r'Reflected power [W]')
ax.set_title('Frequency offset: {:.1f} MHz'.format(F/1e6))


# Magnitude
ax = fig.add_subplot(222)
ax.plot(out.x, np.abs(A_f), label='Finesse')
ax.plot(out.x, np.abs(A_a),'--', label='Analytic')

ax.legend()
ax.grid()
ax.set_yscale('linear')
# ax.set_xlabel(r'Cavity tuning [deg]')
ax.set_ylabel(r'Reflected magnitude [$\sqrt{\mathrm{W}}$]')


# Phase
ax = fig.add_subplot(224)
ax.plot(out.x, np.angle(A_f), label='Finesse')
ax.plot(out.x, np.angle(A_a),'--', label='Analytic')

ax.legend()
ax.grid()
ax.set_yscale('linear')
ax.set_xlabel(r'Cavity tuning [deg]')
ax.set_ylabel(r'Phase [deg]')



plt.show(fig)
```

# Crystal in cavity

All parameters for the analytics are extracted from Finesse. Thus, to change parameters, alter the Finesse model below.

## Base Finesse model


```python
base = pykat.finesse.kat()
base.parse('''

l laser 1 0 10 n0
s s0 0 n0 nBS4

dbs FI nFIa nFIb nFIc nFId

s s1 0 nFIc nIM1

m IM 0.9 0.1 0 nIM1 nIM2
s sCav1 0 nIM2 nNLE1
# nle NLE {0} {1} nNLE1 nNLE2
s sCav 1 nNLE2 nEM1
m EM 1 0 0.0 nEM1 nEM2

s s2 0 nFId nBS1

bs BS 0.5 0.5 0 0 nBS1 nBS2 nBS3 nBS4

hd HD 180 nBS2 nBS3
qhd QHD 180 nBS2 nBS3

pd P1 nBS2
pd P2 nBS3

fsig sig 1

yaxis abs
''')
```

DBS

1 -- > 3

2 ---> 1

3 ---> 4

4 ---> 2


```python
# Extracting parameters that will stay the same from the Finesse-model 
r1 = np.sqrt(base.IM.R.value)
t1 = np.sqrt(base.IM.T.value)
r2 = np.sqrt(base.EM.R.value)
t2 = np.sqrt(base.EM.T.value)
L = base.sCav.L.value

# r_sq = np.log(10**(x_sq/10.0))
```

## Sweeping cavity tuning. No squeeze.

Just for checking if numbers agree without squeezing.


```python
# Single-trip squeeze dB
x_sq = 0
# Squeeze angle
phi_sq = 0
# Frequency
F = 1e6

kat = base.deepcopy()

# Adding non-linear element
kat.parse('nle NLE {0} {1} nNLE1 nNLE2'.format(x_sq, phi_sq))

# Setting frequency
kat.signals.f = F

# Adding xaxis and running Finesse.
kat.noxaxis = False
kat.parse('xaxis EM phi lin -90 90 200')
out = kat.run()


# Extracting parameters from kat-instance (to be used in analytic runs)
F = kat.signals.f.value
P_LO = kat.laser.P.value
phi_LO = kat.laser.phase.value*np.pi/180.0
# phi = kat.IM.phi.value - kat.EM.phi.value
phis = -out.x

f0 = c/kat.lambda0
w0 = 2.0*np.pi*f0


r_sq = np.log(10**(x_sq/10.0))



# Running analytics (Big matrix)
# ------------------
# Converting from squeeze magnitude in dB (Finesse) to squeeze exponent r (analytics). 
# Note that we use 10 instead of 20. Thus, using formula for power dB even though we want this in amplitude dB.

Q1 = np.zeros(len(out.x))
Q2 = np.zeros(len(out.x))

for k, phi in enumerate(phis):
    Q1[k], Q2[k] = nle_in_cav(r1, t1, r2, t2, phi, L, r_sq, phi_sq, P_LO, F)

Q = np.real(np.exp(-1j*phi_LO)*(Q1 + 1j*Q2))
    
# Running analytics 2 (2x2 Matrix)
A0 = np.sqrt(P_LO)
alpha = 2*np.sqrt(hbar*w0/2.0)

Q_theta = np.zeros(len(out.x), dtype=float)

for k, phi in enumerate(phis):
    M = r_cav(r1, t1, r2, t2, phi*np.pi/180.0, L, r_sq, phi_sq*np.pi/180.0, W=2.0*np.pi*F)
    
    Q_theta[k] = A0*alpha*np.sqrt((np.abs(M[0,:])**2).sum()*np.cos(phi_LO)**2 +  
                                  (np.abs(M[1,:])**2).sum()*np.sin(phi_LO)**2 + 
                                  2.0*np.cos(phi_LO)*np.sin(phi_LO)*(np.real(M[0,0].conj()*M[1,0]) + 
                                                                      np.real(M[0,1].conj()*M[1,1])) )

```


```python
# Plotting

fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out.x, out['QHD'], label='Finesse')
# Note that the extra factor of sqrt(2) is needed to make a match.
# ax.plot(out.x, Q,'--', label='Analytic')
ax.plot(out.x, Q_theta,':', label='Analytic 2')

ax.legend()
ax.grid()
ax.set_yscale('log')
ax.set_xlabel(r'Cavity tuning [deg]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
ax.set_title('Frequency: {:.1e} Hz. Squeeze: {} dB'.format(F, x_sq))
plt.show(fig)

# Relative difference
print('Relative difference between Finesse and Analytic:   {:.1e}'
      .format( max((Q - out['QHD'])/out['QHD'])))
print('Relative difference between Finesse and Analytic 2: {:.1e}'
      .format( max((Q_theta - out['QHD'])/out['QHD'])))
print('Relative difference between Analytic and Analytic 2: {:.1e}'
      .format( max((Q - Q_theta)/Q_theta)))
```

Good, numbers agree. Then we go ahaed and switch on the NLE!

## Sweeping squeeze magnitude with cavity on resonance


```python
# Single-trip squeeze dB
x_sq = 0
# Squeeze angle
phi_sq = 33
# Frequency
F = 1e0
# Local oscillator phase
phi_LO = 12
# ------------------------------

kat = base.deepcopy()

# Adding non-linear element
kat.parse('nle NLE {0} {1} nNLE1 nNLE2'.format(x_sq, phi_sq))

# Setting frequency
kat.signals.f = F
kat.EM.phi = 12
kat.IM.phi = 0 

# Setting local oscillator phase
kat.laser.phase = phi_LO

# Running once without squeezing
kat.noxaxis = True
out0 = kat.run()

# Adding xaxis and running Finesse.
kat.noxaxis = False
kat.parse('xaxis NLE r lin -1 1 400')


out = kat.run()


# Extracting parameters from kat-instance (to be used in analytic runs)
F = kat.signals.f.value
P_LO = kat.laser.P.value
phi_LO = -kat.laser.phase.value*np.pi/180.0
phi = kat.IM.phi.value - kat.EM.phi.value

f0 = c/kat.lambda0
w0 = 2.0*np.pi*f0

# Squeeze exponent (note 10 in stead of 20...)
r_sqs = np.log(10**(out.x/10.0))



# Running analytics (Big matrix)
# ------------------
# Converting from squeeze magnitude in dB (Finesse) to squeeze exponent r (analytics). 
# Note that we use 10 instead of 20. Thus, using formula for power dB even though we want this in amplitude dB.

Q1 = np.zeros(len(out.x))
Q2 = np.zeros(len(out.x))

for k, r_sq in enumerate(r_sqs):
    Q1[k], Q2[k] = nle_in_cav(r1, t1, r2, t2, phi, L, r_sq, phi_sq, P_LO, F)
    
Q = np.real(np.exp(-1j*phi_LO)*(Q1 + 1j*Q2))
    
# Running analytics 2 (2x2 Matrix)
A0 = np.sqrt(P_LO)
alpha = 2*np.sqrt(hbar*w0/2.0)

A = np.array([[np.cos(phi_LO)],[np.sin(phi_LO)]])

Q_theta = np.zeros(len(out.x), dtype=float)

for k, r_sq in enumerate(r_sqs):
    M = r_cav(r1, t1, r2, t2, phi*np.pi/180.0, L, r_sq, phi_sq*np.pi/180.0, W=2.0*np.pi*F)
    
    Q_theta[k] = A0*alpha*np.sqrt((np.abs(M[0,:])**2).sum()*np.cos(phi_LO)**2 +  
                                  (np.abs(M[1,:])**2).sum()*np.sin(phi_LO)**2 + 
                                  2.0*np.cos(phi_LO)*np.sin(phi_LO)*(np.real(M[0,0].conj()*M[1,0]) + 
                                                                      np.real(M[0,1].conj()*M[1,1])) )
    
    Q_theta[k] = A0*alpha*np.sqrt((A.T@M@M.T.conj()@A)[0,0])

    
r_sq = 0
M = r_cav(r1, t1, r2, t2, phi*np.pi/180.0, L, r_sq, -phi_sq*np.pi/180.0, W=2.0*np.pi*F)

Q0 = A0*alpha*np.sqrt((np.abs(M[0,:])**2).sum()*np.cos(phi_LO)**2 +  
                              (np.abs(M[1,:])**2).sum()*np.sin(phi_LO)**2 + 
                              2.0*np.cos(phi_LO)*np.sin(phi_LO)*(np.real(M[0,0].conj()*M[1,0]) + 
                                                                  np.real(M[0,1].conj()*M[1,1])) )
```


```python
kat.laser.phase
```


```python
# Plotting

fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out.x, out['QHD'], label='Finesse')
# Note that the extra factor of sqrt(2) is needed to make a match.
ax.plot(out.x, Q,'--', label='Analytic')
ax.plot(out.x, Q_theta,':', label='Analytic 2')

ax.legend()
ax.grid()
ax.set_yscale('log')
ax.set_xlabel(r'Single-pass squeeze magnitude [dB]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
ax.set_title('Frequency: {:.1e} Hz. Cavity detuning: {} deg. Squeeze angle: {} deg'.format(F, phi, phi_sq))
plt.show(fig)

# Relative difference
print('Relative difference between Finesse and Analytic:   {:.1e}'
      .format( max((Q - out['QHD'])/out['QHD'])))
print('Relative difference between Finesse and Analytic 2: {:.1e}'
      .format( max((Q_theta - out['QHD'])/out['QHD'])))
print('Relative difference between Analytic and Analytic 2: {:.1e}'
      .format( max((Q - Q_theta)/Q_theta)))
```


```python
Q /out['QHD']
```


```python
# Plotting

fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out.x, 20*np.log10(out['QHD']/out0['QHD']), label='Finesse')
# Note that the extra factor of sqrt(2) is needed to make a match.
# ax.plot(out.x, Q1*np.sqrt(2),'--', label='Analytic')
ax.plot(out.x, 20*np.log10(Q_theta/Q0),':', label='Analytic 2')

ax.legend()
ax.grid()
ax.set_yscale('linear')
ax.set_xlabel(r'Single-pass squeeze magnitude [dB]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
ax.set_title('Frequency: {:.1e} Hz. Cavity detuning: {} deg. Squeeze angle: {} deg'.format(F, phi, phi_sq))
plt.show(fig)

```


```python
kat.IM.R
```


```python

```


```python

```


```python
# Adding xaxis and running Finesse.
kat = base.deepcopy()
kat.noxaxis = False
kat.parse('xaxis NLE r lin -1 1 400')
out = kat.run()

# Extracting some more parameters
F = kat.signals.f.value
P_LO = kat.laser.P.value
phi_LO = kat.laser.phase.value*np.pi/180.0
phi = kat.IM.phi.value - kat.EM.phi.value

f0 = c/kat.lambda0
w0 = 2.0*np.pi*f0
```


```python
# Running analytics (Big matrix)
# ------------------
# Converting from squeeze magnitude in dB (Finesse) to squeeze exponent r (analytics). 
# Note that we use 10 instead of 20. Thus, using formula for power dB even though we want this in amplitude dB.
r_sqs = np.log(10**(out.x/10.0))

Q1 = np.zeros(len(out.x))
Q2 = np.zeros(len(out.x))

for k, r_sq in enumerate(r_sqs):
    Q1[k], Q2[k] = nle_in_cav(r1, t1, r2, t2, phi, L, r_sq, phi_sq, P_LO, F)
```


```python
# Running analytics 2 (2x2 Matrix)
A0 = np.sqrt(P_LO)
alpha = 2*np.sqrt(hbar*w0/2.0)

Q_theta = np.zeros(len(out.x), dtype=float)

for k, r_sq in enumerate(r_sqs):
    M = r_cav(r1, t1, r2, t2, phi*np.pi/180.0, L, r_sq, phi_sq*np.pi/180.0, W=2.0*np.pi*F)
    
    Q_theta[k] = A0*alpha*np.sqrt((np.abs(M[0,:])**2).sum()*np.cos(phi_LO)**2 +  
                                  (np.abs(M[1,:])**2).sum()*np.sin(phi_LO)**2 + 
                                  2.0*np.cos(phi_LO)*np.sin(phi_LO)*(np.real(M[0,0].conj()*M[1,0]) + 
                                                                      np.real(M[0,1].conj()*M[1,1])) )

```


```python
# Plotting

fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out.x, out['QHD'], label='Finesse')
# Note that the extra factor of sqrt(2) is needed to make a match.
ax.plot(out.x, Q1*np.sqrt(2),'--', label='Analytic')
ax.plot(out.x, Q_theta,':', label='Analytic 2')

ax.legend()
ax.grid()
ax.set_yscale('log')
ax.set_xlabel(r'Single-pass squeeze magnitude [dB]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
ax.set_title('Frequency: {:.1e} Hz. Cavity detuning: {} deg'.format(F, phi))
plt.show(fig)
```

# Sweeping cavity tuning

With frequency and squeeze magnitude set to arbitrary values.


```python
# Single-trip squeeze dB
x_sq = 0.01
# Squeeze angle
phi_sq = 39
# Frequency
F = 1e9
# -------------------------

kat = base.deepcopy()

# Adding non-linear element
kat.parse('nle NLE {0} {1} nNLE1 nNLE2'.format(x_sq, phi_sq))

# Setting frequency
kat.signals.f = F

# Adding xaxis and running Finesse.
kat.noxaxis = False
kat.parse('xaxis EM phi lin -90 90 200')
out = kat.run()


# Extracting parameters from kat-instance (to be used in analytic runs)
F = kat.signals.f.value
P_LO = kat.laser.P.value
phi_LO = -kat.laser.phase.value*np.pi/180.0
# phi = kat.IM.phi.value - kat.EM.phi.value
phis = -out.x

f0 = c/kat.lambda0
w0 = 2.0*np.pi*f0


r_sq = np.log(10**(x_sq/10.0))


# Running analytics (Big matrix)
# ------------------
# Converting from squeeze magnitude in dB (Finesse) to squeeze exponent r (analytics). 
# Note that we use 10 instead of 20. Thus, using formula for power dB even though we want this in amplitude dB.

Q1 = np.zeros(len(out.x))
Q2 = np.zeros(len(out.x))

for k, phi in enumerate(phis):
    Q1[k], Q2[k] = nle_in_cav(r1, t1, r2, t2, phi, L, r_sq, phi_sq, P_LO, F)
    
    
# Running analytics 2 (2x2 Matrix)
A0 = np.sqrt(P_LO)
alpha = 2*np.sqrt(hbar*w0/2.0)

A = np.array([[np.cos(phi_LO)],[np.sin(phi_LO)]])

Q_theta = np.zeros(len(out.x), dtype=float)

for k, phi in enumerate(phis):
    M = r_cav(r1, t1, r2, t2, phi*np.pi/180.0, L, r_sq, phi_sq*np.pi/180.0, W=2.0*np.pi*F)
    
    Q_theta[k] = A0*alpha*np.sqrt((np.abs(M[0,:])**2).sum()*np.cos(phi_LO)**2 +  
                                  (np.abs(M[1,:])**2).sum()*np.sin(phi_LO)**2 + 
                                  2.0*np.cos(phi_LO)*np.sin(phi_LO)*(np.real(M[0,0].conj()*M[1,0]) + 
                                                                      np.real(M[0,1].conj()*M[1,1])) )
    
    #Q_theta[k] = A0*alpha*np.sqrt((A.T@M@M.T.conj()@A)[0,0])
```


```python
# Plotting

fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out.x, out['QHD'], label='Finesse')
# Note that the extra factor of sqrt(2) is needed to make a match.
#ax.plot(out.x, Q1*np.sqrt(2),'--', label='Analytic')
ax.plot(out.x, Q_theta,':', label='Analytic 2')

ax.legend()
ax.grid()
ax.set_yscale('log')
ax.set_xlabel(r'Cavity tuning [deg]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
ax.set_title('Frequency: {:.1e} Hz. Squeeze: {} dB. Angle: {} deg'.format(F, x_sq, phi_sq))
plt.show(fig)

# Relative difference
print('Relative difference between Finesse and Analytic:   {:.1e}'
      .format( max((Q1*np.sqrt(2) - out['QHD'])/out['QHD'])))
print('Relative difference between Finesse and Analytic 2: {:.1e}'
      .format( max((Q_theta - out['QHD'])/out['QHD'])))
print('Relative difference between Analytic and Analytic 2: {:.1e}'
      .format( max((Q1*np.sqrt(2) - Q_theta)/Q_theta)))
```


```python

```


```python

```


```python
kat = base.deepcopy()
kat.noxaxis = False
kat.parse('xaxis EM phi lin -10 10 400')
out = kat.run()

# Extracting parameters
F = kat.signals.f.value
P_LO = kat.laser.P.value
phi_LO = kat.laser.phase.value
r_sq = np.log(10**(x_sq/10.0))
```


```python
# Running analytics
phis = -out.x
Q1 = np.zeros(len(out.x))
Q2 = np.zeros(len(out.x))

# ASDs for the two quadratures. Lopoping over the different tunings used in the Finesse run
for k, phi in enumerate(out.x):
    # Note that we set -phi, as positive phi gives shorter cavity in Finesse, while it 
    # gives longer cavity in the analytics. 
    Q1[k], Q2[k] = nle_in_cav(r1, t1, r2, t2, -phi, L, r_sq, phi_sq, P_LO, F)
```


```python
fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out.x, out['QHD'], label='Finesse')
ax.plot(out.x, Q1*np.sqrt(2),'--', label='Analytic')
# ax.plot(out.x, Q_theta*np.sqrt(2),'--', label='Analytic 2')

ax.legend()
ax.grid()
ax.set_yscale('log')
ax.set_xlabel(r'Cavity tuning [deg]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
ax.set_title('Frequency: {:.1e} Hz. SQZ magnitude: {} dB'.format(F, x_sq))
plt.show(fig)
```


```python

```

# Testing SQZ vs. NLE


```python
base_NLE = pykat.finesse.kat()
base_NLE.parse('''

l laser 1 0 0 n0
s s0 0 n0 nBS4

s s1 0 nQN nNLE1

s s2 0 nNLE2 nBS1

bs BS 0.5 0.5 0 0 nBS1 nBS2 nBS3 nBS4

qhd QHD 180 nBS2 nBS3

fsig sig 1

yaxis abs

printnoises
''')

base_SQ = pykat.finesse.kat()
base_SQ.parse('''

l laser 1 0 0 n0
s s0 0 n0 nBS4

sq SQZ 0 0 0 nSQZ
 
s s2 0 nSQZ nBS1

bs BS 0.5 0.5 0 0 nBS1 nBS2 nBS3 nBS4

qhd QHD 180 nBS2 nBS3

fsig sig 1

yaxis abs
''')
```


```python
x_sq = 1
phi_sq = 0

kat_sq = base_SQ.deepcopy()
kat_nle = base_NLE.deepcopy()

# Setting squeeze magnitude and angle
kat_sq.SQZ.db = x_sq
kat_sq.SQZ.angle = phi_sq
kat_nle.parse('nle NLE {0} {1} nNLE1 nNLE2'.format(x_sq, phi_sq))


# Adding x-axis
kat_sq.parse('xaxis SQZ r lin 0 10 100')
kat_nle.parse('xaxis NLE r lin 0 10 100')


out_sq = kat_sq.run()
out_nle = kat_nle.run()
```


```python
fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out_sq.x, out_sq['QHD'], label='SQ')
ax.plot(out_nle.x, out_nle['QHD'],'--', label='NLE')

ax.legend()
ax.grid()
ax.set_yscale('log')
ax.set_xlabel(r'Squeeze magnitude [dB]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
# ax.set_title('Frequency: {:.1e} Hz. SQZ magnitude: {} dB'.format(F, x_sq))
plt.show(fig)
```


```python
fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out_sq.x, out_sq['QHD']/out_sq['QHD'][0], label='SQ')
ax.plot(out_nle.x, out_nle['QHD']/out_sq['QHD'][0],'--', label='NLE')

ax.legend()
ax.grid()
ax.set_yscale('linear')
ax.set_ylim(0,1)
ax.set_xlabel(r'Input squeeze parameter magnitude [dB]')
ax.set_ylabel(r'Relative quantum noise [N(0 dB)]')
# ax.set_title('Frequency: {:.1e} Hz. SQZ magnitude: {} dB'.format(F, x_sq))
plt.show(fig)
```


```python
fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
ax.plot(out_sq.x, 20*np.log10(out_sq['QHD']/out_sq['QHD'][0]), label='SQ')
ax.plot(out_nle.x, 20*np.log10(out_nle['QHD']/out_sq['QHD'][0]),'--', label='NLE')

ax.legend()
ax.grid()
ax.set_yscale('linear')
ax.set_xlabel(r'Squeeze magnitude [dB]')
ax.set_ylabel(r'Quantum noise [W/$\sqrt{\mathrm{Hz}}$]')
# ax.set_title('Frequency: {:.1e} Hz. SQZ magnitude: {} dB'.format(F, x_sq))
plt.show(fig)
```

# Testing with signal

## vs. frequency


```python
base = pykat.finesse.kat()
base.parse('''

l laser 1 0 0 n0
s s0 0 n0 nFIa

dbs FI nFIa nFIb nFIc nFId

s s1 1 nFIc nNLE1

s s2 0 nNLE2 nM1

m M 1 0 0 nM1 nM2

pd1 S1 $fs nFId 

pd P0 nFId 

ad A1 $mfs nFId
ad A2 $fs nFId


fsig sig M phase 1 0

yaxis abs:deg
''')

x_sq = 10
phi_sq = 15

kat = base.deepcopy()
kat.parse('nle NLE {0} {1} nNLE1 nNLE2'.format(x_sq, phi_sq))
kat.parse('xaxis sig f log 1 1G 200')

out = kat.run()
out.plot()
```


```python
fig = plt.figure(dpi=120)

ax = fig.add_subplot(211)
ax.plot(out.x, np.abs(out['A1']), label='Lower sideband')
ax.plot(out.x, np.abs(out['A2']), '--', label='Upper sideband')

ax.grid()
ax.set_yscale('linear')
ax.set_xscale('log')

ax.set_ylabel(r'Sideband amplitude')
ax.legend()

ax = fig.add_subplot(212)
ax.plot(out.x, np.angle(out['A1'])*180/np.pi, label='Lower sideband')
ax.plot(out.x, np.angle(out['A2'])*180/np.pi, '--', label='Upper sideband')

ax.grid()
ax.set_ylabel(r'Sideband phase')
ax.legend()
ax.set_xscale('log')


ax.set_xlabel(r'Frequency [Hz]')
# ax.set_title('Frequency: {:.1e} Hz. SQZ magnitude: {} dB'.format(F, x_sq))
plt.show(fig)
```

## vs squeeze magnitude


```python
base = pykat.finesse.kat()
base.parse('''

l laser 1 0 0 n0
s s0 0 n0 nFIa

dbs FI nFIa nFIb nFIc nFId

s s1 100 nFIc nNLE1

s s2 0 nNLE2 nM1

m M 1 0 0 nM1 nM2

pd1 S1 $fs nFId 

pd P0 nFId 

ad A1 $mfs nFId
ad A2 $fs nFId


fsig sig M phase 100 0

yaxis abs:deg
''')

x_sq = 10
phi_sq = 0

kat = base.deepcopy()
kat.parse('nle NLE {0} {1} nNLE1 nNLE2'.format(x_sq, phi_sq))
kat.parse('xaxis NLE r lin 0 10 200')

out = kat.run()
out.plot()
```


```python
fig = plt.figure(dpi=120)

ax = fig.add_subplot(211)
ax.plot(out.x, np.abs(out['A1']), label='Lower sideband')
ax.plot(out.x, np.abs(out['A2']), '--', label='Upper sideband')

ax.grid()
ax.set_yscale('linear')
ax.set_ylabel(r'Sideband amplitude')
ax.legend()

ax = fig.add_subplot(212)
ax.plot(out.x, np.angle(out['A1'])*180/np.pi, label='Lower sideband')
ax.plot(out.x, np.angle(out['A2'])*180/np.pi, '--', label='Upper sideband')

ax.grid()
ax.set_ylabel(r'Sideband phase')
ax.legend()


ax.set_xlabel(r'Squeeze magnitude [dB]')
# ax.set_title('Frequency: {:.1e} Hz. SQZ magnitude: {} dB'.format(F, x_sq))
plt.show(fig)
```


```python
fig = plt.figure(dpi=120)

ax = fig.add_subplot(211)
ax.plot(out.x, np.abs(out['A1'])/np.abs(out['A1'][0]), label='Lower sideband')
ax.plot(out.x, np.abs(out['A2'])/np.abs(out['A2'][0]), '--', label='Upper sideband')

ax.grid()
ax.set_yscale('linear')
ax.set_ylabel(r'Relative sideband amplitude')
ax.legend()

ax = fig.add_subplot(212)
ax.plot(out.x, 20*np.log10(np.abs(out['A1'])/np.abs(out['A1'][0])), label='Lower sideband')
ax.plot(out.x, 20*np.log10(np.abs(out['A2'])/np.abs(out['A2'][0])), '--', label='Upper sideband')


ax.grid()
ax.set_ylabel(r'Relative sideband amplitude [dB]')
ax.legend()


ax.set_xlabel(r'Squeeze magnitude [dB]')
# ax.set_title('Frequency: {:.1e} Hz. SQZ magnitude: {} dB'.format(F, x_sq))
plt.show(fig)
```


```python

```


```python

```


```python
# Without NLE

kat = pykat.finesse.kat()
kat.parse('''

l laser 1 0 20 n0
s s0 0 n0 nFIa

dbs FI nFIa nFIb nFIc nFId

s s1 0 nFIc nM1

m M 1 0 0 nM1 nM2

pd1 S1 $fs nFId 

pd P0 nFId 

ad A1 $fs nFId


fsig sig M phase 1 0

yaxis abs:deg
xaxis sig f log 1 1G 200
''')
out = kat.run()
out.plot()
```


```python

```


```python
out['S1']
```

264fc12577e49874f5651a8dc2e8eb6feeb4a024
ipynb
Jupyter Notebook
source/daniel_quantum_optics_numerical_and_analytics.ipynb
daccordeon/aliveKat
61230d8fc10846e43a93e4517425a7b2f50366c9
BSD-3-Clause
2020-07-06T07:26:45.000Z
2020-07-06T07:26:45.000Z
source/daniel_quantum_optics_numerical_and_analytics.ipynb
daccordeon/aliveKat
61230d8fc10846e43a93e4517425a7b2f50366c9
BSD-3-Clause
source/daniel_quantum_optics_numerical_and_analytics.ipynb
daccordeon/aliveKat
61230d8fc10846e43a93e4517425a7b2f50366c9
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES