#!/usr/bin/env python

import pydarm
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker

meas_dir = '/ligo/groups/cal/H1/measurements/'
timestamp = '20230705T192529Z'
IFO = 'H1'


act_L3_sl = meas_dir + f'SUSETMX_L3_SS/SUSETMX_L3_SS_{timestamp}.hdf5'
act_L2_sl = meas_dir + f'SUSETMX_L2_SS/SUSETMX_L2_SS_{timestamp}.hdf5'
act_L1x_sl = meas_dir + f'SUSETMX_L1_SS/SUSETMX_L1_SS_{timestamp}.hdf5'
# act_L1y_sl = meas_dir + f'SUSETMY_L1_SS/SUSETMY_L1_SS_{timestamp}.hdf5'
pcaly_file_sl = meas_dir + f'PCALY2DARM_SS/PCALY2DARM_SS_{timestamp}.hdf5'
# pcalx_file_sl = meas_dir + f'PCALX2DARM_SS/PCALX2DARM_SS_{timestamp}.hdf5'
darm_file_sl = meas_dir + f'DARMOLG_SS/DARMOLG_SS_{timestamp}.hdf5'


act_L3_xml = meas_dir + 'SUSETMX_L3_SS/SUSETMX_L3_SS_20230517T163625Z.xml'
act_L2_xml = meas_dir + 'SUSETMX_L2_SS/SUSETMX_L2_SS_20230517T161119Z.xml'
act_L1x_xml = meas_dir + 'SUSETMX_L1_SS/SUSETMX_L1_SS_20230517T154827Z.xml'
# act_L1y_xml = meas_dir + 'SUSETMY_L1_SS/SUSETMY_L1_SS_20230519T153633Z.xml'
pcaly_file_xml = meas_dir + 'PCALY2DARM_SS/PCALY2DARM_SS_20230621T213239Z.xml'
# pcalx_file_xml = meas_dir + 'PCALX2DARM_SS/PCALY2DARM_SS_20230519T152218Z.xml'
darm_file_xml = meas_dir + 'DARMOLG_SS/DARMOLG_SS_20230621T211522Z.xml'


L3_sl = pydarm.measurement.Measurement(act_L3_sl)
L2_sl = pydarm.measurement.Measurement(act_L2_sl)
L1x_sl = pydarm.measurement.Measurement(act_L1x_sl)
# L1y_sl = pydarm.measurement.Measurement(act_L1y_sl)
pcaly_sl = pydarm.measurement.Measurement(pcaly_file_sl)
# pcalx_sl = pydarm.measurement.Measurement(pcalx_file_sl)
darm_sl = pydarm.measurement.Measurement(darm_file_sl)

L3_xml = pydarm.measurement.Measurement(act_L3_xml)
L2_xml = pydarm.measurement.Measurement(act_L2_xml)
L1x_xml = pydarm.measurement.Measurement(act_L1x_xml)
# L1y_xml = pydarm.measurement.Measurement(act_L1y_xml)
pcaly_xml = pydarm.measurement.Measurement(pcaly_file_xml)
# pcalx_xml = pydarm.measurement.Measurement(pcalx_file_xml)
darm_xml = pydarm.measurement.Measurement(darm_file_xml)

for stage in ['L3', 'L2', 'L1x', 'pcaly', 'darm']: #
    if (stage == 'L2') or (stage == 'L3'):
        print(f'stage: {stage}')
        for meas in ['sl','xml']:
            exec(f'tf_{meas} = {stage}_{meas}.get_raw_tf("{IFO}:SUS-ETMX_{stage}_CAL_EXC","{IFO}:LSC-DARM_IN1_DQ")')
    elif (stage == 'L1x'):
        print(f'stage: {stage}')
        for meas in ['sl','xml']:
            exec(f'tf_{meas} = {stage}_{meas}.get_raw_tf("{IFO}:SUS-ETMX_L1_CAL_EXC","{IFO}:LSC-DARM_IN1_DQ")')
    elif (stage == 'L1y'):
        print(f'stage: {stage}')
        for meas in ['sl','xml']:
            exec(f'tf_{meas} = {stage}_{meas}.get_raw_tf("{IFO}:SUS-ETMY_L1_CAL_EXC","{IFO}:LSC-DARM_IN1_DQ")')
    elif (stage == 'pcalx'):
        print(f'stage: {stage}')
        for meas in ['sl','xml']:
            exec(f'tf_{meas} = {stage}_{meas}.get_raw_tf("{IFO}:CAL-PCALX_RX_PD_OUT","{IFO}:LSC-DARM_IN1_DQ")')
    elif (stage == 'pcaly'):
        continue
        print(f'stage: {stage}')
        for meas in ['sl','xml']:
            exec(f'tf_{meas} = {stage}_{meas}.get_raw_tf("{IFO}:CAL-PCALY_RX_PD_OUT_DQ","{IFO}:LSC-DARM_IN1_DQ")')
    elif (stage == 'darm'):
        # continue
        print(f'stage: {stage}')
        for meas in ['sl','xml']:
            exec(f'tf_{meas} = {stage}_{meas}.get_raw_tf("{IFO}:LSC-DARM1_EXC","{IFO}:LSC-DARM1_IN2")')

    
    fig,axs = plt.subplots(2,1)
    fig.suptitle(f'{stage} sweep, Uncertainties Comparison')
    axs[0].semilogx(tf_sl[0],tf_sl[3],label=f'{stage} Simulines', marker='.')
    axs[0].semilogx(tf_xml[0],tf_xml[3],label=f'{stage} DTT', marker='.')
    axs[0].legend()
    axs[0].set(ylabel='Uncertainty, 1 sigma')
    axs[0].set(xlabel='Frequency [Hz]')
    axs[0].set_ylim(0,0.1)
    axs[0].set_xlim(5,1200)
    
    axs[1].loglog(tf_sl[0],abs(tf_sl[3]/tf_xml[3]),label='SL/DTT', marker='.')
    axs[1].legend()
    axs[1].set(ylabel='SL/DTT Uncert Ratio')
    axs[1].set_ylim(0.02,50)
    axs[1].set_xlim(5,1200)
    axs[1].get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
    axs[1].set_yticks([0.02,0.1,0.2,0.5,1,2,5,10,50])
    axs[1].tick_params(axis='y', which='both', labelsize=10)
    # plt.show()
    plt.savefig(f'./{stage}_uncert_compare.png',bbox_inches='tight', dpi=300)