#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jun  5 13:30:43 2026

@author: leonhardtm
"""


import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path

NN=np.array([20,30,50,100,200,500,1000,5000,100000])

cmap = plt.colormaps['gist_rainbow'] 
C=np.array([0.1,0.18,0.22,0.3,0.53,0.65,0.7,0.85,0.9])
colors = cmap(C)
 
plt.rcParams['figure.figsize'] = 24,5
plt.rcParams['axes.labelsize']=22
fig,axs=plt.subplots(1,3)


axs[0].tick_params(direction='out',labelsize=22)
axs[1].tick_params(direction='out',labelsize=22)
axs[2].tick_params(direction='out',labelsize=22)

### (a) Superradiant non-oscillatory solutions
for i in range(len(NN)):
    X=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG5"/'SR_non-oscillatory_{}.csv'.format(NN[i]))
    X_data=X.to_numpy()
    axs[0].plot(X_data[:,0],X_data[:,1],color=colors[i],label='{}'.format(NN[i]))
    axs[0].fill_between(X_data[:,0],X_data[:,1]-2*X_data[:,2],X_data[:,1]+2*X_data[:,2],alpha=0.3,color=colors[i],edgecolor='none')
    
MF=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG5"/'SR_non-oscillatory_MF.csv')
MF_data=MF.to_numpy()
axs[0].plot(MF_data[:,0],MF_data[:,1],'--',color='black',label=r'$\infty$ (MF)')

### (b) Superradiant oscillatory solutions
for i in range(len(NN)):
    X=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG5"/'SR_oscillations_{}.csv'.format(NN[i]))
    X_data=X.to_numpy()
    axs[1].plot(X_data[:,0],X_data[:,1],color=colors[i],label='{}'.format(NN[i]))
    axs[1].fill_between(X_data[:,0],X_data[:,1]-2*X_data[:,2],X_data[:,1]+2*X_data[:,2],alpha=0.3,color=colors[i],edgecolor='none')
    
MF=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG5"/'SR_oscillations_MF.csv')
MF_data=MF.to_numpy()
axs[1].plot(MF_data[:,0],MF_data[:,1],'--',color='black',label=r'$\infty$ (MF)')


### (c) Normal Phase oscillatory solutions
for i in range(len(NN)):
    X=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG5"/'NP_oscillations_{}.csv'.format(NN[i]))
    X_data=X.to_numpy()
    axs[2].plot(X_data[:,0],X_data[:,1],color=colors[i],label='{}'.format(NN[i]))
    axs[2].fill_between(X_data[:,0],X_data[:,1]-2*X_data[:,2],X_data[:,1]+2*X_data[:,2],alpha=0.3,color=colors[i],edgecolor='none')
    
MF=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG5"/'NP_oscillations_MF.csv')
MF_data=MF.to_numpy()
axs[2].plot(MF_data[:,0],MF_data[:,1],'--',color='black',label=r'$\infty$ (MF)')


axs[1].legend(ncols=10,title=r'Emitters per Cluster $(N)$',fontsize=22,title_fontsize=22,loc=8)

plt.setp(axs[0], xlabel=r'$t(\omega_C^{-1})$')
plt.setp(axs[1], xlabel=r'$t(\omega_C^{-1})$')
plt.setp(axs[2], xlabel=r'$t(\omega_C^{-1})$')
plt.setp(axs[0], ylabel=r'$\langle \hat{S}^x_1 \rangle/N$')
plt.setp(axs[1], ylabel=r'$\langle \hat{S}^x_1 \rangle/N$')
plt.setp(axs[2], ylabel=r'$\langle \hat{S}^x_1 \rangle/N$')





