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

import importlib # this allows me to reload modules in ipython2 that I have changed since importing them the first time eg:  importlib.reload(lo_q)
#import subprocess
# need print function for newline-free printing
import matplotlib
# Force matplotlib to not use any Xwindows backend. NECESSARY FOR HEADLESS.
#matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Use gwpy to fetch data and translate times
import gwpy.timeseries
import gwpy.time
import numpy as np
import scipy.signal
import scipy.constants
import scipy.optimize 
#get this to get sensing function
import cross_corelation

def read_log_file(filename):
    f = open(filename)
    lines = f.readlines()
    logfile = []
    stop_times = np.array([])
    for line in lines[4:]:
        tok = line.split(" ")
        logfile.append([tok[0], tok[1], tok[2]])
        stop_times = np.append(stop_times, tok[2])
    f.close()
    return stop_times

struct_9high = {}
struct_9nom = {}
#since there was a bug in my script that prevented saving log file this time
struct_9high['stop times'] = 1266632350 + np.array([192, 400, 610, 820, 1030, 1350])
struct_9nom['stop times'] = 1266633764 + np.array([192, 400, 610, 820, 1030, 1350])
struct_9high['name'] = '9MHz high'
struct_9nom['name'] = '9MHz nominal'
scans_list = [struct_9high, struct_9nom]

### get the data
slow_chans = ['H1:CAL-CS_TDEP_KAPPA_C_OUTPUT', 'H1:OMC-READOUT_X0_OFFSET', 'H1:OMC-DCPD_SUM_OUT16']
for scan in scans_list:
    scan['kappa c'] = np.array([])
    scan['DCPD sum'] = np.array([])
    scan['x0'] = np.array([])
    scan['optical_gain'] = np.array([])
    for time in scan['stop times']:
        this_data = gwpy.timeseries.TimeSeriesDict.fetch(slow_chans,time-30, time, verbose=False)
        scan['kappa c'] = np.append(scan['kappa c'], np.median(this_data['H1:CAL-CS_TDEP_KAPPA_C_OUTPUT'].value))
        scan['DCPD sum'] = np.append(scan['DCPD sum'], np.median(this_data['H1:OMC-DCPD_SUM_OUT16'].value))
        scan['x0'] = np.append(scan['x0'], np.median(this_data['H1:OMC-READOUT_X0_OFFSET'].value)) 

#get sensing for a nominal time: 
nom_time = 1266624018;
nom_data = gwpy.timeseries.TimeSeriesDict.fetch(slow_chans,nom_time-30, nom_time, verbose=False)
nom_kappac = np.median(nom_data['H1:CAL-CS_TDEP_KAPPA_C_OUTPUT'].value)
nom_x0 = np.median(nom_data['H1:OMC-READOUT_X0_OFFSET'].value)
freq = np.array([40]); # get sensing at 40 Hz
[G, sensing_mA_per_meter] = cross_corelation.get_cal_model(freq, nom_time)

for scan in scans_list:
    #scale sensing by optical gain using pydarm model and kaapa c and x0  
    scan['optical gain'] = abs(sensing_mA_per_meter)*scan['kappa c']*scan['x0']/(nom_kappac*nom_x0)

##make a fit of normalized optical gain vs DCPD power
def DCPD_power(og, alpha, Ijunk):
    return Ijunk + (1/2*alpha)*og**2

for scan in scans_list:
    [scan['popt'], pcov] = scipy.optimize.curve_fit(DCPD_power, scan['optical gain'], scan['DCPD sum'])
    scan['perr'] = np.sqrt(np.diag(pcov))
og_vector = np.linspace(0,2e13, 100)

###### make plots
fig_w = 8  #figure size (for printing to pdf)
fig_h = 15
save_figs = True
File_tag = 'change 9MHz'
fig = plt.figure(figsize=(fig_w, fig_h))
ax = plt.subplot(1,1,1)
for scan in scans_list:
    popt = scan['popt']
    model_label = scan['name'] + ' ' + str(round(scan['popt'][1], 2)) + ' +/-' + str(round(scan['perr'][1], 2))+ ' mW of junk light'
    ax.plot(og_vector, DCPD_power(og_vector, *popt), label = model_label )
    ax.plot(scan['optical gain'], scan['DCPD sum'], 'o', label = scan['name'] + ' data')

ax.set_xlim([0, 7e12])
ax.set_ylim([0, 50])
ax.set_xlabel('optical gain [m/mA]')
ax.set_ylabel('DCPD power (mA)')
ax.legend()
if save_figs:
    fig.savefig(File_tag + 'optical_gain_vs_DARM_offset.pdf', bbox_inches='tight',format='pdf')
ax.loglog(freq, abs(struct_5pm['sensing scaled']), label = '5pm')
ax.loglog(freq, abs(struct_10pm['sensing scaled']), label = '10pm')
ax.loglog(freq, abs(struct_14pm['sensing scaled']), label = '14pm')
ax.set_ylim([1e12, 1e13])
ax.set_xlim([10, 1e4])
ax.legend()
ax.set_ylabel('1/rt Hz')
ax.set_xlabel('Frequency Hz')
fig.suptitle('sensing function')