# -*- coding: utf-8 -*-
import os
import matplotlib.pyplot as plt
plt.ion()
from scipy.optimize import curve_fit
import numpy as np
import code
import re
import multiprocessing
import time
import easygui

version = 1.3
print("Fajho v"+str(version))
print("Run man() to see manual")

fig = None
cid = None
eventX = 0
eventY = 0
eventReceived = False

x,y,heater,Uarr = [],[],[],[]
Ts = []
bTs = False

bBase = False
base_a, base_b = 0,0
base_start, base_end = 0,0
Dbase_a = 0
bMain = False
main_a, main_b = 0,0
Dmain_a = 0
main_start, main_end = 0,0
main_a_corr = 0
Dmain_a_corr = 0
bExp = False
exp_A, exp_b, exp_c = 0,0,0
exp_start, exp_end = 0,0
Dexp_b = 0
bInt = False
bMainCorrected = False

bMainExp = False
mainexp_A, mainexp_b, mainexp_c = 0,0,0
Dmainexp_b = 0
mainexp_start, mainexp_end = 0,0

bBaseExp = False
baseexp_A, baseexp_b, baseexp_c = 0,0,0


bCalib = False
Cp = 0
dCp = 0
alfa = 0
dalfa = 0

bEpszilonVesszo = False
ev = 0

directory = ""


autoshow = True
grid = True
scaletodata = True
dualfit = True

R = 4.60
DR = 0.01


bLegend = False
def legend():
    global bLegend
    bLegend = not bLegend
    if(autoshow):
        show()

def man():
    # work in progress....
    print("Dokumentáció (nem teljes):")
    print("- Függvények:")
    print("    - Fájl megnyitás:")
    print("        new()")
    print("    - Könyvtár váltás:")
    print("        cd()")
    print("")
    print("    - Ábrázolás:")
    print("        s()    [alias: show]")
    print("          - Betöltött adatsor ábrázolása")
    print("        legend()")
    print("          - Ábrafelirat ki/be kapcsolása")
    print("")
    print("    - Illesztések:")
    print("        esz(xmin,xmax)    [alias: eloszakasz]")
    print("          - Előszakasz illesztése xmin és xmax között.")
    print("        fsz(xmin,xmax)    [alias: foszakasz]")
    print("          - Lineáris főszakasz illesztése xmin és xmax között.")
    print("        kfsz(xmin,xmax)")
    print("          - Lineáris főszakasz illesztése xmin és xmax között a korrigált hőmérséklet (T*) görbére.")
    print("        efsz(xmin,xmax)   [alias: expfoszakasz]")
    print("          - [Beejtéses mérésre]. Exponenciális főszakasz illesztése xmin és xmax között. Az eredmény epszilon' paramétert kimenti fájlba. [Feltételek: korrigált hőmérséklet (T*)]")
    print("        usz(xmin,xmax)    [alias: utoszakasz]")
    print("          - Utószakasz exponenciális illesztése xmin és xmax között. Ha a 'dualfit' globális konstans True, akkor automatikusan végrehajtja a df()-et is.")
    print("        df(limit='fixed') ")
    print("          - Kettős illesztés, az előszakasz exponenciális korrekciójához. Ha a limit paraméter 'fixed', akkor az előszakasz számára nem szabad paraméter az (utószakaszból vett) limesz, ha 'free' akkor ez is szabad illesztési paraméter.")
    print("")
    print("        Minden illesztés hívható paraméterek nélkül is, például: esz(). Ekkor az inputról fogja kérni a program a határokat. Minden illesztés eredménye elmentődik a fitlogba.")
    print("")
    print("    - Számolások:")
    print("        korrigaltT()")
    print("          - Korrigált hőmérséklet (T*) számolása. Feltételek: előszakasz, utószakasz.")
    print("")
    print("        - Vízérték mérés adatsora:")
    print("            calib()")
    print("              - Vízérték és hőátadási paraméter számolása (kalibráció). Az eredményeket kimenti fájlba, melyet a program minden új fájlbeolvasásnál megpróbál betölteni. [Feltételek: előszakasz, főszakasz, utószakasz.]")
    print("")
    print("        - Beejtéses mérés adatsora:")
    print("            integral(xmin,xmax)")
    print("              - Fajhő számolása xmin és xmax közti integrálból. [Feltételek: kalibráció, előszakasz.]")
    print("            si(xmin,xmax_low, xmax_high)   [alias: serialintegrate]")
    print("              - Az integrál felső határának variálása, xmax_low és xmax_high között minden értékre integrált és fajhőt számol melyet aztán ábrázol a felső határ függvényében. [Feltételek: kalibráció, előszakasz.]")
    print("        - Ráfűtéses mérés adatsora:")
    print("            rafut()")
    print("              - Fajhő számolása a ráfűtés főszakaszából. [Feltételek: kalibráció, előszakasz, főszakasz]")
    print("            rafutteljes(Ts_start)  ")
    print('              - Fajhő számolása a ráfűtés főszakaszából és utószakaszából. ("A régi módszer.") A Ts_start paraméter a T* korrigált hőmérséklet átlagolásának kezdőpontja. [Feltételek: kalibráció, előszakasz, főszakasz, utószakasz, korrigált hőmérséklet, korábbi beejtéses mérés efsz() illesztése]')





def flin(x, a, b):
    return a*x+b

def fexp(x,A,b,c):
    return A*np.exp(-b*x)+c


def do_fit(x,y,xmin,xmax, func = flin):
    xf = []
    yf = []
    for i in range(len(x)):
        if(xmin <= x[i] <= xmax):
            xf.append(x[i])
            yf.append(y[i])


    return curve_fit(func, xf, yf)


