import numpy as np
import math
import sympy as sp
from scipy.integrate import quad
import tkinter as tk
from tkinter import filedialog

root = tk.Tk()
root.withdraw()

file_path = filedialog.askopenfilename()

msmt= np.genfromtxt(file_path,
                    dtype = float,
                    invalid_raise = False,
                    usemask = False,
                    delimiter = " ",
                    usecols = (0,1,2,3,4,5,6,7,8,9,10),
                    skip_header = 3,
                    missing_values = 0.0)

F1_64 = 1/64
F63_64 = 63/64
F63_16 = 63/16
F1_64_20 = 1/64/20


import matplotlib.pyplot as plt
#1st col
t = msmt[:,0]
len_data = len(t)

#2nd col
el = msmt[:,1]/1000

#3rd col
et = msmt[:,2]/1000

#4th col
fl = msmt[:,3]/1000

#5th col
ft = msmt[:,4]/1000

#6th col
sf = msmt[:,5]

#7th col
eu = msmt[:,6]
for i in range(len_data):
    eu[i]=(eu[i]-eu[0])*10

#8th col
t1  = msmt[:,7]

#9th col
t2  = msmt[:,8]

#10th col
t3  = msmt[:,9]

#11th col
t4  = msmt[:,10]

d = 9.3e-3
pi = math.pi
a0=pi / 4 * d**2

ed=np.zeros(len_data)
fd=np.zeros(len_data)
es=np.zeros(len_data)
fs=np.zeros(len_data)
for i in range(len_data):
    ed[i] = (el[i] - et[i]) * 2 / 3
    fd[i] = (fl[i] - ft[i]) * 2 / 3
    es[i] = (el[i] + 2 * et[i]) / 3
    fs[i] = (fl[i] + 2 * ft[i]) / 3

#color vector to plot dots (prepared for cycles)
colors_dot = ['b.', 'g.','r.','c.','m.','y.','k.']
#color vector to plot big dots (prepared for cycles)
colors_o = ['bo', 'go','ro','co','mo','yo','ko']
#color vector to plot lines (prepared for cycles)
colors_line = ['b', 'g','r','c','m','y','k']





import matplotlib as mpl
mpl.rc('font',family='Palatino Linotype')
start = 29
end = 3097
o = range(len_data)

# figure 0 begin
oo = o[start:end]
ar_p = .005*sf[start-1:end-1]
al_p = el[start-1:end-1]
at_p = et[start-1:end-1]
#ap_p = -(et[start-1:end-1]-et[start-1]) / (el[start-1:end-1]-el[start-1])
#ad_p = ed[start-1:end-1]
#as_p = es[start-1:end-1]
plt.figure(0, figsize=(10,5))

plt.plot(oo,ar_p,colors_line[6],label='ar')
plt.plot(oo,al_p,colors_line[2],label='al')
plt.plot(oo,at_p,colors_line[1],label='at')

plt.title('', fontsize = 15)
plt.ylabel('')
plt.xlabel('Measurement point number')
plt.axis([start, end, -30, 50])
plt.xticks(size = 10, rotation = 30)
plt.yticks(size = 10)
plt.grid(True)
plt.legend(loc=1)
plt.savefig("keplekenyedeshez.jpg", dpi=300)

# figure 0 end

#gömbi illesztés
ee=np.zeros(len_data)
ss=np.zeros(len_data)
for i in range(len_data):
    ee[i] = es[i]
    ss[i] = 1e-6/a0 * 1 / 3 * sf[i] / (1 + (et[i] - et[0]) / 1000)**2



u = sp.Symbol('u')
u2 = u*u
u3 = u2*u
u4 = u3*u
u5 = u4*u
u6 = u5*u
u7 = u6*u
u8 = u7*u
u9 = u8*u
u10 = u9*u
P0 = F1_64 * (63*u - 42*u3 + 18*u7 - 7*u9)
P1 = F1_64_20 * (1260*u + 630*u2 - 840*u3 - 630*u4 + 360*u7 + 315*u8 - 140*u9 - 126*u10)
DP0 = F63_64 * (-2*u2 + 2*u6 - u8)
DP1 = F1_64 * (-126*u2 - 84*u3 + 126*u6 + 108*u7 - 63*u8 - 56*u9)
D2P0 = - F63_16 * u * (2*u2 + 1) * (u2 - 1)**2
D2P1 = F63_64 * (-4*u - 2*u2 + 12*u5 + 10*u6 - 8*u7 - 7*u8)

#procedure

#input parameters
jint_start = 29
jint_end = 58
J = [1,2,3,4,5] #[tau_sph, E0_sph, E1_sph, E2_sph, delta_sph]
DNi = 2

