def mean(values, axis=0):
    """Returns mean values and their mean errors of a given array. Return value will be a unp.uarray
    Args:
            values:     (list)  Array containing numbers whose mean is desired.
            axis:       (int)   Axis along which the means are computed. The default is to compute the mean of the flattened array.
    """
    return unp.uarray((np.mean(noms(values), axis=axis), scipy.stats.sem(noms(values), axis=axis)))

def table(name, data):
	j=0
	i=0
	f = np.zeros(len(data))
	for i in range(len(data)):
		if(type(data[i][0]) == type(dummy) or type(data[i][0]) == type(dummyarray[1]) or type(data) == type(udummyarray)):
			f[i] = True
		else:
			f[i] = False
	print(f)
	#Runden
	for i in range(len(data)):
		if(f[i]):
			for j in range(data[0]):
				sdevs(data[i][j]) = rts(sdevs(data[i][j]))
				noms(data[i][j]) = round((noms(data[i][j]), -int(m.floor(m.log10(abs(x))))))
							
			
	output = open(name, 'w')
	output.write(r'\begin{table}[h]' + '\n' + r'\centering' + '\n' + r'\caption{CAPTION}' + '\n' +r'\sisetup{%uncertainty-seperator = {\,},'+'\n'+r'table-number-alignment = center,'+'\n'+'table-unit-alignment = center,'+'\n'+'%table-figures-integer = 1,'+'\n'+'%table-figures-decimal = 1,'+'\n'+'table-auto-round = true'+'\n'+'}'+'\n'+ r'\begin{tabular}{ ')
	for i in range(len(data)):
		if(f[i]):
			output.write(r'S[table-format= 3.1]'+'\n'+' @{\,$\pm{}$\,} '+'\n' + r' S[table-format= 3.1] ')
		else:
			output.write(r' S[table-format= 3.1] '+'\n')
	output.write(r'}' + '\n' + r'\toprule' + '\n')
	
	for i in range(len(data)):
		if(i < (len(data)-1)): 
			if(f[i]):
				output.write(r'\multicolumn{2}{c}{TITLE}'+'\n'+'&')
			else:
				output.write(r'{$\text{Title}$}'+'\n'+'&')
		else:
			if(f[i]):
				output.write(r'\multicolumn{2}{c}{TITLE} \\'+'\n')
			else:
				output.write(r'{$\text{Title}$} \\'+'\n')
	output.write(r' \midrule' + '\n')
	
	#Tabelle
	for j in range(len(data[0])):
		i = 0
		while i <= len(data)-1:
				if(f[i]):
					if(i == len(data)-1):
						output.write(str(data[i][j].n) + '&' + str(data[i][j].s) + r'\\' + '\n')
					else:
						output.write(str(data[i][j].n)+ '&' + str(data[i][j].s) + '&')
					i = i+1
				else:
					if(i == len(data)-1):
						output.write(str(data[i][j]) + r'\\' + '\n')
					else:
						output.write(str(data[i][j]) + '&')
					i = i+1
	#Tabelle Ende
	output.write(r'\bottomrule' + '\n' + r'\end{tabular}' + '\n' + r'\label{tab:LABEL}' + '\n' + r'\end{table}')
	output.close()

def make_table(columns, figures=None):
    assert hasattr(columns[0],'__iter__'), "Wenn nur eine Zeile von Daten vorliegt, funktioniert zip nicht mehr; die Elemente von columns müssen Listen sein, auch wenn sie ihrerseits nur ein Element enthalten."

    if figures is None:
        figures = [None] * len(columns)

    cols = []
    for column, figure in zip(columns, figures):
        assert (type(column) != str), "Hier ist ein einzelner String übergeben worden. Baue daraus eine Liste und alles ist gut ( column = [string] )."
        if (type(column) == list):
            col = zip(*zip(column))     # hard to find this kind of code... this will unzip the list column, ref: https://docs.python.org/3/library/functions.html#zip
            cols.extend(col)
        elif np.any(stds(column)):
            if figure is None:
                figure = ''
            col = list(zip(*['{0:.{1:}uf}'.format(x, figure).split('+/-') for x in column]))
        else:
            col = list(zip(*[['{0:.{1:}f}'.format(x, figure)] for x in noms(column)]))
        cols.extend(col)

    max_lens = [max(len(s) for s in col) for col in cols]
    cols = [['{0:<{1:}}'.format(s, ml) for s in col] for col, ml in zip(cols, max_lens)]

    rows = list(itertools.zip_longest(*cols))

    return (r' \\' + '\n').join([' & '.join(s for s in row if s is not None) for row in rows]) + r' \\'

