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

#for reading in files
import os
#for making plots
import matplotlib
import matplotlib.pyplot as plt
#for filters
from scipy import signal
import numpy as np

def open_data(filename):
    f = open(filename)
    lines = f.readlines()
    freq = np.array([])
    spec = np.array([])
    for line in lines[16:816]:  #cut out header lines
            tok = line.split("\t")
            #print(tok)
            freq = np.append(freq, float(tok[0]))
            spec = np.append(spec, float(tok[1]))
    return freq , spec

freq_dark, spec_dark = open_data('CC002.TXT')#CC001 is dark noise from PD1 connector on D1700367, CC002 is from PD2
freq_nosqz, spec_nosqz = open_data('BB001.TXT')
#convert units
spec_V_rtHz = np.sqrt((1e-3/50)*10**(spec_nosqz/20))#50 Ohms, converting dBm to V^2/Hz, take square root to get ASD units (V/rtHz)
spec_dark_V_rtHz = np.sqrt((1e-3/50)*10**(spec_dark/20))
#create a filter to immitate D1700376, mulitply spectrum by 1/mag
w_vector = np.array(freq_nosqz)*2*np.pi
filt = signal.ZerosPolesGain([100,1],[-290e3, -265e3, -60e6],[10*60e6])
w, mag_split_dB, phase_split = signal.bode(filt, w_vector)
mag_split = 10**(mag_split_dB/20)

#from Lee's alog 54651, a fit to OMC DCPD transimpendance
omc_ti = signal.ZerosPolesGain(np.array([-1.09367517e+06+9.16935186e+04j, -1.09367517e+06-9.16935186e+04j,
       -4.54321270e+01+3.94631525e-01j, -4.54321270e+01-3.94631525e-01j,
       -2.85835003e+07+0.00000000e+00j]), np.array([-1.12584306e+07+1.59926553e+07j, -1.12584306e+07-1.59926553e+07j,
       -3.11667017e+07+1.49827313e+08j, -3.11667017e+07-1.49827313e+08j,
       -4.41594538e+07+5.59042545e+07j, -4.41594538e+07-5.59042545e+07j,
       -1.51046460e+07+2.52263715e+07j, -1.51046460e+07-2.52263715e+07j,
       -1.02822684e+05+0.00000000e+00j, -9.54273138e+04+0.00000000e+00j,
       -5.08417603e+02+0.00000000e+00j, -4.91824766e+02+0.00000000e+00j]), 2.71124765615596e+56)

w, mag_TI_dB, phase_TI = signal.bode(omc_ti, w_vector)
mag_TI = 10**(mag_TI_dB/20)

#cal_spec_sqz = spec_sqz -mag_split_dB - mag_TI_dB
spec_mArtHz_nosqz = spec_V_rtHz/(mag_split*mag_TI)
spec_mArtHz_dark = spec_dark_V_rtHz/(mag_split*mag_TI)

quad_diff = np.sqrt(spec_mArtHz_nosqz**2 - spec_mArtHz_dark**2)
#also esimate shot noise
electron_charge = 1.6e-19  #columbs
shot_noise = np.sqrt(2*electron_charge*10e-3)#10mA DC photocurrent

### make figures
fig_w = 12  #figure size (for printing to pdf)
fig_h = 8

fig = plt.figure(figsize=(fig_w, fig_h))
ax = plt.subplot(1,1,1)
#ax.plot(freq_sqz, spec_sqz, label = 'squeezer on')
ax.semilogy(freq_nosqz, spec_mArtHz_nosqz, label = 'squeezer blocked')
ax.semilogy(freq_nosqz, spec_mArtHz_dark, label = 'dark')
ax.semilogy(freq_nosqz, quad_diff, label = 'quadrature difference')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('mA/rtHz')
ax.set_ylim([1e-10, 1e-8])
ax.grid()
ax.legend()
fig.savefig('DCPD_4MHz.pdf', bbox_inches='tight',format='pdf', dpi=fig.dpi)



fig = plt.figure(figsize=(fig_w, fig_h))
ax = plt.subplot(2,1,1)
ax.semilogx(w, mag_split_dB, label ='D1700376')
ax.semilogx(w, mag_TI_dB, label ='OMC TI')
ax.set_xlim([100,100e6])
ax.grid()
ax.legend()
ax = plt.subplot(2,1,2)
ax.semilogx(w, phase_split, label ='D1700376')
ax.semilogx(w, phase_TI, label ='OMC TI')
ax.set_xlim([100,100e6])
ax.grid()

fig = plt.figure(figsize=(fig_w, fig_h))
ax = plt.subplot(1,1,1)
ax.plot(freq_sqz, cal_spec_sqz, label = 'squeezer on')
ax.plot(freq_nosqz, cal_spec_nosqz, label = 'squeezer blocked')
#ax.plot(freq_nosqz, spec_nosqz, label = 'squeezer blocked uncalibrated')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('dBm with split off filter calibrated out')
ax.set_ylim([-170, -80])
ax.grid()
ax.legend()
plt.show()