def getK():
    if(not bEpszilonVesszo or not bCalib):
        return False
    return ev/(ev - alfa/Cp)




def fitreport(params, covs):
    perr = np.sqrt(np.diag(covs))
    print("Final set of parameters            Asymptotic Standard Error: ")
    print("=======================            ==========================")
    for i in range(len(params)):
        print("Param-"+str(i)+":   "+str(params[i])+"      +/-   "+str(perr[i])+"      ("+ f'{abs(perr[i]/params[i]):.3e}'+"%)")
    with open(directory+"OUTfajho_fitlog.dat", "a") as f:
        print("Final set of parameters            Asymptotic Standard Error: ",file=f)
        print("=======================            ==========================",file=f)
        for i in range(len(params)):
            print("Param-"+str(i)+":   "+str(params[i])+"      +/-   "+str(perr[i])+"      ("+ f'{abs(perr[i]/params[i]):.3e}'+"%)",file=f)




def onclick(event):
    global eventX, eventY, eventReceived
    global fig, cid
    
    if(event.dblclick):
        eventX = event.xdata
        eventY = event.ydata
        fig.canvas.mpl_disconnect(cid)
        eventReceived = True

        #endEvent("click")
        return True
    return False


def rangeselectorKeyboard(text):
    global eventX, eventY, eventReceived
    value = input(text)

    eventX = value
    eventReceived = True

    #endEvent("keyboard")

def rangeselectorMouseActivate():
    global cid
    global fig
    cid = fig.canvas.mpl_connect('button_press_event', onclick)


#def endEvent(source):
    #if(source == "click"):



def createRangeEvent(text):
    global eventReceived

    import threading
    rangeselectorMouseActivate()
    proc = threading.Thread(target=rangeselectorKeyboard, args=(text,))
    #proc = multiprocessing.Process(target=rangeselectorKeyboard, args=(text,))
    #proc = multiprocessing.Process(target=rangeselectorKeyboard, args=(text))
    proc.start()


    while(not eventReceived):
        time.sleep(0.1)

    if(proc.is_alive()):
        proc.terminate()
    
    xmin = eventX
    eventReceived = False
    return xmin




def fitbase(*args):
    global bBase
    global base_a
    global Dbase_a
    global base_b
    global base_start
    global base_end

    if(len(args) == 2):
        xmin = args[0]
        xmax = args[1]
    else:
        xmin = float(input("Fit start: "))
        xmax = float(input("Fit end: "))
    
    base_start = xmin
    base_end = xmax
    params, covs = do_fit(x,y,xmin,xmax)
    log("\nBase fit in x["+str(xmin)+"; "+str(xmax)+"],   y = p0*x + p1")
    fitreport(params, covs)

    bBase = True
    base_a = params[0]
    base_b = params[1]

    perr = np.sqrt(np.diag(covs))
    Dbase_a = perr[0]

    if(autoshow):
        show()
fb = fitbase
eloszakasz = fitbase
esz = fitbase



def fitmain(*args):
    global bMain
    global main_a
    global Dmain_a
    global main_b
    global main_start
    global main_end
    global main_a_corr
    global Dmain_a_corr

    if(not bBase):
        print("Error: missing base fit")
        return

    if(bMainCorrected):
        sure = input("Are you sure? (yes/no): ")
        if(sure != "yes"):
            return

    if(len(args) == 2):
        xmin = args[0]
        xmax = args[1]
    else:
        xmin = float(input("Fit start: "))
        xmax = float(input("Fit end: "))
    
    main_start = xmin
    main_end = xmax
    params, covs = do_fit(x,y,xmin,xmax)
    log("\nMain fit in x["+str(xmin)+"; "+str(xmax)+"],   y = p0*x + p1")
    fitreport(params, covs)

    bMain = True
    main_a = params[0]
    main_b = params[1]

    perr = np.sqrt(np.diag(covs))
    Dmain_a = perr[0]

    if(bBase):
        main_a_corr = main_a - base_a
        Dmain_a_corr = Dmain_a + Dbase_a


    if(autoshow):
        show()
fm = fitmain
foszakasz = fitmain
fsz = fitmain



def fitmaincorrected(*args):
    global bMainCorrected
    global bMain
    global main_a
    global Dmain_a
    global main_b
    global main_start
    global main_end
    global main_a_corr
    global Dmain_a_corr

    if(not bTs):
        print("Error: missing T* fit")
        return

    if(len(args) == 2):
        xmin = args[0]
        xmax = args[1]
    else:
        xmin = float(input("Fit start: "))
        xmax = float(input("Fit end: "))
    
    main_start = xmin
    main_end = xmax
    params, covs = do_fit(x,Ts,xmin,xmax)
    log("\nMain fit in x["+str(xmin)+"; "+str(xmax)+"] onto T*,   y = p0*x + p1")
    fitreport(params, covs)

    bMain = True
    bMainCorrected = True
    main_a = params[0]
    main_b = params[1]

    perr = np.sqrt(np.diag(covs))
    Dmain_a = perr[0]

    if(bBase):
        main_a_corr = main_a - base_a
        Dmain_a_corr = Dmain_a + Dbase_a

    if(autoshow):
        show()
kfsz = fitmaincorrected
fszk = fitmaincorrected
fmc = fitmaincorrected
fcm = fitmaincorrected




