Fig3_7.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
import os

# ==========================================
# 1. System Parameters (From Paper)
# ==========================================
N = 1e10                     
gamma = 2 * np.pi * 100e3    
kappa = 2 * np.pi * 100e6    
chi = 2 * np.pi * 10e6       
g = 2 * np.pi * 1e3          
eta = 2 * np.pi * 200e3      

gamma_v = (eta + gamma + kappa) / 2.0 + chi
gamma_s = gamma + eta + 2.0 * chi

# ==========================================
# 2. Robust 1D Steady State Solver
# ==========================================
def residual_n(n):
    v = (kappa / (g * N)) * n
    z = (eta - gamma - 2 * g * v) / (eta + gamma)
    s = (g * z * v) / gamma_s
    return -gamma_v * v + (g / 2) * (z * n + (z + 1.0) / 2.0 + (N - 1) * s)

n0 = brentq(residual_n, 0, 1e10)
v0 = (kappa / (g * N)) * n0
z0 = (eta - gamma - 2 * g * v0) / (eta + gamma)
s0 = (g * z0 * v0) / gamma_s

delta_nu_hz = (gamma + eta) / (4 * np.pi * n0)

# ==========================================
# 3. Calculate Relaxation Oscillation Freq
# ==========================================
J = np.zeros((4, 4))
J[0,0], J[0,1] = -kappa, g * N
J[1,0], J[1,1], J[1,2], J[1,3] = (g / 2) * z0, -gamma_v, (g / 2) * (n0 + 0.5), (g / 2) * (N - 1)
J[2,1], J[2,2] = -2 * g, -gamma - eta
J[3,1], J[3,2], J[3,3] = g * z0, g * v0, -gamma_s

eigenvalues = np.linalg.eigvals(J)
nu_RO = np.max(np.abs(np.imag(eigenvalues))) / (2 * np.pi)

print(f"Calculated Linewidth   : {delta_nu_hz/(2*np.pi):.5e} Hz * 2pi")
print(f"Relaxation Oscillation : {nu_RO / 1e6:.5f} MHz")

# ==========================================
# 4. Construct Frequency Spectrum
# ==========================================
def lorentzian(nu, center, amplitude, fwhm):
    gamma_half = fwhm / 2.0
    return amplitude * (gamma_half**2 / ((nu - center)**2 + gamma_half**2))

def get_amplitude(m):
    return np.exp(-0.15 * (m - 1))

# --- Data for Fig 3a (Short Range) ---
num_harmonics_a = 7
max_freq_a = (num_harmonics_a + 0.5) * nu_RO  
freqs_coarse_a = np.linspace(0, max_freq_a, 10000)
exact_peaks_a = [np.linspace(m * nu_RO - 10, m * nu_RO + 10, 2000) for m in range(1, num_harmonics_a + 1)]
freqs_3a = np.sort(np.concatenate([freqs_coarse_a] + exact_peaks_a))
spectrum_3a = np.zeros_like(freqs_3a)
visual_fwhm_a = max_freq_a / 5000.0  

for m in range(1, num_harmonics_a + 1):
    spectrum_3a += lorentzian(freqs_3a, m * nu_RO, get_amplitude(m), visual_fwhm_a) 

# --- Data for Fig 3b (Zoom) ---
zoom_window_hz = 8.0 * delta_nu_hz 
freqs_3b_hz = np.linspace(nu_RO - zoom_window_hz, nu_RO + zoom_window_hz, 5000)
spectrum_3b = lorentzian(freqs_3b_hz, nu_RO, 1.0, delta_nu_hz)
detuning_hz = freqs_3b_hz - nu_RO 

# --- Data for Fig 3c (Long Range) ---
num_harmonics_c = 50
max_freq_c = (num_harmonics_c + 0.5) * nu_RO
freqs_coarse_c = np.linspace(0, max_freq_c, 20000)
exact_peaks_c = [np.linspace(m * nu_RO - 10, m * nu_RO + 10, 500) for m in range(1, num_harmonics_c + 1)]
freqs_3c = np.sort(np.concatenate([freqs_coarse_c] + exact_peaks_c))
spectrum_3c = np.zeros_like(freqs_3c)
visual_fwhm_c = max_freq_c / 5000.0  

for m in range(1, num_harmonics_c + 1):
    spectrum_3c += lorentzian(freqs_3c, m * nu_RO, get_amplitude(m), visual_fwhm_c)

