#!/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(50*(1e-3)*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(50*(1e-3)*10**(spec_dark/20))
#create a filter to immitate D1700376, mulitply spectrum by 1/mag
w_vector = np.array(freq_nosqz)*2*np.pi  #np.logspace(-1, 7, 500)*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)
w, mag_TI_3MHz_dB, phase_TI = signal.bode(omc_ti, 2*np.pi*3.125e6)
mag_TI = 10**(mag_TI_dB/20)
mag_TI_3MHz = 10**(mag_TI_3MHz_dB/20)
#cal_spec_sqz = spec_sqz -mag_split_dB - mag_TI_dB
spec_ArtHz_nosqz = spec_V_rtHz/(mag_split*mag_TI)
spec_ArtHz_dark = spec_dark_V_rtHz/(mag_split*mag_TI)

quad_diff = np.sqrt(abs(spec_ArtHz_nosqz**2 - spec_ArtHz_dark**2))
#also esimate shot noise
electron_charge = 1.6e-19  #columbs
shot_noise = np.sqrt(2*electron_charge*10e-3)#10mA DC photocurrent
#what kind of amplitude noise around 3MHz can we handle?
#estimate the amplitude of CLF in Nov 2019 (we saw some improvement when we reduced this), there was -20dBm at demod, before 20dB gain from splitter chasis
responsivity = 0.858#A/W
A_3MHz = np.sqrt(50*1e-3*10**(-40/20))/(mag_TI_3MHz)
E_carrier = np.sqrt(20e-3/responsivity) #proportional to the carrier amplitude
E_clf = A_3MHz/(E_carrier*responsivity)
noise_we_need_to_measure_now = shot_noise*abs(E_carrier/E_clf)
noise_we_need_to_measure_filter_cavity = shot_noise*abs(E_carrier/(E_clf*np.sqrt(20)))#20 times more CLF power for filter cavity
noise_to_equal_shot = noise_we_need_to_measure_now*np.ones(800)
### 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_ArtHz_nosqz, label = 'squeezer blocked')
ax.semilogy(freq_nosqz, spec_ArtHz_dark, label = 'dark')
ax.semilogy(freq_nosqz, noise_to_equal_shot, label = 'amplitude noise to downconvert to shot noise level')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('A/rtHz')
ax.set_ylim([1e-9, 1e-5])
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.loglog(w, mag_split, label ='D1700376')
ax.loglog(w, mag_TI, label ='OMC TI')
ax.set_xlim([0.1,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([0.1,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()