def fitexp(*args):
    global bExp
    global exp_A
    global exp_b
    global Dexp_b
    global exp_c
    global exp_start
    global exp_end


    if(not bBase):
        print("Error: missing base fit")
        return


    if(len(args) == 2):
        xmin = args[0]
        xmax = args[1]
    else:
        xmin = float(input("Fit start: "))
        xmax = float(input("Fit end: "))
    
    exp_start = xmin
    exp_end = xmax


    xf = []
    yf = []
    for i in range(len(x)):
        if(xmin <= x[i] <= xmax):
            xf.append(x[i])
            yf.append(y[i])
    a_guess = 3.6
    b_guess = 0.001
    c_guess = 20
    params, covs = curve_fit(lambda t, a, b, c: a * np.exp(-b * t) + c, xf, yf, p0=(a_guess, b_guess, c_guess))

    log("\nExponential fit in x["+str(xmin)+"; "+str(xmax)+"],   y = p0*exp(-p1 * x) + p2")
    fitreport(params, covs)
    bExp = True
    exp_A = params[0]
    exp_b = params[1]
    exp_c = params[2]

    perr = np.sqrt(np.diag(covs))
    Dexp_b = perr[1]

    if(dualfit):
        eszusz()

    if(autoshow):
        show()
fe = fitexp
utoszakasz = fitexp
usz = fitexp


'''
# Amugy is global fit kellene ha ilyet akarunk: https://mail.python.org/pipermail/scipy-user/2013-April/034403.html
def globalfit(*args):
    global bExp
    global exp_A
    global exp_b
    global Dexp_b
    global exp_c

    global y #ez meg minek kell

    # fixme: this doesn't work yet


    if(not bBase):
        print("Error: missing base fit")
        return
    if(len(args) == 2):
        xmin = args[0]
        xmax = args[1]
    else:
        xmin = float(input("Fit start: "))
        xmax = float(input("Fit end: "))


    # create datasets and xf sets
    yf_exp = []
    yf_lin = []
    xf_exp = []
    xf_lin = []
    for i in range(len(x)):
        if(base_start <= x[i] <= base_end):
            yf_lin.append(y[i])
            xf_lin.append(x[i])
        if(xmin <= x[i] <= xmax):
            yf_exp.append(y[i])
            xf_exp.append(x[i])


    data = np.array([xf_exp,yf_exp,xf_lin,yf_lin])

    # get starting value for fit
    Tk, DTk = getTk()

    from lmfit import minimize, Parameters, report_fit

    def exponential(x, A, b, c):
        return A * np.exp(-b * x) + c
    def linconst(x, c):
        return c

    def func_dataset(params, func, x):
        amp = params['amp_%i' % (i+1)].value
        tau = params['tau_%i' % (i+1)].value
        lim = params['lim_%i' % (i+1)].value
        if(func == "exp"):
            return exponential(x, amp, tau, lim)
        return linconst(x,lim)

    def objective(params, data):
        resid1 = 0.0*data[1]
        resid2 = 0.0*data[3]
        # make residual per data set
        resid1[:] = data[1, :] - func_dataset(params, "exp", data[0])
        resid2[:] = data[3, :] - func_dataset(params, "lin", data[2])

        return resid1+resid2



    # create sets of parameters, one per data set
    fit_params = Parameters()
    for iy in range(2):
        fit_params.add( 'amp_%i' % (iy+1), value=3.6)
        fit_params.add( 'tau_%i' % (iy+1), value=0.001)
        fit_params.add( 'lim_%i' % (iy+1), value=Tk)

    # but now constrain all values of lim to have the same value
    #for iy in (2, 3, 4, 5):
    fit_params['lim_%i' % 2].expr='lim_1'

    # run the global fit to all the data sets
    result = minimize(objective, fit_params, args=(data))
    report_fit(result)
'''

def getTs(*args):
    global Ts
    global bTs
    if(not bBase):
        print("Error: base fit not found")
        return
    if(not bExp):
        print("Error: exponential fit not found")
        return
    
    if(len(args) == 1):
        startpoint = args[0]
    else:
        startpoint = base_end

    Ts = []

    
    for i in range(len(x)):
        if(x[i] < startpoint):
            Ts.append(y[i])
        else:
            t = x[i]
            xf = []
            yf = []
            for j in range(len(x)):
                if(startpoint <= x[j] <= t):
                    xf.append(x[j])
                    yf.append(y[j])


            Tk, DTk = getTk(x[j])
            # ez igy garantalja hogy ha van dualfit akkor azt kapjuk, ha nem akkor pedig az exp_limeszt kapjuk vissza

            for j in range(len(xf)):
                yf[j] -= Tk

            Ts.append(y[i] + exp_b*trapcalc(xf,yf))
    bTs = True
    if(autoshow):
        show()
korrigaltT = getTs
kT = getTs


def fitmainexp(*args):
    global bMainExp
    global mainexp_A
    global mainexp_b
    global Dmainexp_b
    global mainexp_c
    global mainexp_start
    global mainexp_end
    global bEpszilonVesszo
    global ev

    if(not bTs):
        print("Main exp must be fit onto T*!")
        return


    if(len(args) == 2):
        xmin = args[0]
        xmax = args[1]
    else:
        xmin = float(input("Fit start: "))
        xmax = float(input("Fit end: "))
    
    mainexp_start = xmin
    mainexp_end = xmax


    xf = []
    yf = []
    for i in range(len(x)):
        if(xmin <= x[i] <= xmax):
            xf.append(x[i])
            yf.append(Ts[i])
    a_guess = -2.44586e+06
    b_guess = 0.0426355
    c_guess = 22.0029
    params, covs = curve_fit(lambda t, a, b, c: a * np.exp(-b * t) + c, xf, yf, p0=(a_guess, b_guess, c_guess))

    log("\nExponential fit (main) in x["+str(xmin)+"; "+str(xmax)+"],   y = p0*exp(-p1 * x) + p2")
    fitreport(params, covs)
    bMainExp = True
    mainexp_A = params[0]
    mainexp_b = params[1]
    mainexp_c = params[2]

    perr = np.sqrt(np.diag(covs))
    Dmainexp_b = perr[1]


    with open(directory+"OUTfajho_beejtepszilon.dat", "w") as f:
        print("# b (ev) [1/s]", file=f)
        print(mainexp_b,file=f)

    bEpszilonVesszo = True
    ev = mainexp_b

    if(autoshow):
        show()