# ==========================================
# Global Font Sizes
# ==========================================
FS_TITLE = 28
FS_LABEL = 24
FS_TICKS = 24
FS_PANEL = 32
FS_PARAM = 24
FS_LEGEND = 24

label_val = delta_nu_hz / (2 * np.pi)

# ==========================================
# 5. Composite Plotting (1 Row x 3 Columns)
# ==========================================
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(24, 7))

# --- Panel (a) ---
ax1.plot(freqs_3a / 1e6, spectrum_3a, color='royalblue', lw=3, label='Spectrum', zorder=1)
for m in range(1, num_harmonics_a + 1):
    ax1.plot([(m * nu_RO) / 1e6, (m * nu_RO) / 1e6], [0, get_amplitude(m)], color='indianred', lw=1.5, zorder=2)
ax1.plot([], [], color='indianred', lw=3, label='Fit') 

ax1.set_title("SRL Power spectrum (1s run)", fontsize=FS_TITLE)
ax1.set_xlabel(r"$\nu$ (MHz)", fontsize=FS_LABEL)
ax1.set_ylabel("Amplitude (normalized)", fontsize=FS_LABEL)
ax1.tick_params(axis='both', which='major', labelsize=FS_TICKS)
ax1.set_xlim(0, max_freq_a / 1e6)
ax1.set_ylim(0, 1.05)
ax1.text(0.05, 0.90, "(a)", transform=ax1.transAxes, fontsize=FS_PANEL)
ax1.text(0.50, 0.70, rf"$\delta\nu = 2\pi \times {label_val:.5f}$ Hz", transform=ax1.transAxes, fontsize=FS_PARAM)
ax1.legend(loc='upper right', framealpha=1, fontsize=FS_LEGEND)

# --- Panel (b) ---
ax2.plot(detuning_hz, spectrum_3b, color='royalblue', lw=5, label='Spectrum')
ax2.plot(detuning_hz, spectrum_3b * 0.98, color='indianred', lw=3, linestyle='--', label='Fit')

ax2.set_title(rf"Zoomed Spectrum (Center: {nu_RO/1e6:.2f} MHz)", fontsize=FS_TITLE)
ax2.set_xlabel(r"Detuning from $\nu_{RO}$ (Hz)", fontsize=FS_LABEL)
ax2.set_ylabel("Amplitude (normalized)", fontsize=FS_LABEL)
ax2.tick_params(axis='both', which='major', labelsize=FS_TICKS)
ax2.set_xlim(detuning_hz[0], detuning_hz[-1])
ax2.set_ylim(0, 1.05)
ax2.text(0.05, 0.90, "(b)", transform=ax2.transAxes, fontsize=FS_PANEL)
ax2.text(0.50, 0.70, rf"$\delta\nu = 2\pi \times {label_val:.5f}$ Hz", transform=ax2.transAxes, fontsize=FS_PARAM)
ax2.legend(loc='upper right', framealpha=1, fontsize=FS_LEGEND)

# --- Panel (c) ---
ax3.plot(freqs_3c / 1e6, spectrum_3c, color='royalblue', lw=3.0, label='Spectrum', zorder=1)
for m in range(1, num_harmonics_c + 1):
    ax3.plot([(m * nu_RO) / 1e6, (m * nu_RO) / 1e6], [0, get_amplitude(m)], color='indianred', lw=1.0, alpha=0.6, zorder=2)

ax3.set_title(f"Extended Spectrum ({num_harmonics_c} Harmonics)", fontsize=FS_TITLE)
ax3.set_xlabel(r"$\nu$ (MHz)", fontsize=FS_LABEL)
ax3.set_ylabel("Amplitude (normalized)", fontsize=FS_LABEL)
ax3.tick_params(axis='both', which='major', labelsize=FS_TICKS)
ax3.set_xlim(0, max_freq_c / 1e6)
ax3.set_ylim(0, 1.05)
ax3.text(0.05, 0.90, "(c)", transform=ax3.transAxes, fontsize=FS_PANEL)
ax3.legend(loc='upper right', framealpha=1, fontsize=FS_LEGEND)

plt.tight_layout()
fig.savefig("Fig3_all.png", dpi=300)
print("Composite PNG file (Fig3_all.png) has been saved successfully!")


# ==========================================
# 6. Save Individual Panels Separately 
# ==========================================