def make_table(columns, figures=None):
    if figures is None:
        figures = [None] * len(columns)

    cols = []
    for column, figure in zip(columns, figures):
        if np.any(stds(column)):
            if figure is None:
                figure = ''
            col = list(zip(*['{0:.{1:}uf}'.format(x, figure).split('+/-') for x in column]))
        else:
            col = list(zip(*[['{0:.{1:}f}'.format(x, figure)] for x in noms(column)]))
        cols.extend(col)

    max_lens = [max(len(s) for s in col) for col in cols]
    cols = [['{0:<{1:}}'.format(s, ml) for s in col] for col, ml in zip(cols, max_lens)]

    rows = list(itertools.zip_longest(*cols))

    return (r' \\' + '\n').join([' & '.join(s for s in row if s is not None) for row in rows]) + r' \\'

params_gitterkonstante = ucurve_fit(
    reg_linear, sin_phi_helium, lambda_helium)

g, offset = params_gitterkonstante                  # g in m, offset Einheitenfrei
write('build/gitterkonstante.tex', make_SI(g * 1e9, r'\nano\meter', figures=1))
write('build/offset.tex', make_SI(offset * 1e9, r'\nano\meter', figures=1))
write('build/Tabelle_messdaten_kalium.tex', make_table([phi_kalium*180/np.pi],[1]))
write('build/Tabelle_messdaten_natrium.tex', make_table([phi_natrium*180/np.pi],[1]))
write('build/Tabelle_messdaten_rubidium.tex', make_table([phi_rubidium*180/np.pi],[1]))

##### PLOT lineare Regression #####
t_plot = np.linspace(np.amin(sin_phi_helium), np.amax(sin_phi_helium), 2)
plt.xlim(t_plot[0] - 1 / np.size(sin_phi_helium) * (t_plot[-1] - t_plot[0]),
         t_plot[-1] + 1 / np.size(sin_phi_helium) * (t_plot[-1] - t_plot[0]))
plt.plot(t_plot,
         reg_linear(t_plot, *noms(params_gitterkonstante))* 1e9,
         'b-', label='Fit')
plt.plot(sin_phi_helium,
         lambda_helium * 1e9,
         'rx', label='Messdaten')
plt.ylabel(r'$\lambda \:/\: \si{\nano\meter}$')
plt.xlabel(r'$\sin(\varphi)$')
plt.legend(loc='best')
plt.tight_layout(pad=0, h_pad=1.08, w_pad=1.08)
plt.savefig('build/aufgabenteil_a_plot.pdf')
plt.clf()
#### Ende Plot ####

