# -*- coding: utf-8 -*-
"""
Created on Wed Aug 31 13:24:08 2016

@author: jbrody
"""

import numpy as np, matplotlib.pyplot as plt, pandas as pd
from scipy.optimize import curve_fit, brentq
from sklearn.metrics import r2_score
x=.021#Change this to height from bottom!
ta=20#Change this to room temperature!
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;
print ('Your time and temperature data should appear as columns in Excel\n')
input("Copy your time data to clipboard (CTRL-C), click to the right of the period at the end of this sentence, and press Enter.")
millisec=pd.read_clipboard(header=None)
sec=np.array(millisec).T[0]/1000
input("Copy your temperature data to clipboard (CTRL-C), click to the right of the period at the end of this sentence, and press Enter.")
rawtemp=pd.read_clipboard(header=None)
rawtemp=np.array(rawtemp).T[0]
meastemp=rawtemp+ta-rawtemp[0]#calibration to make initial measurement ta
mu=np.zeros(terms)
def transcendental(mu,k,h):
    return 1/np.tan(mu)-k*mu/h/l
def func(sec,alpha,m,h):
    k=alpha*rho*c
    temp=ta-h*ta*np.cosh(m*(l-x))/(h*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,h))
    lambd=mu/l
    for i in np.arange(terms):
        temp=temp+2*ta*lambd[i]*h/(lambd[i]*l+np.sin(lambd[i]*l)*np.cos(lambd[i]*l))/(m**2+lambd[i]**2)/(h*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
popt, pcov = curve_fit(func, sec, meastemp,bounds=([1e-8,1,10],[1e-3,100,10000]))
plt.plot(sec/60,meastemp,'o',label="measured data")
tempfit=func(sec,*popt)
plt.plot(sec/60,tempfit,label='theoretical fit ({0} terms)'.format(terms))
plt.xlabel('t (minutes)')
plt.ylabel('T (\xb0C)')
perr = np.sqrt(np.diag(pcov))
r2=r2_score(meastemp, tempfit)
print("The fitting parameters are",popt)
print("The uncertainties in the fitting parameters are",perr)
print("The coefficient of determination is", r2)
"""Uncomment to see only one, two, and three terms in sum
def funcplot(sec,alpha,m,h):
    k=alpha*rho*c
    temp=ta-h*ta*np.cosh(m*(l-x))/(h*np.cosh(m*l)+m*k*np.sinh(m*l))
    
    for n in np.arange(3):
        mu[n]=brentq(transcendental,n*np.pi+1e-10,(n+1)*np.pi-1e-10,args=(k,h))
    lambd=mu/l
    for i in np.arange(3):
        temp=temp+2*ta*lambd[i]*h/(lambd[i]*l+np.sin(lambd[i]*l)*np.cos(lambd[i]*l))/(m**2+lambd[i]**2)/(h*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)
        plt.plot(sec/60,temp,label='{0} terms'.format(i+1))
    return
funcplot(sec,*popt)"""    
plt.legend()