# analysis od the A1 JPL PD diode
# Aidan Brooks - 10-Sept-2021

import cdsutils
import numpy as np
import matplotlib.pyplot as plt
import os, glob
import scipy.signal


# function to get slow channels from frames
def getSlowChannels(QEstart):

    # calibration of channels
    REFcal = 0.862E-3 # https://nodus.ligo.caltech.edu:8081/QIL/2651
    JPLcal = 0.5/1000 # factor of 2 in filter and 1000 from preamp
    
    # get the slow data from frames
    # [REF PD, JPL PD, PD BIAS, LASER SWITCH']
    chans = ['C4:TST-FM29_OUT16',   # REFPD
             'C4:TST-FM30_OUT16',   # JPL PD
             'C4:TST-FM13_OFFSET',  # PD BIAS
             'C4:TST-FM12_OFFSET']  # LASER SW
    duration = 156
    dataSlow = cdsutils.getdata(chans, duration, QEstart+1)
    
    # get the time
    time = np.arange(0, duration, 1/np.float(16))
    REFPD =  REFcal*(dataSlow[0].data - np.min(dataSlow[0].data))
    JPLPD = JPLcal*dataSlow[1].data
    bias = dataSlow[2].data
    laserSW = dataSlow[3].data

    return REFPD, JPLPD, bias, laserSW



# get the Keithley and conditions files from the current measurement directory
def getDirTextFiles(dirIN):
    condFile = dirIN + '_conditions.txt'
    dataCond = np.genfromtxt(condFile, comments='PD_gain') # skip the PD_gain line
    tempList = dataCond[0:4,1]
    specStart = dataCond[4,1]
    QEstart = dataCond[5,1]
    
    
    # get the Keithley filenames
    darkFile = dirIN + '_darkcurrent.txt'
    photoLowFile = dirIN + '_photocurrent_20mA.txt'
    photoFullFile = dirIN + '_photocurrent_full.txt'
    
    # get the Keithley data
    dataDark = np.loadtxt(darkFile)
    dataLow = np.loadtxt(photoLowFile)
    dataFull = np.loadtxt(photoFullFile)

    return dataDark, dataLow, dataFull, QEstart, specStart, tempList


# function to calculate the QE
def getQEFromPDTimeSeries(REFPD, JPLPD):

    #-------------------------------------------
    # constants
    h = 6.63E-34
    c = 3.0E8
    e = 1.6E-19
    L = 2004E-9
    QEsc = (h*c)/(L*e)
    #-------------------------------------------
    
    # loop through biases for QE measurement
    biasList = [-100, 0, 300, 600, 800, 1000]
    QEmean = np.array([])
    for jj in range(0, len(biasList)):
        # get the indices where results are at the correct bias level
        ind0 = np.argwhere(bias == biasList[jj])
        # subtract off the zero light voltage of the PD
        minPD = np.mean(JPLPD[ind0[-160:-128]])
        PD = JPLPD[ind0]-minPD
        REF = REFPD[ind0]

        # get values where REFPD level is greater than 5% of maximum
        ind1 = np.argwhere(REF>(0.03*np.max(REF)))
        ind1 = ind1[:,0]
        #print(ind1)
        #print(PD[ind1].shape)
        #print(REF[ind1].shape)
        QE = QEsc*PD[ind1]/REF[ind1]
        QEmean = np.append(QEmean, np.mean(QE[4:12]))

    return QEmean


# function to get the dark noise spectra
def getDarkNoiseSpectra(specStart):

    # AC channel data calibration
    countsToVolts = 20/np.float(65536)
    SR560 = 1E3  # gain of 1000
    preamp = 1E3 # transimpedance of 1000 V/A
    ACscale = countsToVolts/(SR560*preamp)   # count to amps
    
    # get the spectra data
    chanSpec = 'C4:TST-FM31_OUT_DQ'
    duration = 300
    dataSpec = cdsutils.getdata(chanSpec, duration, specStart+1)
    dataBias = cdsutils.getdata('C4:TST-FM13_OFFSET', duration, specStart+1)
    
    # get the times
    time0 = np.arange(0,duration,1/np.float(16384))
    time1 = np.arange(0,duration,1/np.float(16))

    biasList = [-100, 0, 300, 600, 1000]
    t0 = 0
    for jj in range(0, 5):
        # work out times where bias is desired value
        # (BIAS = CORRECT VALUE) AND (TIME > PREVIOUS BIAS LEVEL)
        ind0 = np.argwhere( (dataBias.data == biasList[jj]) * (time1>t0))
        timeSub = time1[ind0]
        t0 = np.min(timeSub)
        t1 = np.max(timeSub)
        # get corresponding PD data (when PD data is within same time frame)
        ind1 = np.argwhere( (time0 >= t0)*(time0 < t1) )
        dataACSub = ACscale*dataSpec.data[ind1]
        # bypass the first 10s
        dataSub = dataACSub[10*16384:]
        dataSub = dataSub.reshape(len(dataSub)) # reshape into a row for WELCH            
        f, Pxx = scipy.signal.welch(dataSub, dataSpec.sample_rate,
                                    window='hann', nperseg=65536/4,
                                    detrend='linear')


    return f, Pxx

#******************************************************************




#----------------------------------------------
# get the files
files = sorted(glob.glob('../A1/131*-1*'))
print(len(files))

arr = np.array([[1,1,1,1, 1]])
QEList = np.array([])
# go through all the measurements
#for ii in range(0, len(files)):
for ii in range(0, 10):
    print(ii)

    dirIN = files[ii] + '/'
    t0 = dirIN[6:16]
    dirIN = dirIN + t0

    try:
        # import data from the stored text files
        dataDark, dataLow, dataFull, QEstart, specStart, tempList = getDirTextFiles(dirIN)
        
        # get the slow data from frames
        REFPD, JPLPD, bias, laserSW = getSlowChannels(QEstart)

        # extract the QE from the slow data time series
        QEmean = getQEFromPDTimeSeries(REFPD, JPLPD)

        # get the data noise spectra
        f, Pxx = getDarkNoiseSpectra(specStart)
        plt.loglog(f, np.sqrt(Pxx))
                
        print(dirIN)
 
    except Exception as e:
        print(e)


print(f.shape)
print(Pxx.shape)
plt.grid(True)
plt.show()