fme = fitmainexp
fem = fitmainexp
expfoszakasz = fitmainexp
foszakaszexp = fitmainexp
efsz = fitmainexp
fsze = fitmainexp





def eszusz(limit="fixed"):
    global bBaseExp, baseexp_A, baseexp_b, baseexp_c

    if(not bBase):
        print("Error eszusz(): missing base fit")
        return
    if(not bExp):
        print("Error eszusz(): missing exp fit")
        return
    
    xf = []
    yf = []
    for i in range(len(x)):
        if(base_start <= x[i] <= base_end):
            xf.append(x[i])
            yf.append(y[i])

        if(i > len(x) - 10):
            xf.append(x[i])
            yf.append(exp_c)


    mp = (exp_c - base_b)/base_a

    a_guess = 3.6
    b_guess = 0.001
    c_guess = exp_c
    #if(mp < base_end):
    #    getav = []
    #    for i in range(len(x)):
    #        if(base_end - 10 <= x[i] <= base_end + 1):
    #            getav.append(y[i])
    #    c_guess = np.average(getav)
    if(limit == "free"):
        params, covs = curve_fit(lambda t, a, b, c: a * np.exp(-b * t) + c, xf, yf, p0=(a_guess, b_guess, c_guess))
        log("\nDual fit: Exponential base correction,   y = p0*exp(-p1 * x) + c")
        print("Dual fit: Exponential base correction,   y = p0*exp(-p1 * x) + c")
    elif(limit == "fixed"):
        params, covs = curve_fit(lambda t, a, b: a * np.exp(-b * t) + exp_c, xf, yf, p0=(a_guess, b_guess))
        log("\nDual fit: Exponential base correction [fixed limit],   y = p0*exp(-p1 * x) + "+str(exp_c))
        print("Dual fit: Exponential base correction [fixed limit],   y = p0*exp(-p1 * x) + "+str(exp_c))
    else:
        print("bad argument to 'limit' (free/fixed)")
        return
    fitreport(params, covs)
    bBaseExp = True
    baseexp_A = params[0]
    baseexp_b = params[1]
    baseexp_c = params[2] if limit == "free" else exp_c


    if(autoshow):
        show()
df = eszusz

def calib(save=True):
    global bCalib
    global Cp
    global dCp
    global alfa
    global dalfa

    if(not bBase):
        print("Error: missing base")
        return
    if(not bMain):
        print("Error: missing main")
        return
    if(not bExp):
        print("Error: missing exp")
        return

    print("R = ("+str(R)+" +/- "+str(DR)+") Ohm")


    # get U
    Uheat = []
    for i in range(len(x)):
        if(main_start <= x[i] <= main_end):
            Uheat.append(Uarr[i])
    U = np.average(Uheat)
    DU = np.std(Uheat)
    print("U = ("+str(U)+" +/- "+str(DU)+") V")

    print("a_corr = ("+str(main_a_corr)+" +/- "+str(Dmain_a_corr)+") K/s")

    print("Beta = ("+str(exp_b)+" +/- "+str(Dexp_b)+") 1/s")
    # get Cp
    Cp = U**2/R/main_a_corr
    dCp = 2*DU/U + DR/R + Dmain_a_corr/main_a_corr
    print("=======================")
    print("Cp = ("+str(Cp)+" +/- "+str(dCp*Cp)+") J/K")

    alfa = exp_b * Cp
    dalfa = dCp + Dexp_b/exp_b
    print("Alfa = ("+str(alfa)+" +/- "+str(dalfa*alfa)+") W/K")
    bCalib = True

    log("\nCalibration data:")
    log("R = ("+str(R)+" +/- "+str(DR)+") Ohm")
    log("U = ("+str(U)+" +/- "+str(DU)+") V")
    log("a_corr = ("+str(main_a_corr)+" +/- "+str(Dmain_a_corr)+") K/s")
    log("Beta = ("+str(exp_b)+" +/- "+str(Dexp_b)+") 1/s")
    log("=======================")
    log("Cp = ("+str(Cp)+" +/- "+str(dCp*Cp)+") J/K")
    log("Alfa = ("+str(alfa)+" +/- "+str(dalfa*alfa)+") W/K\n")


    if(save):
        with open(directory+"OUTfajho_kalibracio.dat", "w") as f:
            print("# Cp [J/K], D(Cp) [J/K], alfa [W/K], D(alfa) [W/K]", file=f)
            print(Cp, dCp*Cp, alfa, dalfa*alfa,file=f)
c = calib
k = calib
kalib = calib
kalibracio = calib
calibrate = calib



def trapcalc(x,y):
    sum = 0
    for i in range(len(x)-1):
        if(y[i] < y[i+1]):
            sum += y[i] * (x[i+1]-x[i])  +  (y[i+1]-y[i]) * (x[i+1]-x[i])/2
        else:
            sum += y[i+1] * (x[i+1]-x[i])  +  (y[i]-y[i+1]) * (x[i+1]-x[i])/2
    return sum



