from scipy.optimize import curve_fit
from scipy.signal import decimate
import cdsutils as cds
import numpy as np
import matplotlib.pyplot as plt

#extract data
channel = ['C1:PEM-SEIS_EX_TEMP_OUT_DQ']
tstart = 1200962230
tend = 1201041084
data = cds.getdata(channel,tend-tstart,tstart)

temp = data[0].data

time = []
for i in range(len(temp)):
    time.append(i)
    
temp_new = []
for i in range(len(temp)):
    if i % 100 == 0:
        temp_new.append(temp[i])

time_new = []
for i in range(len(time)):
    if i % 100 == 0:
        time_new.append(time[i])

time_new1 = np.linspace(0,len(temp)/32.,(len(data[0].data+1)))*3.125/3600
time_new1 = time_new1[::100]/32.

#fit an exponential to the data
def func(x,a,b,c):
    return a*np.exp(-b*x)+c

popt, pcov = curve_fit(func,np.array(time_new1),np.array(temp_new),p0=(-500,0.3,-4800))

#errors
errors = temp_new - func(time_new1, *popt)

#filtering and fitting the data, along with approximate calibration:
y = decimate(np.array(temp_new)*-1./198.4,4,ftype='fir')
x = []
for i in range(len(y)):
    x.append(i)
x = np.array(x[8:])*(3.125/3600./3.)
y = np.array(y[8:])
popt1, pcov1 = curve_fit(func, x, y, p0=(1,0.3,1))
errors_1 = y - func(x, *popt1)

f, (ax1, ax2) = plt.subplots(2)
ax1.plot(x, y, 'b.', label = 'Measured')
ax1.plot(x, func(x, *popt1), 'r', label = 'Fit')
ax1.legend()
ax1.set(title = 'Cooling Fit', ylabel = 'Temp(C)')
ax2.plot(x, errors_1)
ax2.set(title = 'Errors', xlabel = 'Time (hr)', ylabel = 'Temp (C)')
plt.show()