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()