#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jun  4 17:34:52 2026

@author: leonhardtm
"""

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



##### (c) Standard Dicke Model Plot
## Parameters
# wA=0.1
# wC=1
# kappa=1
# gc=0.5*np.sqrt(wA*(kappa**2+wC**2)/wC)

NP=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG1"/"StandardDicke_NormalPhase.csv")
SR=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG1"/"StandardDicke_Superradiant.csv")

NP_data=NP.to_numpy()
SR_data=SR.to_numpy()

plt.figure()
plt.plot(NP_data[:,0],NP_data[:,1],color='C0')
plt.plot(SR_data[:,0],SR_data[:,1],color='C0')



##### (d) nc=6 stable steady states
## Parameters
# wA=0.1
# wC=1
# kappa=1
MS=['.','s',"$o$",'v','d','*','p','v']
C=['C0','C1','C2','gold','C9','C4','C8','C5','C3','C7','C6']

plt.figure()
nc=6
X=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG1"/'StableSteadyStates_nc_{}.csv'.format(nc))
X_data=X.to_numpy()
n_col=np.shape(X_data)[1]
    
for i in range(n_col-1):
    if i==0:
        plt.plot(X_data[:,0],X_data[:,i+1],'.',color=C[i],markersize=5,label=X.columns[i+1],linestyle='None')
    else:
        plt.plot(X_data[:,0],X_data[:,i+1],marker=MS[i],color=C[i],markersize=3,label=X.columns[i+1],linestyle='None')

plt.xlabel=r'$g(\omega_C)$'
plt.ylabel=r'$Z_\mathrm{tot}/n_c$'


###### (e) Normal Phase Orbit
##parameters
#A=0.6
#g=0.1
# wA=0.1
# wC=1
# kappa=1

A=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG1"/'NormalPhase_Orbit_1.csv')
A_data=A.to_numpy()

b1=qt.Bloch()
b1.view=[-60,30]
b1.add_points([2*A_data[:,0],2*A_data[:,1],2*A_data[:,2]],meth='l')
b1.show()


###### (f) Superradiant Phase Orbit
##parameters
#A=0.6
#g=0.3
# wA=0.1
# wC=1
# kappa=1

B=pd.read_csv(Path(__file__).resolve().parent.parent/"DATA"/"FIG1"/'Superradiant_Orbit_1.csv')
B_data=B.to_numpy()

b2=qt.Bloch()
b2.view=[-60,30]
b2.add_points([2*B_data[:,0],2*B_data[:,1],2*B_data[:,2]],meth='l')
b2.show()

