#!/usr/bin/env python
# -*- coding: utf-8 -*-

#before running this, in command line : export NDS2_CLIENT_GAP_HANDLER=STATIC_HANDLER_NEG_INF


import gwpy.timeseries
import numpy as np
import matplotlib.pyplot as plt
import os
from scipy.optimize import minimize
#import h5py

os.environ["NDS2_CLIENT_GAP_HANDLER"] = "STATIC_HANDLER_NEG_INF"


chans = ['H1:IOO-OFI_THERMISTOR_2_TEMPERATURE.mean,m-trend','H1:IOO-OFI_PD_A_DC_POWER.mean,m-trend' , 'H1:ASC-AS_C_NSUM_OUT16.mean,m-trend' ]#, 'L1:GRD-ISC_LOCK_STATE_N.min,m-trend']
gpsStart = 1446844666
gpsStop = 1446864156
data = gwpy.timeseries.TimeSeriesDict.fetch(chans,gpsStart, gpsStop, verbose=True, pad=float("-inf"))

#OFI PD, calibrated into mW (1e-3) is behind 1% transmission beam splitter
ratio = data['H1:IOO-OFI_PD_A_DC_POWER.mean,m-trend'].value*100*1e-3/ data['H1:ASC-AS_C_NSUM_OUT16.mean,m-trend'].value

#want to take an average once the temp has settled. 
#temp_hist = np.histogram(data['H1:IOO-OFI_THERMISTOR_2_TEMPERATURE.mean,m-trend'].value)
#manually pick out times we want averages for, times just before a temp step
avg_stop_idx = [48, 108, 168, 228, 287, 324]
avg_temp = np.array([])
avg_ratio = np.array([])
temp_vector = data['H1:IOO-OFI_THERMISTOR_2_TEMPERATURE.mean,m-trend'].value
for step in avg_stop_idx:
    #average for 10 minutes before you hand picked times above
    avg_temp = np.append(avg_temp, np.mean(temp_vector[step-10:step]))
    avg_ratio =  np.append(avg_ratio, np.mean(ratio[step-10:step]))

def leakage_power(mu, temp):
    #construct a model to fit to the data, this is a model of the ratio of power heading to HAM7 over HAM6 (1 would be all the power)
    constant_leakage, angle_coeff , best_temp = mu
    return constant_leakage + (np.sin((np.pi/180)*angle_coeff* (temp- best_temp)))**2

def fit_leakage(meas, guess = [0.0033, 0.2, 27.5]):
    temp , ratio = meas

    def cost(mu):
        return np.sum((leakage_power(mu, temp) - ratio)**2)
    
    res = minimize(
            cost,
            x0=guess,
            method='Nelder-Mead',
            options={"maxiter":1500},
            )
    return res.x

fit_params = fit_leakage([avg_temp, avg_ratio])

fig, [ax1,ax2, ax3, ax4]  = plt.subplots(4,1,figsize = [10,8])
ax1.plot(data['H1:IOO-OFI_THERMISTOR_2_TEMPERATURE.mean,m-trend'].value , label = f'OFI thermistor 2 [C]')
ax2.plot(data['H1:IOO-OFI_PD_A_DC_POWER.mean,m-trend'].value , label = f'OFI PD A power [mW]')
ax3.plot(data['H1:ASC-AS_C_NSUM_OUT16.mean,m-trend'].value , label = f'AS port power arriving in HAM6[W]')
ax4.plot(ratio*100, label ='% of power towards HAM7/ power towards HAM6')


for ax in [ax1,ax2,ax3,ax4]:
    ax.legend()
ax4.set_xlabel('time [minutes]')
plt.savefig('OFI_temp_changes.png')

temp = np.linspace(20,35,100)
plt_title = f'Fit parameters: {fit_params[0]*100:0.2f}% constant leakage,{fit_params[1]:0.3f} degrees rotation / degrees C, \n'+\
     f'{fit_params[2]:0.1f} degrees C minimizes leakage'

fig, ax = plt.subplots(1,1,figsize = [10,8])
ax.plot(avg_temp, avg_ratio*100,'o', label = 'average of last 10 minutes of data for each step')
ax.plot(temp, 100*leakage_power(fit_params,temp), label = f'fit model')
ax.set_xlabel('OFI temperature [C]')
ax.set_ylabel('% of power to HAM7')
ax.legend()
ax.set_title(plt_title, fontsize='medium')
plt.savefig('ratio_vs_temp.png')