def integrate(*args):
    if(not bCalib):
        print("Error: missing calib")
        return
    if(len(args) == 0):
        xmin = float(input("Integral start: "))
        xmax = float(input("Integral end: "))
    elif(len(args) == 2):
        xmin = args[0]
        xmax = args[1]
    else:
        print("Insufficient arguments")
        return

    xf = []
    yf = []
    for i in range(len(x)):
        if(xmin <= x[i] <= xmax):
            xf.append(x[i])
            yf.append(y[i])


    Tt = yf[-1]

    if(not bBase):
        fitbase(0,xmin)


    Tm, DTm = getTm(xmin)

    Tk, DTk = getTk()

    for i in range(len(xf)):
        if(bBaseExp):
            sway = fexp(xf[i], baseexp_A, baseexp_b, baseexp_c)
        else:
            sway = flin(xf[i], base_a, base_b)
        yf[i] -= sway

    integ = trapcalc(xf,yf)

    cm = (Cp*(Tt - Tk)+alfa*integ)/( Tm-Tt )

    # assuming DTk deviation for T(t) as well, and assuming DI ~ 0
    Dcm = abs((Tt - Tk)/(Tm-Tt))* (dCp*Cp) + abs((alfa*integ + Cp*(Tm - Tk)/(Tm - Tt)**2))*DTk + abs(-Cp/(Tm - Tt))*DTk + abs(integ/(Tm - Tt))*(dalfa*alfa) + abs((Cp*(Tk - Tt) - alfa*integ)/(Tm - Tt)**2)* DTm

    print("Tm = ("+str(Tm)+" +/- "+str(DTm)+") C")
    print("T0 = ("+str(Tk)+" +/- "+str(DTk)+") C")
    print("T(t) = "+str(Tt)+" C")
    print("Integral = "+str(integ)+" Cs")
    print("=======================")
    print("Cm = ("+str(cm)+" +/- "+str(Dcm)+") J/K")

    with open(directory+"OUTfajho_beejtesEgyszeruIntegral.dat", "w") as f:
        print("Tm = ("+str(Tm)+" +/- "+str(DTm)+") C", file=f)
        print("T0 = ("+str(Tk)+" +/- "+str(DTk)+") C", file=f)
        print("T(t) = "+str(Tt)+" C", file=f)
        print("Integral = "+str(integ)+" Cs", file=f)
        print("=======================", file=f)
        print("Cm = ("+str(cm)+" +/- "+str(Dcm)+") J/K", file=f)
integral = integrate





def getTk(*args):
    # ha kell Tk, ehhez kell at least base

    # fv viselkedese:
    # DTk szorast eloszakaszbol veszi
    # ha van beadva x pont neki, akkor ott veszi Tk-t, ha nincs akkor a base_end kornyeken
    # ha van dualfit, akkor azt hasznalja mindenhol.
    # ha nincs dualfit, de van exponencialis, akkor
    #   ha van beadott x-pont, akkor az exp limeszt adja vissza
    #   ha nincs (azaz base_end) nezzuk, akkor ott exp limeszt es linearis base fitet atlagol
    # ha nincs exp se, akkor a linearis base futet hasznalja mindenhol

    if(not bBase):
        print("Error in getTk(): missing base fit")
        return False

    aTk = []
    for i in range(len(x)):
        if(base_start <= x[i] <= base_end):
            aTk.append(y[i])
    DTk = np.std(aTk)


    if(len(args) == 1):
        xp = args[0]
        if(bBaseExp):
            return fexp(xp,baseexp_A,baseexp_b,baseexp_c), DTk
        if(bExp):
            return exp_c, DTk
        return flin(xp,base_a,base_b), DTk


    if(bBaseExp):
        aTk = []
        for i in range(len(x)):
            if(base_end - 10 <= x[i] <= base_end + 10):
                aTk.append( fexp(x[i], baseexp_A, baseexp_b, baseexp_c) )
        return np.average(aTk), DTk
    if(bExp):
        aTk = []
        for i in range(len(x)):
            if(base_end - 10 <= x[i] <= base_end + 10):
                aTk.append( flin(x[i], base_a, base_b) )
        return (np.average(aTk) + exp_c)/2.0, DTk
    
    aTk = []
    for i in range(len(x)):
        if(base_end - 10 <= x[i] <= base_end + 10):
            aTk.append( flin(x[i], base_a, base_b) )
    return np.average(aTk), DTk

    '''
    if(len(args) == 0):
        if(Tkmode == "dualfit"):
            if(bBaseExp):
                return baseexp_c, DTk
            return Tk, DTk
        elif(Tkmode == "base"):
            return Tk, DTk
        else: # "exp"
            return exp_c, DTk
    elif(len(args) == 1):
        xp = args[0]
        if(Tkmode == "dualfit"):
            if(bBaseExp):
                return fexp(xp,baseexp_A,baseexp_b,baseexp_c), DTk
            return flin(xp,base_a,base_b), DTk
        elif(Tkmode == "base"):
            return flin(xp,base_a,base_b), DTk
        else: # "exp"
            return exp_c, DTk

    # somethings wrong
    return False
    '''


def getTm(position):
    Tm_l = []
    for i in range(len(x)):
        if(position - 10 <= x[i] <= position + 5):
            Tm_l.append(heater[ i ])
    return np.average(Tm_l), np.std(Tm_l)