#global parameters
PP = np.zeros(5)
QQ = np.zeros(5)
sj = np.zeros(len_data)
ej = np.zeros(len_data)

#local parameters
j_num = jint_end - jint_start
jmin = jint_start - 5
jmax = jmin + j_num
oo = range(len_data)[jmin:jmax+1]
kappa = 1 / 3
# global: Ip, PP,QQ sj, ej, P0, DP0 D2P0, P1, DP1, D2P1, A1, A2, A3, A4, A5, B1, B2, B3, B4, B5
# local: jmin, jmax, oo, shift, Ni, Np, intl, i, j1, j2, t1, t2, Dt, FDt_2, F2_Dt, F2_Dt2, n, ta,
#        tb, ua, ub, dP0, dDP0, dD2P0, dP1, dDP1, dD2P1, fa, fb, c1, c00, s0, s1, s2, s3, s4, L, 
#        j, k iA, P, S, l, s, sss, eee, st, b0, b2, b2, b3, kappa


#progress
sj[jmin] = ss[jmin]
ej[jmin] = ee[jmin]
ej[jmin] = ee[jmin]
sss = 0
eee = 0
L = 0
if J[L]==1:
    L=L+1
if J[L]==2:
    L=L+1
if J[L]==3:
    L=L+1
if J[L]==4:
    L=L+1
if J[L]==5:
    L=L+1
Ni = L + DNi
shift = math.trunc((j_num+1)/Ni)
Np = (jmax-jmin) - (Ni-1) * shift
intl = np.zeros((Ni,2))
for i in range(Ni):
    intl[i] = [jmin + i*shift, jmin + i*shift + Np]
Ip = np.zeros((6,Ni))
for i in range(Ni):
    j1 = intl[i,0]
    j2 = intl[i,1]
    t1 = t[math.trunc(j1)]
    t2 = t[math.trunc(j2)]
    Dt = t2 - t1
    FDt_2 = Dt/2
    F2_Dt = 2/Dt
    F2_Dt2 = F2_Dt**2
    s0 = 0
    s1 = 0
    s2 = 0
    s3 = 0
    s4 = 0
    for k in range(math.trunc(j1),math.trunc(j2)):
        ta = t[k]
        tb = t[k + 1]
        ua = -1 + F2_Dt*(ta - t1)
        ub = -1 + F2_Dt*(tb - t1)

        P0_ub = P0.subs(u,ub).doit()
        P0_ua = P0.subs(u,ua).doit()
        dP0 = P0_ub - P0_ua
        
        DP0_ub = DP0.subs(u,ub).doit()
        DP0_ua = DP0.subs(u,ua).doit()
        dDP0 = -F2_Dt*(DP0_ub - DP0_ua)

        D2P0_ub = D2P0.subs(u,ub).doit()
        D2P0_ua = D2P0.subs(u,ua).doit()
        dD2P0 = F2_Dt2*(D2P0_ub - D2P0_ua)

        P1_ub = P1.subs(u,ub).doit()
        P1_ua = P1.subs(u,ua).doit()
        dP1 = FDt_2*(P1_ub - P1_ua)

        DP1_ub = DP1.subs(u,ub).doit()
        DP1_ua = DP1.subs(u,ua).doit()
        dDP1 = -(DP1_ub - DP1_ua)

        D2P1_ub = D2P1.subs(u,ub).doit()
        D2P1_ua = D2P1.subs(u,ua).doit()
        dD2P1 = F2_Dt*(D2P1_ub - D2P1_ua)

        fa = ss[k]
        fb = ss[k+1]
        c1 = (fb - fa) / (tb - ta)
        c00 = c1 * t1 + (tb * fa - ta * fb) / (tb - ta)
        s0 = s0 + c00 * dP0 + c1 * dP1
        s1 = s1 - c00 * dDP0 - c1 * dDP1
        fa = ee[k]
        fb = ee[k+1]
        c1 = (fb - fa) / (tb - ta)
        c00 = c1 * t1 + (tb * fa - ta * fb) / (tb - ta)
        s2 = s2 + c00 * dP0 + c1 * dP1
        s3 = s3 + c00 * dDP0 + c1 * dDP1
        s4 = s4 + c00 * dD2P0 + c1 * dD2P1
    Ip[0,i] = s0
    Ip[1,i] = s1
    Ip[2,i] = s2
    Ip[3,i] = s3
    Ip[4,i] = s4
    Ip[5,i] = 1

A = np.zeros((L,L))
B = np.zeros(L)
for k in range(L):
    for j in range(L):
        for jj in range(Ni):
            A[k,j] += Ip[J[k],jj] * Ip[J[j],jj]
    for jj in range(Ni):
        B[k] += Ip[0,jj] * Ip[J[k],jj]
