cavity reflection with class

import numpy as np
import matplotlib.pyplot as plt
from cavity_class import Device

c = 299792458

# Emitter parameters
gamma = 2*np.pi*30e6
fe = 484.115e12 # Hz

# Cavity parameters
Q = 12000
kin_factor = 0.8
Ace = 50e9
fc = fe-Ace

# QED
C = 100


device = Device(fc, Q, kin_factor, fe, gamma, C)
f_high, R_high, f_low, R_low = device.find_extrema()
Ass = device.find_overlapping_splitting()
device2 = Device(fc, Q, kin_factor, fe+Ass*1e9, gamma, C)
f_high2, R_high2, f_low2, R_low2 = device2.find_extrema()

fp = np.linspace(fc-device.span, fc+device.span, 1000000)
x_axis = (fp-fc)/1e9

# Plotting reflection amplitude
plt.figure(1)
plt.plot(x_axis, device.R(fp))
plt.plot(x_axis, np.zeros(np.shape(fp)), 'k--')
plt.plot(x_axis, np.ones(np.shape(fp)), 'k--')
plt.scatter(f_high, R_high, c='k')
plt.scatter(f_low, R_low, c='k')
plt.plot(x_axis, device2.R(fp), '--')
plt.scatter(f_high2, R_high2, c='k')
plt.scatter(f_low2, R_low2, c='k')
plt.xlabel('Frequency offset wrt cavity resonance (GHz)')
plt.ylabel('Reflectivity')

# Plotting reflection phase
plt.figure(2)
plt.plot((fp-fc)/1e9, device.phase(fp)*180/np.pi)
plt.plot((fp-fc)/1e9, 180*np.ones(np.shape(fp)), 'k--')
plt.plot((fp-fc)/1e9, -180*np.ones(np.shape(fp)), 'k--')
plt.xlabel('Frequency offset wrt cavity resonance (GHz)')
plt.ylabel('Angle (deg)')
plt.show()