def serialintegrate(*args):
    if(not bCalib):
        print("Error: missing calib")
        return
    if(not bBase):
        fitbase(0,xmin)

    if(len(args) == 3):
        xmin = args[0]
        xmax1 = args[1]
        xmax2 = args[2]
    else:
        xmin = float(input("Integral start: "))
        xmax1 = float(input("Integral end low: "))
        xmax2 = float(input("Integral end high: "))


    Tm, DTm = getTm(xmin)

    # kell Tk
    Tk, DTk = getTk()

    integrals = []
    uplims = []
    Tts = []

    for i in range(len(x)):
        if(xmax1 <= x[i] <= xmax2):
            uplims.append(x[i])
            xf = []
            yf = []
            for j in range(len(x)):
                if(xmin <= x[j] <= x[i]):
                    xf.append(x[j])
                    yf.append(y[j])

            # base correction
            for j in range(len(xf)):
                if(bBaseExp):
                    sway = fexp(xf[j], baseexp_A, baseexp_b, baseexp_c)
                else:
                    sway = flin(xf[j], base_a, base_b)
                yf[j] -=sway

            integ = trapcalc(xf,yf)
            integrals.append(integ)
            Tts.append(y[i])  # yes this is i

    cms = []
    for i in range(len(integrals)):
        cms.append( (Cp*(Tts[i] - Tk)+alfa*integrals[i])/( Tm-Tts[i] ) )


    plt.close()
    plt.plot(uplims,cms)
    plt.xlabel("Integral upper limit (s)")
    plt.ylabel(r'$C_{m}$ $(J/K)$')

    with open(directory+"OUTfajho_beejtesSorozatosIntegral.dat", "w") as f:
        print("# t [s], Cm [J/K]", file=f)
        for i in range(len(uplims)):
            print(uplims[i], cms[i], file=f)

ismeteltintegral = serialintegrate
repeatintegrate = serialintegrate
si = serialintegrate
ii = serialintegrate



def beejtes(*args):
    if(not bCalib):
        print("Error: missing calibration")
        return
    if(not bBase):
        print("Error: missing base fit")
        return
    if(not bMainExp):
        print("Error: missing main fit")
        return
    if(not bExp):
        print("Error: missing exp fit")
        return
    if(not bTs):
        print("Error: missing T*")
        return

    if(len(args) == 1):
        Tsint_start = args[0]
    else:
        Tsint_start = exp_start
    Tsav, DTsav = getAvTs(Tsint_start)


    Tk, DTk = getTk()
    Tm, DTm = getTm(mainexp_start)

    K = getK()
    Tms = Tk + K*(Tsav - Tk)
    # deltaK = e0/ev*(delta e0 + delta ev) = e0/ev*6%
    deltaK = (alfa/Cp) / ev * 0.06 
    DTms = K*DTsav + (K - 1)*DTk + (Tsav - Tk)*deltaK
    cm = Cp*(Tsav- Tk)/(Tm - Tms)

    dcm = dCp + (DTsav + DTk)/(Tsav - Tk) + (DTm + DTsav)/(Tm - Tsav)

    print("Tm = ("+str(Tm)+" +/- "+str(DTm)+") C")
    print("Tm* = ("+str(Tms)+" +/- "+str(DTms)+") C")
    print("Tk = ("+str(Tk)+" +/- "+str(DTk)+") C")
    print("T* = ("+str(Tsav)+" +/- "+str(DTsav)+") C")
    print("=======================")
    print("Cm = ("+str(cm)+" +/- "+str(dcm * cm)+") J/K")

    with open(directory+"OUTfajho_beejtesRegi.dat", "w") as f:
        print("Tm = ("+str(Tm)+" +/- "+str(DTm)+") C", file=f)
        print("Tm* = ("+str(Tms)+" +/- "+str(DTms)+") C", file=f)
        print("Tk = ("+str(Tk)+" +/- "+str(DTk)+") C", file=f)
        print("T* = ("+str(Tsav)+" +/- "+str(DTsav)+") C", file=f)
        print("=======================")
        print("Cm = ("+str(cm)+" +/- "+str(dcm * cm)+") J/K", file=f)
beejt = beejtes


def rafut():
    if(not bCalib):
        print("Error: missing calib")
        return
    if(not bBase):
        print("Error: missing base fit")
        return
    if(not bMain):
        print("Error: missing main fit")
        return

    print("R = ("+str(R)+" +/- "+str(DR)+") Ohm")


    # get U
    t, U, DU = getheatingdata()
    print("\nCalculating Cm data:")
    print("R = ("+str(R)+" +/- "+str(DR)+") Ohm")
    print("U = ("+str(U)+" +/- "+str(DU)+") V")
    print("a_corr = ("+str(main_a_corr)+" +/- "+str(Dmain_a_corr)+") K/s")
    # get Ck
    Ck = U**2/R/main_a_corr
    dCk = 2*DU/U + DR/R + Dmain_a_corr/main_a_corr
    print("=======================")
    print("Ck = ("+str(Ck)+" +/- "+str(dCk*Ck)+") J/K")
    print("Cm = ("+str(Ck-Cp)+" +/- "+str(dCk*Ck + dCp*Cp)+") J/K")

    log("\nCalculating Cm data:")
    log("R = ("+str(R)+" +/- "+str(DR)+") Ohm")
    log("U = ("+str(U)+" +/- "+str(DU)+") V")
    log("a_corr = ("+str(main_a_corr)+" +/- "+str(Dmain_a_corr)+") K/s")
    log("=======================")
    log("Ck = ("+str(Ck)+" +/- "+str(dCk*Ck)+") J/K")
    log("Cm = ("+str(Ck-Cp)+" +/- "+str(dCk*Ck + dCp*Cp)+") J/K")


    with open(directory+"OUTfajho_rafutes.dat", "w") as f:
        print("# Ck [J/K], D(Ck) [J/K], Cm [J/K], D(Cm) [J/K]", file=f)
        print(Ck, dCk*Ck, Ck-Cp, dCk*Ck + dCp*Cp,file=f)
