# -*- coding: utf-8 -*-
"""
Created on Mon Aug 15 15:05:53 2016

@author: jbrody
"""

import serial
import numpy as np, matplotlib.pyplot as plt
from scipy.optimize import curve_fit, brentq
from sklearn.metrics import r2_score
ta=20#Change this to room temperature!
minutes=5
points=minutes*12
time1=np.zeros(points)
temp1=np.zeros(points)
adj1=np.zeros(points)
index=np.arange(points)
plt.ylim((0,25))
plt.xlabel('t (minutes)')
plt.ylabel('T (\xb0C)')
a=serial.Serial('COM12',9600)
for i in index:
    data=a.readline()
    splitdata=data.split()
    time1[i]=int(splitdata[0])/1000
    temp1[i]=float(splitdata[1])
    adj1[i]=temp1[i]+ta-temp1[0]#calibration to make initial measurement ta
    plt.plot(time1[i]/60,adj1[i],'bs')
    plt.pause(0.05)
a.close()
#%%Run this cell to use old data with new constants below
x1=.021#Change this to height of sensor!
l=.0254*6#Change this to the length of the rod!
side=.0254*.25
p=4*side#Change this to the perimeter of the rod!
a=side**2#Change this to the cross-sectional area of the rod!
rho=8030#Change this to density of the rod!
c=500#Change this to the specific heat of the rod!
terms=10;
mu=np.zeros(terms)
def transcendental(mu,k,h0):
    return 1/np.tan(mu)-k*mu/h0/l
def func(sec,alpha,m,h0):
    k=alpha*rho*c
    temp=ta-h0*ta*np.cosh(m*(l-x))/(h0*np.cosh(m*l)+m*k*np.sinh(m*l))#This is the steady-state temperature
    for n in np.arange(terms):
        mu[n]=brentq(transcendental,n*np.pi+1e-10,(n+1)*np.pi-1e-10,args=(k,h0))
    lambd=mu/l
    for i in np.arange(terms):
        temp=temp+2*ta*lambd[i]*h0/(lambd[i]*l+np.sin(lambd[i]*l)*np.cos(lambd[i]*l))/(m**2+lambd[i]**2)/(h0*np.cosh(m*l)+m*k*np.sinh(m*l))*(m*np.cos(lambd[i]*l)+lambd[i]*np.sin(lambd[i]*l))*np.sinh(m*l)*np.cos(lambd[i]*(l-x))*np.exp(-alpha*(m**2+lambd[i]**2)*sec)
    return temp
def fitdata(xin,time,temp):
    global x
    x=xin
    popt, pcov = curve_fit(func, time, temp,bounds=([1e-8,1,10],[1e-3,100,10000]),p0=[4e-6,10,100])
    fit=func(time,*popt)
    plt.plot(time/60,fit)
    perr = np.sqrt(np.diag(pcov))
    r2=r2_score(temp, fit)
    print("The fitting parameters are",popt)
    print("The uncertainties in the fitting parameters are",perr)
    print("The coefficient of determination is", r2) 
fitdata(x1,time1,adj1)