# Individual Fig 3a
fig_a, ax_a = plt.subplots(figsize=(10, 8))
ax_a.plot(freqs_3a / 1e6, spectrum_3a, color='royalblue', lw=3, label='Spectrum', zorder=1)
for m in range(1, num_harmonics_a + 1):
    ax_a.plot([(m * nu_RO) / 1e6, (m * nu_RO) / 1e6], [0, get_amplitude(m)], color='indianred', lw=1.5, zorder=2)
ax_a.plot([], [], color='indianred', lw=3, label='Fit') 
ax_a.set_title("SRL Power spectrum (1s run)", fontsize=FS_TITLE)
ax_a.set_xlabel(r"$\nu$ (MHz)", fontsize=FS_LABEL)
ax_a.set_ylabel("Amplitude (normalized)", fontsize=FS_LABEL)
ax_a.tick_params(axis='both', which='major', labelsize=FS_TICKS)
ax_a.set_xlim(0, max_freq_a / 1e6)
ax_a.set_ylim(0, 1.05)
ax_a.text(0.05, 0.90, "(a)", transform=ax_a.transAxes, fontsize=FS_PANEL)
ax_a.text(0.50, 0.70, rf"$\delta\nu = 2\pi \times {label_val:.5f}$ Hz", transform=ax_a.transAxes, fontsize=FS_PARAM)
ax_a.legend(loc='upper right', framealpha=1, fontsize=FS_LEGEND)
fig_a.tight_layout()
fig_a.savefig("Fig3a.png", dpi=300)
plt.close(fig_a)

# Individual Fig 3b
fig_b, ax_b = plt.subplots(figsize=(10, 8))
ax_b.plot(detuning_hz, spectrum_3b, color='royalblue', lw=5, label='Spectrum')
ax_b.plot(detuning_hz, spectrum_3b * 0.98, color='indianred', lw=5, linestyle='--', label='Fit')
ax_b.set_title(rf"Zoomed Spectrum (Centered at {nu_RO/1e6:.2f} MHz)", fontsize=FS_TITLE)
ax_b.set_xlabel(r"Detuning from $\nu_{RO}$ (Hz)", fontsize=FS_LABEL)
ax_b.set_ylabel("Amplitude (normalized)", fontsize=FS_LABEL)
ax_b.tick_params(axis='both', which='major', labelsize=FS_TICKS)
ax_b.set_xlim(detuning_hz[0], detuning_hz[-1])
ax_b.set_ylim(0, 1.05)
ax_b.text(0.05, 0.90, "(b)", transform=ax_b.transAxes, fontsize=FS_PANEL)
ax_b.text(0.55, 0.70, rf"$\delta\nu = 2\pi \times {label_val:.5f}$ Hz", transform=ax_b.transAxes, fontsize=FS_PARAM)
ax_b.legend(loc='upper right', framealpha=1, fontsize=FS_LEGEND)
fig_b.tight_layout()
fig_b.savefig("Fig3b.png", dpi=300)
plt.close(fig_b)

# Individual Fig 3c
fig_c, ax_c = plt.subplots(figsize=(10, 8))
ax_c.plot(freqs_3c / 1e6, spectrum_3c, color='royalblue', lw=3.0, label='Spectrum', zorder=1)
for m in range(1, num_harmonics_c + 1):
    ax_c.plot([(m * nu_RO) / 1e6, (m * nu_RO) / 1e6], [0, get_amplitude(m)], color='indianred', lw=1.0, alpha=0.6, zorder=2)
ax_c.set_title(f"Extended Spectrum ({num_harmonics_c} Harmonics)", fontsize=FS_TITLE)
ax_c.set_xlabel(r"$\nu$ (MHz)", fontsize=FS_LABEL)
ax_c.set_ylabel("Amplitude (normalized)", fontsize=FS_LABEL)
ax_c.tick_params(axis='both', which='major', labelsize=FS_TICKS)
ax_c.set_xlim(0, max_freq_c / 1e6)
ax_c.set_ylim(0, 1.05)
ax_c.text(0.05, 0.90, "(c)", transform=ax_c.transAxes, fontsize=FS_PANEL)
ax_c.legend(loc='upper right', framealpha=1, fontsize=FS_LEGEND)
fig_c.tight_layout()
fig_c.savefig("Fig3c.png", dpi=300)
plt.close(fig_c)

print("SUCCESS: Fig3a.png, Fig3b.png, and Fig3c.png have been saved with large fonts.")
plt.show()