r = rafut
rafutes = rafut
rafuteses = rafut



def getheatingdata():
    act = False
    start = 0
    end = 0
    Uact = []
    for i in range(len(x)):
        if Uarr[i] < 0.5 and act:
            end = x[i]
            break

        if Uarr[i] > 0.5 and not act:
            start = x[i]
            act = True

        if act:
            Uact.append(Uarr[i])

    # turnon workaround
    counts = 0
    av = np.average(Uact)
    retake = True
    while(retake):
        for i in range(len(Uact)):
            if(Uact[i] < av - 0.1 or Uact[i] > av + 0.1):
                counts += 1
                del Uact[i]
                break
        else:
            retake = False


    return end - start, np.average(Uact), np.std(Uact)


def rafutteljes(*args):
    if(not bCalib):
        print("Error: missing calib")
        return
    if(not bEpszilonVesszo):
        print("Error: epsilon' not found")
        return
    if(not bBase):
        print("Error: base fit not found")
        return
    if(not bTs):
        print("Error: T* fit not found")
        return
    
    if(len(args) == 1):
        Tsint_start = args[0]
    else:
        Tsint_start = exp_start


    # cm = z/K
    #K = ev/(ev-e0),  e0: b kalib, ev: b beejtes
    # z = Q(Ts-Tk) - v

    t, U, DU = getheatingdata()
    Q = U**2/R *t
    dQ = 2*DU/U + DR/R
    Tk, DTk = getTk()

    Tsav, DTsav = getAvTs(Tsint_start)

    z = Q/(Tsav - Tk) - Cp

    Dtort = Q/(Tsav - Tk) * (dQ + (DTsav + DTk)/(Tsav - Tk))
    Dz = Dtort + dCp*Cp

    dK = 0.06 * (alfa/Cp)/ev 
    cm = z/getK()
    dcm = Dz/z + dK
    

    print("\nCalculating Cm data:")
    print("R = ("+str(R)+" +/- "+str(DR)+") Ohm")
    print("U = ("+str(U)+" +/- "+str(DU)+") V")
    print("t = "+str(t)+" s")
    print("Q = ("+str(Q)+" +/- "+str(dQ*Q)+") J")
    print("Tk = ("+str(Tk)+" +/- "+str(DTk)+") C")
    print("T* = ("+str(Tsav)+" +/- "+str(DTsav)+") C")
    print("=======================")
    print("Cm = ("+str(cm)+" +/- "+str(dcm*cm)+") J/K")
    print("Cm [K=1] = ("+str(z)+" +/- "+str(Dz)+") J/K")


    log("\nCalculating Cm data:")
    log("R = ("+str(R)+" +/- "+str(DR)+") Ohm")
    log("U = ("+str(U)+" +/- "+str(DU)+") V")
    log("t = "+str(t)+" s")
    log("Q = ("+str(Q)+" +/- "+str(dQ*Q)+") J")
    log("Tk = ("+str(Tk)+" +/- "+str(DTk)+") C")
    log("T* = ("+str(Tsav)+" +/- "+str(DTsav)+") C")
    log("=======================")
    log("Cm = ("+str(cm)+" +/- "+str(dcm*cm)+") J/K")
    log("Cm [K=1] = ("+str(z)+" +/- "+str(Dz)+") J/K")



def getAvTs(start):
    Tsavarr = []
    for i in range(len(x)):
        if( start < x[i]):
            Tsavarr.append(Ts[i])
    return np.average(Tsavarr), np.std(Tsavarr)




def process_data(data_in):
    t = []
    T = []
    heater = []
    Uarr = []
    for i in range(len(data_in)):
        t.append(data_in[i][0])
        T.append(data_in[i][1])
        heater.append(data_in[i][3])
        Uarr.append(data_in[i][4])
    return t,T,heater,Uarr



colors = {}
colors["data"] = "blue"
colors["base"] = "dimgrey"
colors["main"] = "darkorange"
colors["exp"] = "crimson"
colors["T*"] = "green"



def show():
    global fig
    if(not fig):
        fig = plt.figure()
    plt.clf()
    plt.plot(x,y, color=colors["data"], label='Kaloriméter hőmérséklet')
    plt.xlabel("t (s)")
    plt.ylabel(r'T $(^{\circ}C)$')


    auto_ylim = plt.gca().get_ylim()
    auto_xlim = plt.gca().get_xlim()


    if(grid):
        plt.grid("true",linestyle="--")

    xs = np.linspace(x[0], x[-1], 100)
    if(bBase):
        #plt.autoscale(False)
        if(bBaseExp):
            plt.plot(xs,fexp(xs,baseexp_A,baseexp_b,baseexp_c),color=colors["base"], label="Előszakasz illesztés (exp. korr.)")
        else:
            plt.plot(xs,flin(xs,base_a,base_b),color=colors["base"], label="Előszakasz illesztés")

    if(bTs):
        #plt.autoscale(True)
        plt.plot(x,Ts, color=colors["T*"], label="Korrigált hőmérséklet")
        auto_ylim = plt.gca().get_ylim()

    if(bMain):
        #plt.autoscale(False)
        plt.plot(xs,flin(xs,main_a,main_b), color=colors["main"], label=("Főszakasz illesztés (korr.)" if bMainCorrected else "Főszakasz illesztés") )
    if(bMainExp):
        xs2 = np.linspace(mainexp_start, mainexp_end, 100)
        #plt.autoscale(False)
        plt.plot(xs2,fexp(xs2,mainexp_A,mainexp_b,mainexp_c), color=colors["main"], label="Exp. főszakasz illesztés")
    if(bExp):
        #plt.autoscale(False)
        plt.plot(xs,fexp(xs,exp_A,exp_b,exp_c), color=colors["exp"], label="Utószakasz illesztés")
        plt.plot(xs,flin(xs,0,exp_c), color=colors["exp"], linestyle="--")

    if(bLegend):
        plt.legend()

    if(scaletodata):
        plt.xlim(auto_xlim)
        plt.ylim(auto_ylim)