print(A)
print(B)
iA = np.linalg.inv(A)
P = iA @ B
S = 0
for i in range(Ni):
    add_aux = 0
    for k in range(L):
        add_aux += Ip[J[k],i] * P[k]
    S += (Ip[0,i] - add_aux)**2
s = math.sqrt(S / Ni)
PP = np.zeros(5)
QQ = np.zeros(5)
for k in range(L):
    PP[J[k]-1] = P[k]
    QQ[J[k]-1] = s * np.sqrt(iA[k,k])
sj[jmin] = ss[jmin]
sss = 0
print("\n")
k = 0
for j in range(jmin + 1, jmax + 1):
    k += 1
    sj[j] = (PP[1]/2*(ee[j] + ee[j-1]) + PP[2]*(ee[j] - ee[j-1])/1/(t[j] - t[j-1]) + PP[3]/2*( ( (ee[j+1] - ee[j])/(t[j+1] - t[j]) - (ee[j] - ee[j-1])/(t[j] - t[j-1]) )/( (t[j+1] + t[j])/2 - (t[j] + t[j-1])/2 ) + ( (ee[j] - ee[j-1])/(t[j] - t[j-1]) - (ee[j-1] - ee[j-2])/(t[j-1] - t[j-2]) )/( (t[j] + t[j-1])/2 - (t[j-1] + t[j-2])/2 ) ) + PP[4] + ( PP[0]/(t[j] - t[j-1]) - 1/2)*sj[j-1] ) / ( 1/2 + PP[0]/(t[j] - t[j-1]) )
    sss = sss + math.fabs(sj[j] - ss[j])
    print(str(t[j]) + " " + str(ee[j]) + " " + str(ss[j]) + " " + str(sj[j]))
sss = sss / ((jmax - 1) - (jmin + 1) + 1)

ej[jmin] = ee[jmin]
ej[jmin + 1] = ee[jmin + 1]
eee = 0
kappa = 1/3
for j in range(jmin + 2, jmax + 1):
    ej[j] = ( ( 2*PP[3]/(t[j] - t[j-2]) * ( 1/(t[j] - t[j-1]) + 1/(t[j-1] - t[j-2]) ) - PP[1]*kappa ) * ej[j - 1] + ( PP[2]/(t[j] - t[j-2]) - PP[3]*2/(t[j-1] - t[j-2])/(t[j] - t[j-2]) - PP[1]*(1 - kappa)/2 ) * ej[j-2] + kappa*ss[j-1] + (1 - kappa)/2 * ( ss[j] + ss[j-2] ) + 1/(t[j] - t[j-2]) * ( PP[0] * (ss[j] - ss[j-2]) ) - PP[4] ) / ( PP[1]*(1 - kappa)/2 + PP[2]/(t[j] - t[j-2]) + PP[3]*2/(t[j-1] - t[j-2])/(t[j] - t[j-2]) )

print("\n                      T               U               V               W               X             hiba")
print("Jóslat:           " + str(round(PP[0],6)) + "        " + str(round(PP[1],6)) + "        " + str(round(PP[2],6)) + "        " + str(round(PP[3],6)) + "        " + str(round(PP[4],6)))
print("Absz. hibaja:     " + str(round(QQ[0],6)) + "        " + str(round(QQ[1],6)) + "        " + str(round(QQ[2],6)) + "        " + str(round(QQ[3],6)) + "        " + str(round(QQ[4],6)) + "       " + str(round(sss,6)))


# figure 2 begin

plt.figure(2, figsize=(10,5))
plt.plot(oo,ss[jmin:jmax + 1],colors_dot[3],label='ss')
plt.plot(range(jmin,jmax+1),sj[jmin:jmax+1],colors_line[3],label='sj')
plt.plot(oo,2.8 * ee[jmin:jmax+1],colors_o[4],label='ee')
plt.plot(range(jmin,jmax+1),2.8 * ej[jmin:jmax+1],colors_line[4],label='ss')
plt.title('', fontsize = 15)
plt.ylabel('')
plt.xlabel('Measurement point number')
my_xticks = range(jmin, jmax + 1)
my_xticks_label = ["" for x in range(j_num + 1)]
for i in range(j_num + 1):
    my_xticks_label[i] = str(my_xticks[i] + 5)
plt.axis([jmin + 1, jmax - 1, 0, 7])
plt.xticks(my_xticks, my_xticks_label, size = 10, rotation = 30)
plt.yticks(size = 10)
plt.grid(True)
plt.legend(loc=0)
plt.savefig("joslat.jpg", dpi=300)
#plt.savefig("diagram_1_01.jpg", dpi=300)
# figure 2 end



plt.show()