#### TABELLE ####
write('build/Tabelle_a.tex', make_table([lambda_helium*1e9, -phi_helium, -sin_phi_helium],[1, 3, 3]))
write('build/Tabelle_a_texformat.tex', make_full_table(

Begin der Auswertung zum Dopplereffekt

"""

# Laden der Daten zur Bestimmung der Geschwindigkeit
## Laden der Strecke
l = np.loadtxt("Messdaten/Strecke.txt", unpack=True)
l_err = np.loadtxt("Messdaten/Fehler_Strecke.txt")

### Fehler behaftete Messwerte
ul = unp.uarray(l, [l_err]*len(l))

### Mittelwert
ul_avr = Umean(ul)
ul_avr = ufloat(noms(ul_avr), stds(ul_avr))
ul_avr *= 1e-02  # [cm] --> [m]

## Laden der Zeiten in den verschiedenen Gängen
G, t_h1, t_h2, t_r1, t_r2 = np.loadtxt("Messdaten/Zeiten.txt", unpack=True)
t_err = np.loadtxt("Messdaten/Fehler_Zeiten.txt")

### Fehlerbehaftete Messwerte
ut_h1 = unp.uarray(t_h1, [t_err]*len(t_h1))
ut_h2 = unp.uarray(t_h2, [t_err]*len(t_h2))
ut_r1 = unp.uarray(t_r1, [t_err]*len(t_r1))
ut_r2 = unp.uarray(t_r2, [t_err]*len(t_r2))

### Mittelwerte der Zeiten

uT_h_avr = unp.uarray(np.zeros(len(G)), np.zeros(len(G)))

def mean(values, axis=0):
    return unp.uarray((np.mean(noms(values), axis=axis), scipy.stats.sem(noms(values), axis=axis)))

def do_job_a(filename, error, j, filename_out=None):
    # Einlesen der Messdaten
    P, Delta_f_30, Delta_f_15, Delta_f_60 = np.genfromtxt(filename, unpack=True)

    #
    di = [7, 10, 16]
    colors = ["rx", "bx", "gx"]

    Delta_f_30_error = Delta_f_30 * error
    Delta_f_30 = unp.uarray(Delta_f_30, Delta_f_30_error)
    Delta_f_15_error = Delta_f_15 * error
    Delta_f_15 = unp.uarray(Delta_f_15, Delta_f_15_error)
    Delta_f_60_error = Delta_f_60 * error
    Delta_f_60 = unp.uarray(Delta_f_60, Delta_f_60_error)

    v = unp.uarray(np.zeros(3), np.zeros(3))
    v[0] = c_L / 2 / nu_0 * Delta_f_30 / np.cos(alpha[0])
    v[1] = c_L / 2 / nu_0 * Delta_f_15 / np.cos(alpha[1])
    v[2] = c_L / 2 / nu_0 * Delta_f_60 / np.cos(alpha[2])

    v_mean = mean([v[0], v[1], v[2]], 0)

    # TABLES
    write(
        "build/Tabelle_a_" + str(di[j]) + ".tex",
        make_table([P, Delta_f_30, Delta_f_15, Delta_f_60, v[0], v[1], v[2], v_mean], [0, 1, 1, 1, 1, 1, 1, 1]),
    )
    write(
        "build/Tabelle_a_" + str(di[j]) + "_texformat.tex",
        make_full_table(
            r"Messdaten und daraus errechnete Geschwindikgiet für $\d_i = $" + str(di[j]) + r"$\si{\milli\meter}$.",
            "table:A" + str(j),
            "build/Tabelle_a_" + str(di[j]) + ".tex",
            [1, 2, 3, 4, 5, 6, 7],
            [
                r"$\frac{P}{P_\text{max}} \:/\: \si{\percent}$",
                r"$\Delta f_{30°} \:/\: \si{\hertz}$",
                r"$\Delta f_{15°} \:/\: \si{\hertz}$",
                r"$\Delta f_{60°} \:/\: \si{\hertz}$",
                r"$v_{30°} \:/\: \si{\meter\per\second}$",
                r"$v_{15°} \:/\: \si{\meter\per\second}$",
                r"$v_{60°} \:/\: \si{\meter\per\second}$",
                r"$\overline{v} \:/\: \si{\meter\per\second}$",
            ],
        ),
    )

    # Plotting
    plt.figure(1)
    y = Delta_f_30 / np.cos(alpha[0])
    plt.errorbar(
        noms(v[0]),
        noms(y),
        fmt=colors[j],
        xerr=stds(v[0]),
        yerr=stds(y),
        label=r"$d_i = " + str(di[j]) + r"\si{\milli\meter}$",
    )

    plt.figure(2)
    y = Delta_f_15 / np.cos(alpha[1])
    plt.errorbar(
        noms(v[1]),
        noms(y),
        fmt=colors[j],
        xerr=stds(v[1]),
        yerr=stds(y),
        label=r"$d_i = " + str(di[j]) + r"\si{\milli\meter}$",
    )

    plt.figure(3)
    y = Delta_f_60 / np.cos(alpha[2])
    plt.errorbar(
        noms(v[2]),
        noms(y),
        fmt=colors[j],
        xerr=stds(v[2]),
        yerr=stds(y),
        label=r"$d_i = " + str(di[j]) + r"\si{\milli\meter}$",
    )

    i = 1
    if filename_out:
        for name in filename_out:
            plt.figure(i)
            plt.xlabel(r"$v \:/\: \si{\meter\per\second}$")
            plt.ylabel(r"$\Delta\nu / \cos{\alpha} \:/\: \si{\kilo\volt}$")
            plt.legend(loc="best")
            plt.tight_layout(pad=0, h_pad=1.08, w_pad=1.08)
            plt.savefig(name)
            i += 1

print("Eta, bestimmt mit der kleinen Kugel, ist in Pa*sec:")
print(etaK)
print("Mit diesem Eta ist die Apparaturkonstante in willkuerlichen Einheiten:")
print(apparatG)
print("")
print("Viskositaet in Pa*sec:")
print(etaG_2) 
print("")
print("Viskositaet in mPa*sec, vgl. Literatur:")
print(etaG_2*1000) 
print("")
print("Geschwindigkeiten in m/s:")
print((0.1/fallzeit))
print("")
print("Reynoldzahlen:")
print(Re)

np.savetxt("Auswertung/Viskositaeten.txt", np.array([noms(etaG_2),stds(etaG_2)]).T)
np.savetxt("Auswertung/Reynoldzahlen.txt", np.array([noms(Re),stds(Re)]).T)
np.savetxt("Auswertung/Geschwindigkeiten.txt", np.array([noms((0.1/fallzeit)),stds((0.1/fallzeit))]).T)
np.savetxt("Auswertung/Zeiten.txt", np.array([noms((fallzeit)),stds((fallzeit))]).T)

def geradeF(x, b):
    return noms(param_VI_g_M) * x + b

v_8_0 = 88 # s/mm
U_8 = 185

v_9_auf, v_9_ab = np.genfromtxt('messdaten/9_Tropfen.txt', unpack=True)
v_9_auf_mittel = ufloat(np.mean(v_9_auf), np.std(v_9_auf))  # s/mm
v_9_ab_mittel = ufloat(np.mean(v_9_ab), np.std(v_9_ab))     # s/mm
v_9_0 = 112 # s/mm
U_9 = 184

v_10_auf, v_10_ab = np.genfromtxt('messdaten/10_Tropfen.txt', unpack=True)
v_10_auf_mittel = ufloat(np.mean(v_10_auf), np.std(v_10_auf))  # s/mm
v_10_ab_mittel = ufloat(np.mean(v_10_ab), np.std(v_10_ab))     # s/mm
v_10_0 = 59.4 # s/mm
U_10 = 286

v_auf_mittel_nom=[ noms(v_2_auf_mittel), noms(v_3_auf_mittel), noms(v_4_auf_mittel), noms(v_5_auf_mittel), noms(v_6_auf_mittel), noms(v_7_auf_mittel), noms(v_8_auf_mittel), noms(v_9_auf_mittel), noms(v_10_auf_mittel)]

v_auf_mittel_stds=[ stds(v_2_auf_mittel), stds(v_3_auf_mittel), stds(v_4_auf_mittel), stds(v_5_auf_mittel), stds(v_6_auf_mittel), stds(v_7_auf_mittel), stds(v_8_auf_mittel), stds(v_9_auf_mittel), stds(v_10_auf_mittel)]

v_ab_mittel_nom=[ noms(v_2_ab_mittel), noms(v_3_ab_mittel), noms(v_4_ab_mittel), noms(v_5_ab_mittel), noms(v_6_ab_mittel), noms(v_7_ab_mittel), noms(v_8_ab_mittel), noms(v_9_ab_mittel), noms(v_10_ab_mittel)]

v_ab_mittel_stds=[ stds(v_2_ab_mittel), stds(v_3_ab_mittel), stds(v_4_ab_mittel), stds(v_5_ab_mittel), stds(v_6_ab_mittel), stds(v_7_ab_mittel), stds(v_8_ab_mittel), stds(v_9_ab_mittel), stds(v_10_ab_mittel)]

# v_0 = array(1/v_2_0, 1/v_3_0, 1/v_4_0, 1/v_5_0, 1/v_6_0, 1/v_7_0, 1/v_8_0, 1/v_9_0, 1/v_10_0)
v_0 = np.genfromtxt('messdaten/V_0.txt', unpack=True)
v_0 = 1/v_0
U = np.genfromtxt('messdaten/Spannung.txt', unpack=True)
v_auf = unp.uarray(v_auf_mittel_nom, v_auf_mittel_stds)
v_auf = 1/v_auf
v_ab  = unp.uarray(v_ab_mittel_nom, v_ab_mittel_stds)
v_ab = 1/v_ab

#
# plt.plot(t_plot, reg.reg_linear(t_plot, *noms(params)), 'b-', label='Fit')
# plt.xlim(t_plot[0], t_plot[-1])
# # plt.xlabel(r'$t \:/\: \si{\milli\second}$')
# # plt.ylabel(r'$U \:/\: \si{\kilo\volt}$')
# plt.legend(loc='best')
# plt.tight_layout(pad=0, h_pad=1.08, w_pad=1.08)
# plt.savefig('build/test-plot.pdf')

# Ablesen der Grenzfrequenzen und Umrechnen
f_gr = 9
f_gr1 = 5.60
f_gr2 = 14.73
f_gr3 = 18.15

f_gr = np.exp(reg.reg_linear(f_gr, noms(m1), noms(b1)))
f_gr1 = np.exp(reg.reg_linear(f_gr1, noms(m2), noms(b2)))
f_gr2 = np.exp(reg.reg_linear(f_gr2, noms(m2), noms(b2)))
f_gr3 = np.exp(reg.reg_linear(f_gr3, noms(m2), noms(b2)))

write('build/Z_w_gr.tex', make_SI(Wellenwiderstand(2*np.pi*f_gr), r'\ohm', figures=0))
write('build/f_mess.tex', make_SI(f_gr*1e-3, r'\kilo\hertz', 'e-3',figures=1))
write('build/f1_mess.tex', make_SI(f_gr1*1e-3, r'\kilo\hertz', 'e-3',figures=1))
write('build/f2_mess.tex', make_SI(f_gr2*1e-3, r'\kilo\hertz', 'e-3',figures=1))
write('build/f3_mess.tex', make_SI(f_gr3*1e-3, r'\kilo\hertz', 'e-3',figures=1))

# Theoriewerte der Grenzfrequenzen
w_th = 2 / np.sqrt(L*C1)
w1_th = np.sqrt(2/(L*C1))
w2_th = np.sqrt(2/(L*C2))
w3_th = np.sqrt( 2/L * (C1+C2)/(C1*C2) )

#plt.show()
plt.clf()


# Subtraktion der Langlebigen Zerfälle
N_Rh_lang = np.array(np.exp(-popt_Rh_2[0]*T_Rh[:10])*exp(popt_Rh_2[1]))
#print(N_Rh_lang)

# Fehler der Werte
n_Rh_lang_err = np.array([m.sqrt(x) for x in N_Rh_lang])
# Fehlerbehaftete Größe
N_Rh_lang_err = unp.uarray(N_Rh_lang, n_Rh_lang_err)

N_Rh_kurz_err= np.array(np.subtract(N_Rh_err[:10], N_Rh_lang_err))
#print(N_Rh_kurz_err)
lnN_Rh_kurz = np.array([m.log(x) for x in noms(N_Rh_kurz_err)])
#print(lnN_Rh_kurz)



# Regression kurzlebig
popt_Rh_1, pcov_Rh_1 = curve_fit(f_gerade, T_Rh[:10], lnN_Rh_kurz) #sigma=stds(N_Rh_kurz_err)/noms(N_Rh_kurz_err))
errors_Rh_1 = np.sqrt(np.diag(pcov_Rh_1))
param_a_Rh_1 = ufloat(popt_Rh_1[0], errors_Rh_1[0])
param_b_Rh_1 = ufloat(popt_Rh_1[1], errors_Rh_1[1])

print("Parameter Exp. Regression a,b :", param_a_Rh_1, param_b_Rh_1)



# Plot der Messwerte für t < t*

uP_2 = unp.uarray(P_2, P_dim * [P_2_err])

#%%

## Plots der Temperaturverläufe

# 16 t- Werte in Sekunden (0, 90, 180, ..., 1350)
t = np.linspace(0, (T_dim - 1) * DELTA_T, num=T_dim)

# "Kontinuierliche" Werte von -10 bis 2000
x = np.linspace(-10, 2000, num=1000)


# Bearbeitung von T_1
    #Berechnung der Fit-Parameter
poptI_T1, pcovI_T1 = curve_fit(T_FitI, t, noms(uT_1), sigma=stds(uT_1))
poptII_T1, pcovII_T1 = curve_fit(T_FitII, t, noms(uT_1), sigma=stds(uT_1))
poptIII_T1, pcovIII_T1 = curve_fit(T_FitIII, t, noms(uT_1), sigma=stds(uT_1))

    # Berechnung der Fit-Fehler
errorI_T1 = np.sqrt(np.diag(pcovI_T1))
#errorII_T1 = np.sqrt(np.diag(pcovII_T1))
#errorIII_T1 = np.sqrt(np.diag(pcovIII_T1))

uA_1 = ufloat(poptI_T1[0], errorI_T1[0]) 
uB_1 = ufloat(poptI_T1[1], errorI_T1[1]) 
uC_1 = ufloat(poptI_T1[2], errorI_T1[2]) 



    # Plot Einstellungen

Ip_calc = uY_p/R
Im_calc = uY_m/R

n = 1
for c in uC3:
    plt.clf()
    plt.grid()
    plt.tick_params("both", labelsize=16)
    plt.xlabel("Frequenz $f\,[\mathrm{kHz}]$", fontsize=16)
    plt.xlim(2e04, 6e04)
    plt.gca().xaxis.set_major_formatter(mpl.ticker.FuncFormatter
                                       (lambda x, _: float(x * 1e-03)))
    plt.gca().yaxis.set_major_formatter(mpl.ticker.FuncFormatter
                                       (lambda x, _: float(x * 1e03)))
    plt.ylabel("Stromstärke $I\,[\mathrm{mA}]$", fontsize=16)
    plt.plot(Xrange, noms(Strom(4, Xrange * 2 * const.pi, uC, c, R, uL)))
    plt.savefig("Grafiken/Stromverlauf{}.pdf".format(str(n)))
    n += 1

plt.clf()
plt.grid()
plt.xlabel("Frequenz $f\,[\mathrm{kHz}]$")
plt.xlim(2e04, 6e04)
plt.gca().xaxis.set_major_formatter(mpl.ticker.FuncFormatter
                                    (lambda x, _: float(x * 1e-03)))
plt.ylabel("Stromstärke $I\,[\mathrm{A}]$")
plt.plot(Xrange, noms(Strom(4, Xrange * 2 * const.pi, uC, uC3[0], R, uL)))


### Noch ein paar Plots
plt.clf()

#N_zu

I_s1=3.110
I_s2=2.440
I_s3=1.294
I_s4=0.721
I_s5=0.250


Austritt_1=-unp.log((I_s1/f)*(const.h**3)/(const.e*const.m_e*(const.k**2)*(T_w1**2)))*const.k*T_w1/const.e
Austritt_2=-unp.log((I_s2/f)*(const.h**3)/(const.e*const.m_e*(const.k**2)*(T_w2**2)))*const.k*T_w2/const.e
Austritt_3=-unp.log((I_s3/f)*(const.h**3)/(const.e*const.m_e*(const.k**2)*(T_w3**2)))*const.k*T_w3/const.e
Austritt_4=-unp.log((I_s4/f)*(const.h**3)/(const.e*const.m_e*(const.k**2)*(T_w4**2)))*const.k*T_w4/const.e
Austritt_5=-unp.log((I_s5/f)*(const.h**3)/(const.e*const.m_e*(const.k**2)*(T_w5**2)))*const.k*T_w5/const.e




#Austritt_Leistung=-unp.log((I_s1/f)*(const.h**3)/(const.e*const.m_e*(const.k**2)*(T_w1**2)))*const.k*T_w1/const.e

print('Austritt_1',Austritt_1)
print('Austritt_2',Austritt_2)
print('Austritt_3',Austritt_3)
print('Austritt_4',Austritt_4)
print('Austritt_5',Austritt_5)


print('Mittelwert normal',np.mean([Austritt_1,Austritt_2,Austritt_3,Austritt_4,Austritt_5]),'+-',np.std([noms(Austritt_1),noms(Austritt_2),noms(Austritt_3),noms(Austritt_4),noms(Austritt_5)]))
#print('Mittelwert np',np.mean([Austritt_2,Austritt_1]),'+-',np.std([noms(Austritt_2),noms(Austritt_1)]))

# write('build1/offset2.tex', make_SI(1e-5*g(373.1, *noms(params_2))-1, r'\bar', '', 1))
# p2 += 1e5-g(373.1, *noms(params_2))
# params_2 = ucurve_fit(g, T2, p2)

a2, b2, c2, d2 = params_2
d2 += 48*1e3    # Korrektur Offset (laut Protokoll)
write('build1/a2.tex', make_SI(a2 * 1e-5, r'\bar\per\kelvin\tothe{3}', '', 1))
write('build1/b2.tex', make_SI(b2 * 1e-5, r'\bar\per\kelvin\tothe{2}', '', 1))
write('build1/c2.tex', make_SI(c2 * 1e-5, r'\bar\per\kelvin\tothe{1}', '', 1))
write('build1/d2.tex', make_SI(d2 * 1e-5, r'\bar', '', 1))

T_hilf = np.linspace(np.amin(T2), np.amax(T2), 5)
T_plot = np.linspace((T_hilf[0]-1/np.size(T2)*(T_hilf[-1]-T_hilf[0])), (T_hilf[-1]+1/np.size(T2)*(T_hilf[-1]-T_hilf[0])), 10)

plt.clf()
plt.plot(T_plot, g(T_plot, *noms(params_2))*1e-5, 'b-', label='Fit')
plt.plot(T2, p2*1e-5, '.r', label='Messdaten')

plt.xlim(np.amin(T_plot), np.amax(T_plot))
plt.xlabel(r'$T \:/\: \si{\kelvin}$')
plt.ylabel(r'$p \:/\: \si{\bar}$')
plt.legend(loc='best')
plt.tight_layout(pad=0, h_pad=1.08, w_pad=1.08)
plt.savefig('build1/plot2.pdf')

plt.clf()
T_plot = np.linspace(200,600,100)
plt.plot(T_plot, L(0.9, R[0], T_plot, *noms(params_2)), 'b-', label='Fit numpy curve_fit')
plt.plot(T_plot, L(0.9, R[0], T_plot, 4.32*1e-1, -4.07*1e2, 1.262*1e5, -127.82*1e5+48*1e3), 'r-', label='Fit Origin 8')
#plt.plot(T2, p2*1e-5, '.r', label='Messdaten')

write('build/Abschirmkonstante_Zink.tex', make_SI(sigma_zink, r' ', figures=2))
write('build/Abschirmkonstante_Zirkonium.tex', make_SI(sigma_zir, r' ', figures=2))

#Moseley-Diagramm

E_k = (E_zink, E_ger, E_zir)
Z   = (30,32,40) # Zn, Ge, Zr
E_k_wurzel = np.sqrt(E_k)
params = ucurve_fit(reg_linear, Z, E_k_wurzel)
m,b = params
write('build/hcRydbergonstante.tex', make_SI(4/3*m**2, r'\electronvolt', figures=1))
write('build/Rydbergonstante.tex', make_SI(4/3*m**2/(h*c), r'\per\meter', figures=1))

plt.clf
t_plot = np.linspace(25,45, 100)
plt.plot(t_plot , reg_linear(t_plot, *noms(params)), 'b-', label='Fit')
plt.plot(Z, E_k_wurzel, 'rx', label='Messdaten')
plt.xlabel(r'Kernladungszahl  $Z$')
plt.ylabel(r'$\sqrt{E_\textrm{k} \:/\: \si{\kilo\electronvolt}}$')
plt.legend(loc='best')
plt.savefig("build/Moseley_Diagramm.pdf")
plt.close

################################ FREQUENTLY USED CODE ################################
#
########## IMPORT ##########
# t, U, U_err = np.genfromtxt('data.txt', unpack=True)
# t *= 1e-3


########## ERRORS ##########

# Umrechnung der Messwerte
L_1_off_err = L_1_err[:]
L_1_err = AbstandOhneOffset(L_1_err)
#print("Abstände ohne", L_1_err)

# Laden der Intensitäten, l nicht gebraucht
l , J_1 = np.loadtxt("Messdaten/Abstand_Intensitaet_1.txt", unpack=True)

# Fehlerbehaftete Intensitäten
J_1_err = unp.uarray(J_1, [j_err1 if J_1[i] < 10 else j_err2 for i in range(len(J_1))])


# Auftragen des Stroms gegen die Intensität
# linearer Fit der Messwerte
func_gerade = lambda x,a,b: a*x+b
popt_1, pcov_1 = curve_fit(func_gerade, noms(J_1_err), noms(I_k_1_err))
error_1 = np.sqrt(np.diag(pcov_1))
param_a_1 = ufloat(popt_1[0], error_1[0])
param_b_1 = ufloat(popt_1[1], error_1[1])
print("Kurzschlussstrom-Fit:")
print("Steigung:", param_a_1)
print("Y-Achsenabschnitt:", param_b_1)


# Plot der Wertepaare (J/I)
plt.plot(noms(J_1_err), noms(I_k_1_err), "xr", label="Messwerte")

# Plot der Fit-Gerade
X = np.linspace(1,30, 300)
plt.plot(X, func_gerade(X, *popt_1), color="gray", label="Regressionsgerade")

#Messung 2

def Brech(z,lamb,T,p0,b,T0,deltaP):
    n=1+z*lamb*T*p0/(2*b*T0*deltaP)
    return n

p1,p2,ticks =np.genfromtxt('Messung2.txt', unpack = True)
deltaP=p2-p1
p0=1.0132
T0=273.15
b=50e-3
T=293.15


n=unp.uarray(noms(Brech(ticks,L,T,p0,b,T0,deltaP)),stds(Brech(ticks,L,T,p0,b,T0,deltaP)))

# n=Brech(ticks,L,T,p0,b,T0,deltaP)
print(n,'n Luft')
print(np.mean(noms(n)),stats.sem(noms(n)),'Luft np.mean(noms(n)),stats.sem(noms(n))')
nL=ufloat(np.mean(noms(n)),stats.sem(noms(n)))




#Messung3
p1,p2,ticks =np.genfromtxt('Messung3.txt', unpack = True)
deltaP=p2-p1

n=unp.uarray(noms(Brech(ticks,L,T,p0,b,T0,deltaP)),stds(Brech(ticks,L,T,p0,b,T0,deltaP)))
print(n,'n 1-Butylen')