s = show
p = show


def log(text):
    with open(directory+"OUTfajho_fitlog.dat", "a") as f:
        print(text, file=f)


def loadcalib():
    global bCalib
    global Cp
    global dCp
    global alfa
    global dalfa
    path = easygui.fileopenbox()
    if(path):
        calin = np.loadtxt(path,skiprows = 1)
        Cp = calin[0]
        dCp = calin[1] / Cp
        alfa = calin[2]
        dalfa = calin[3] / alfa
        bCalib = True
        print("Calibration loaded:")
        print("Cp = ("+str(Cp)+" +/- "+str(dCp*Cp)+") J/K")
        print("Alfa = ("+str(alfa)+" +/- "+str(dalfa*alfa)+") W/K")
    else:
        print("Error: no calibration file found in path")





def cd(*args):
    global directory
    if(len(args)==0):
        path = easygui.diropenbox()
    else:
        path = args[0]

    try:
        os.chdir(path)
        print("Current working directory: {0}".format(os.getcwd()))
        directory = path
    except FileNotFoundError:
        print("Directory: {0} does not exist".format(path))
    except NotADirectoryError:
        print("{0} is not a directory".format(path))
    except PermissionError:
        print("You do not have permissions to change to {0}".format(path))




def png():
    plt.rcParams["savefig.format"] = "png"
def pdf():
    plt.rcParams["savefig.format"] = "pdf"






def new(*args):
    global directory

    global bBase
    global bMain
    global bExp
    global bInt
    global bTs
    global bMainExp
    global bBaseExp
    global bMainCorrected

    global x
    global y
    global heater
    global Uarr
    global R

    global bCalib
    global Cp
    global dCp
    global alfa
    global dalfa

    global bEpszilonVesszo
    global ev

    if(len(args)==0):
        path = easygui.fileopenbox(filetypes=[["*","All files"],["*.*","All files"]]) 
    else:
        path = args[0]
    if(path):

        # get the directory
        pathtrunc = path
        while(True):
            if(pathtrunc[-1] == '/'):
                directory = pathtrunc
                break
            else:
                pathtrunc = pathtrunc[:-1]
            
            if(len(pathtrunc) <= 0):
                print("Error while recursing for directory")
                return

        # workaround: first log message doesn't seem to go out, if log does not exist, create it
        if(not os.path.exists(directory+"OUTfajho_fitlog.dat")):
            with open(directory+"OUTfajho_fitlog.dat", "a") as f:
                print("", file=f)
            
        bBase = False
        bMain = False
        bMainCorrected = False
        bExp = False
        bInt = False
        bTs = False
        bMainExp = False
        bCalib = False
        bEpszilonVesszo = False
        bBaseExp = False
        data0 = np.loadtxt(path,skiprows=2)

        with open(path) as f:
            R_line = f.readline()
            R = float(re.sub("[^\d\.]", "", R_line))

        t = time.localtime()
        current_time = time.strftime("%H:%M:%S", t)
        log("\nOpened file: "+path+" at ["+str(current_time)+"]\n")

        try:
            os.chdir(directory)
            #print("Current working directory: {0}".format(os.getcwd()))
        except FileNotFoundError:
            print("Directory: {0} does not exist".format(directory))
        except NotADirectoryError:
            print("{0} is not a directory".format(directory))
        except PermissionError:
            print("You do not have permissions to change to {0}".format(directory))

        if(os.path.exists(directory+"OUTfajho_kalibracio.dat")):
            calin = np.loadtxt(directory+"OUTfajho_kalibracio.dat",skiprows = 1)
            Cp = calin[0]
            dCp = calin[1] / Cp
            alfa = calin[2]
            dalfa = calin[3] / alfa
            bCalib = True
            print("Existing calibration found:")
            print("Cp = ("+str(Cp)+" +/- "+str(dCp*Cp)+") J/K")
            print("Alfa = ("+str(alfa)+" +/- "+str(dalfa*alfa)+") W/K")

        if(os.path.exists(directory+"OUTfajho_beejtepszilon.dat")):
            evin = np.loadtxt(directory+"OUTfajho_beejtepszilon.dat",skiprows = 1)
            ev = evin
            bEpszilonVesszo = True
            print("Existing epszilon' found:")
            print("ev = "+str(ev)+" 1/s")


        x,y,heater,Uarr = process_data(data0)
        s()
    else:
        print("Error: no file found in path")
o = new
uj = new


new()

import matplotlib as mpl 
mpl.rcParams["savefig.directory"] = ""
plt.rcParams["savefig.format"] = "pdf"

code.interact(local=locals())