import matplotlib.pyplot as plt
import numpy as np
import scipy
from gwpy.timeseries import TimeSeries
from gwpy.plot import Plot

import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

###

import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (8,6)
mpl.rcParams['xtick.labelsize'] = 22
mpl.rcParams['ytick.labelsize'] = 22
mpl.rcParams['axes.grid'] = True
mpl.rcParams['grid.linestyle'] = ':'
mpl.rcParams['grid.color'] = 'grey'
mpl.rcParams['lines.linewidth'] = 2
mpl.rcParams['axes.labelsize'] = 24
mpl.rcParams['legend.handlelength'] = 3
mpl.rcParams['legend.fontsize'] = 22

from matplotlib import rc
rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
rc('text', usetex=True)

import seaborn as sns
sea = sns.color_palette("Set1")

###


## Reading data 

print("Getting H1 data for O4a")
# Dec 22, 2023 (165 MPc)
tstart = '2023-12-22 13:00:00'
tend = '2023-12-22 14:00:00'
h1_O4a = TimeSeries.get('H1:GDS-CALIB_STRAIN_CLEAN', tstart, tend, verbose=True)
h1asd_O4a = h1_O4a.asd(8,4)
h1asd_O4a.write('h1O4a_231222_1300_1400_asd.txt')

print("Getting L1 data for O4a")
# Dec 31, 2023 (165 MPc)
tstart = '2023-12-31 05:50:00'
tend = '2023-12-31 06:50:00'
l1_O4a = TimeSeries.get('L1:GDS-CALIB_STRAIN_CLEAN', tstart, tend, verbose=True)
l1asd_O4a = l1_O4a.asd(8,4)
l1asd_O4a.write('l1O4a_231231_0550_0650_asd.txt')

print("Getting H1 data for O3")
#H1, an hour without large glitches on March 19 2020
tstart = '2020-03-19 03:15:00'
tend = '2020-03-19 04:15:00'
h1_O3 = TimeSeries.fetch_open_data('H1', tstart, tend, sample_rate=16384)
h1asd_O3 = h1_O3.asd(8,4)
h1asd_O3.write('h1O3b_200319_0315_0414_asd.txt')

print("Getting L1 data for O3")
#L1, Jan 4 2020, an hour without glitches
tstart = '2020-01-04 00:00:00'
tend = '2020-01-04 01:00:00'
l1_O3 = TimeSeries.fetch_open_data('L1', tstart, tend, sample_rate=16384)
l1asd_O3 = l1_O3.asd(8,4)
l1asd_O3.write('l1O3b_200104_0000_0100_asd.txt')

# Load O3b asds
asd_H_O3 = TimeSeries.read('h1O3b_200319_0315_0414_asd.txt')
asd_L_O3 = TimeSeries.read('l1O3b_200104_0000_0100_asd.txt')

#Load O4a asds
asd_H_O4a = TimeSeries.read('h1O4a_231222_1300_1400_asd.txt')
asd_L_O4a = TimeSeries.read('l1O4a_231231_0550_0650_asd.txt')

def retrieve_data(timeseries):
    data = timeseries.value  
    times = timeseries.times.value  
    return times, data

times_h1O3b, data_h1O3b = retrieve_data(h1O3b)
times_l1O3b, data_l1O3b = retrieve_data(l1O3b)
times_h1O4a, data_h1O4a = retrieve_data(h1O4a)
times_l1O4a, data_l1O4a = retrieve_data(l1O4a)

asd_H_O4a_interp = scipy.interpolate.interp1d(times_h1O4a, data_h1O4a)(times_h1O3b)
asd_L_O4a_interp = scipy.interpolate.interp1d(times_l1O4a, data_l1O4a)(times_l1O3b)

fig, ax = plt.subplots(1,1)
ax.plot(times_h1O3b, data_h1O3b/asd_H_O4a_interp, label = 'LHO', color='salmon')
ax.plot(times_l1O3b, data_l1O3b/asd_L_O4a_interp, label = 'LLO', color='dodgerblue')
plt.xlabel('Frequency [Hz]')
plt.ylabel(r"ASD Ratio ( O3/O4a )")
ax.set_xlim(10,2000)
ax.set_ylim(1e-3,1e3)
ax.set_xscale('log')
ax.set_yscale('log')
ax.minorticks_on()
#plt.grid(True, which="both", ls="-")
plt.legend()
fig.tight_layout()
plt.savefig("../paper_figures/ASD_ratio.pdf", dpi = 500, bbox_inches = "tight")
plt.savefig("../paper_figures/ASD_ratio.png", dpi = 500, bbox_inches = "tight")

exit()