def plot_u(cal_file, fm_file, offset_file, description, accidental_offset,
        results_file):
    M = np.genfromtxt(cal_file)
    N = np.genfromtxt(fm_file)
    O = np.genfromtxt(offset_file)

    i_helm = M[:,1] #current applied to helmholtz for calibration measurement
    b_helm = M[:,2] #field applied to helmholtz coil for calibration measurement
    p, cov = np.polyfit(i_helm, b_helm, 1,  cov=True) #fit a line to calibration measurement so that we get a calibration

    i_fm = N[:,1] #current applied to helmmholtz for shielding measurement
    b_fm = unumpy.uarray(N[:,2],0.0005) - accidental_offset #field measured inside of ferromagnet shield


    B_earth = np.polyval(p,0) #We get the Earths magnetic field from i=0 of the Helmholtz calibration
    B_fm_no_i = unumpy.uarray(np.mean(O[:,2]), np.std(O[:,2])) #Get average and error for initial magnetization

    mag = B_fm_no_i - B_earth #initial magnetization is the field inside of the ferromagnet before any field is applied minus the earths magnetic field 

    Bin = b_fm - mag #internal magnetization is the measured internal field minus the initial magnetization. This correction might not be necessary for a soft ferromagnet
    
    Bext = unumpy.uarray(np.polyval(p,i_fm), 0.0005) #external field
    Bext_nom = unumpy.nominal_values(Bext)
    Bext_err = unumpy.std_devs(Bext)

    B = Bext/Bin
    c = a/b
    u=(-2*B + c**2 - 2*unumpy.sqrt(B**2 - B*c**2 - B + c**2) + 1)/(c**2 - 1)

    u_nom = unumpy.nominal_values(u)
    u_err = unumpy.std_devs(u)

    #cakculate uerr with just point to point uncertainties. I define this as just
    #uncertainty from the field measurements
    u_pp=(-2*B + c.nominal_value**2 - 2*unumpy.sqrt(B**2 - B*c.nominal_value**2 - B + c.nominal_value**2) + 1)/(c.nominal_value**2 - 1)
    
    #calculate uerr from just geometry uncertainties
    u_geom=(-2*unumpy.nominal_values(B) + c**2 - 2*unumpy.sqrt(unumpy.nominal_values(B)**2 - unumpy.nominal_values(B)*c**2 - unumpy.nominal_values(B) + c**2) + 1)/(c**2 - 1)

    ##obtain uncertainties from field
    u_err_pp = unumpy.std_devs(u_pp)

    ##obtain uncertainties from geometry
    u_err_geom = unumpy.std_devs(u_geom)

    with open(results_file, "w") as myfile:
        myfile.write('#Bext, sig_Bext, ur, sig_ur, sig_ur_pp, sig_ur_corr\n')
        for j in range(0, len(u_nom)):
            myfile.write('%s\t%s\t%s\t%s\t%s\t%s\n' %(
                Bext_nom[j], Bext_err[j],
                u_nom[j], u_err[j], u_err_pp[j], u_err_geom[j]))
    

    plt.errorbar(Bext_nom, u_nom, u_err, marker = '.', label = description)

def plot_u(cal_file, fm_file,  description, accidental_offset,
        results_file):
    M = np.genfromtxt(cal_file) #turn calibration file into a matrix
    N = np.genfromtxt(fm_file) #turn fm_scan file into a matrix


    i_helm = M[:,1] #current applied to helmholtz for calibration measurement
    b_helm = M[:,2] #field applied to helmholtz coil for calibration measurement
    p, cov = np.polyfit(i_helm, b_helm, 1,  cov=True) #fit a line to calibration measurement so that we get a calibration

    
    i_fm = N[:,1] #current applied to helmmholtz for shielding measurement
    Bin = unumpy.uarray(N[:,2],0.0005) - accidental_offset #field measured inside of ferromagnet shield


    Bin_nom = unumpy.nominal_values(Bin) 
    Bin_err = unumpy.std_devs(Bin)

    Bext = unumpy.uarray(np.polyval(p,i_fm), 0.0005) #external field
    Bext_nom = unumpy.nominal_values(Bext)
    Bext_err = unumpy.std_devs(Bext)

    B = Bin/Bext #ratio of internal to external field

    #calculate permeability
    u = (B*c**2 + B -2 -2*unumpy.sqrt(B**2*c**2 - B*c**2 - B + 1))/(B*c**2-B) 
    print(u)
    u_nom = unumpy.nominal_values(u)
    u_err = unumpy.std_devs(u)

    #cakculate uerr with just point to point uncertainties. I define this as just
    #uncertainty from the field measurements
    u_pp = (B*c.n**2 + B -2 -2*unumpy.sqrt(B**2*c.n**2 - B*c.n**2 - B +
        1))/(B*c.n**2-B)

    u_err_pp = unumpy.std_devs(u_pp)

    #calculate uerr from just geometry uncertainties
    u_geom = (unumpy.nominal_values(B)*c**2 + unumpy.nominal_values(B) -2 -2*unumpy.sqrt(unumpy.nominal_values(B)**2*c**2 - unumpy.nominal_values(B)*c**2 - unumpy.nominal_values(B) +
        1))/(unumpy.nominal_values(B)*c**2-unumpy.nominal_values(B))

    u_err_geom = unumpy.std_devs(u_geom)


    #write results onto a text file
    with open(results_file, "w") as myfile:
        myfile.write('#Bext, sig_Bext, Bi, sig_Bi, ur, sig_ur, sig_ur_pp, sig_ur_corr\n')
        for j in range(0, len(u_nom)):
            myfile.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' %(
                Bext_nom[j], Bext_err[j], Bin_nom[j], Bin_err[j],
                u_nom[j], u_err[j], u_err_pp[j], u_err_geom[j]))
    
    plt.errorbar(Bext_nom, u_nom, u_err, marker = '.', label = description)

def calc_mu(Bin, Bout, radius_inner, radius_outer):

    # ratio of inner and outer radius
    radius_ratio = radius_inner / radius_outer
    print( radius_ratio )
    # ratio of internal to external field
    B_ratio = Bin / Bout
    # convert from Series object to Numpy Array object
    B_ratio = B_ratio.values

    # If value under square root becomes neagtive, set B_ratio to NaN
    # (this seems to happen with Argonne MRI measurements at low fields)
    B_ratio[ (B_ratio**2) * (radius_ratio**2) - B_ratio * (radius_ratio**2) - B_ratio + 1 < 0 ] = np.nan

    # Calculate permeability. Here, use both uncertainties from field measurements and uncertainties from geometry, i.e. radius measurements.
    mu = ( B_ratio * (radius_ratio**2)
           + B_ratio
           - 2
           -2 * unumpy.sqrt( (B_ratio**2) * (radius_ratio**2) - B_ratio * (radius_ratio**2) - B_ratio + 1 )
           ) / ( B_ratio * (radius_ratio**2) - B_ratio )

    # store nominal values of mu in separate array
    mu_val = unumpy.nominal_values(mu)

    # store combined uncertainties of mu in separate array
    mu_err = unumpy.std_devs(mu)

    # Calculate uncertainties of mu values from just field measurement uncertainties (= point-to-point fluctuations). Ignore geometry (=radius) uncertainties.
    mu_pp = ( B_ratio * (radius_ratio.n**2)
              + B_ratio
              - 2
              -2 * unumpy.sqrt( (B_ratio**2) * (radius_ratio.n**2) - B_ratio * (radius_ratio.n**2) - B_ratio + 1 )
              ) / ( B_ratio * (radius_ratio.n**2) - B_ratio )

    # store point-to-point uncertainties of mu in separate array
    mu_err_pp = unumpy.std_devs(mu_pp)

    # Calculate uncertainties of mu values from just geometry uncertainties (= systematic uncertainty, i.e. all points move together). Ignore field uncertainties.
    B_ratio_n = unumpy.nominal_values( B_ratio )
    mu_geom = ( B_ratio_n * (radius_ratio**2)
                + B_ratio_n
                - 2
                -2 * unumpy.sqrt( (B_ratio_n**2) * (radius_ratio**2) - B_ratio_n * (radius_ratio**2) - B_ratio_n + 1 )
                ) / ( B_ratio_n * (radius_ratio**2) - B_ratio_n )

    # store geometric uncertainties of mu in separate array
    mu_err_geom = unumpy.std_devs(mu_geom)

    return( mu_val, mu_err, mu_err_pp, mu_err_geom )

def radius(drops, viscosity):
    radi = []
    for i  in range(len(drops)):
        radius = unp.sqrt(9*viscosity*(drops[i][2]-drops[i][3])/(2*d.g*(d.density_oil-d.density_air)))
        radi.append(radius)
#    print("UNICORNS 4 EVVAA:", len(radi))
    return radi

def stress_from_strain_inplane(ex,ey,exy=0,ec=None,dex=None):
    if dex is None:
        if ec is None:
            raise 
        else:
            E,v = ec
    else:
        E,v = __Ev_from_DEX(*dex)
        
    f = E/(1-v**2)

    stress_x=f*(ex+v*ey)    
    stress_y=f*(ey+v*ex)
    stress_xy=f*(1-v)*exy
    
    smean = (stress_x+stress_y)/2.
    sdiff = unumpy.sqrt(((stress_x-stress_y)/2)**2.+stress_xy**2)
    stress_1 = smean+sdiff
    stress_2 = smean-sdiff
    rotation = unumpy.arctan2(stress_xy,sdiff)/2
    
    #s1,s2,srot = __principal_inplane(nominal_value(stress_x),nominal_value(stress_y),nominal_value(stress_xy))  
    
    return dict(stress_x=stress_x,
                stress_y=stress_y,
                stress_xy=stress_xy,
                stress_1=stress_1,
                stress_2=stress_2,
                stress_rotation=rotation)

def charge(E, drops, viscosity):
    q_charge = []
    for i in range(len(drops)):
        charge = 3 * np.pi * viscosity * unp.sqrt(9*viscosity*(drops[i][2]-drops[i][3])/(4*d.g*(d.density_oil-d.density_air))) *(drops[i][2]+drops[i][3])/E[i]
        q_charge.append(charge)
#    print("HEEEEELPP:", len(q_charge))
    return q_charge

def analyze_spektrallinien(fileprefix, figindex, crstl, sl, d=None, y=None):

    data = np.append(np.loadtxt(fileprefix+'.b.1.txt', skiprows=1), np.loadtxt(fileprefix+'.b.2.txt', skiprows=1), axis=0)

    b, n = data[:,0], data[:,1]
    n = unp.uarray(n, np.sqrt(n*20)/20)
    
    sl = [ [(b >= bounds[0]) & (b <= bounds[1]) for bounds in sl_row] for sl_row in sl]

    def fit_gauss(x, m, s, A, n_0):
        return A/np.sqrt(2*const.pi)/s*np.exp(-((x-m)**2)/2/(s**2))+n_0
    
    r = []
    
    plt.clf()
    papstats.plot_data(b,n)
    papstats.savefig_a4('3.'+str(figindex)+'.a.png')

    plt.clf()
    plt.suptitle('Diagramm 3.'+str(figindex)+u': Spektrallinien von Molybdän bei Vermessung mit einem '+crstl+'-Kristall')
    for i in range(2):
        r.append([])
        # Linie
        for k in range(2):
            # Ordnung
            b_k = b[sl[i][k]]
            n_k = n[sl[i][k]]
            xspace = np.linspace(b_k[0], b_k[-1], num=1000)
            plt.subplot(2,2,i*2+k+1)
            plt.xlim(xspace[0], xspace[-1])
            if i==1:
                plt.xlabel(u'Bestrahlungswinkel '+r'$\beta \, [^\circ]$')
            if k==0:
                plt.ylabel(u'Zählrate '+r'$n \, [\frac{Ereignisse}{s}]$')
            plt.title('$K_{'+(r'\alpha' if i==0 else r'\beta')+'}$ ('+str(k+1)+'. Ordnung)')
            papstats.plot_data(b_k, n_k)
            # Gauss-Fit
            popt, pstats = papstats.curve_fit(fit_gauss, b_k, n_k, p0=[b_k[0]+(b_k[-1]-b_k[0])/2, (b_k[-1]-b_k[0])/4, np.sum(n_k).n, n_k[0].n])
            plt.fill_between(b_k, 0, unp.nominal_values(n_k), color='g', alpha=0.2)
            FWHM = popt[1]*2*unp.sqrt(2*unp.log(2))
            plt.hlines(popt[3].n+(fit_gauss(xspace, *unp.nominal_values(popt)).max()-popt[3].n)/2, popt[0].n-FWHM.n/2, popt[0].n+FWHM.n/2, color='black', lw=2, label='$'+papstats.pformat(FWHM, label='FWHM', unit=r'^\circ')+'$')
            papstats.plot_fit(fit_gauss, popt, xspace=xspace, plabels=[r'\mu', r'\sigma', 'A', 'n_0'], punits=['^\circ', '^\circ', 's^{-1}', 's^{-1}'])
            plt.ylim(unp.nominal_values(n_k).min()-n_k[unp.nominal_values(n_k).argmin()].s, unp.nominal_values(n_k).max()+(unp.nominal_values(n_k).max()-unp.nominal_values(n_k).min()))
            plt.legend(loc='upper center', prop={'size':10})

            b_S = unc.ufloat(popt[0].n, np.abs(popt[1].n))
            print "Winkel:", papstats.pformat(b_S, unit='°', format='.2u')
            if y is None:
                r[i].append(y_bragg(b_S, n=k+1))
                print "Wellenlänge der Linie:", papstats.pformat(r[i][k]/const.pico, label='y', unit='pm', format='.2u')
            if d is None:
                r[i].append((k+1)*y[i][k]/unc.umath.sin(b_S*const.degree))
                print "Gitterkonstante:", papstats.pformat(r[i][k]/const.pico, label='a', unit='pm', format='.2u')

    papstats.savefig_a4('3.'+str(figindex)+'.png')

    return r

def bethe(E, n, z, I, m):
    print(E, n, z, I, m)
    v = unp.sqrt((2 * E / m).to_base_units().magnitude) * (u.meter/u.second)
    print(v.to('m/s'))

    ln = unp.log(((2 * m_e * v**2) / I).to_base_units().magnitude)

    a = (4 * np.pi * n * z**2) / (m_e * v**2)
    b = ((c.e * u.coulomb)**2 / (4 * np.pi * epsilon_0))**2

    return (a * b * ln).to('MeV / cm')

def experimental_form2(q_bin, dG, m, m_P, m_D, m_Pr, t0):
    Return = []
    Q = np.array([np.mean(i) for i in q_bin])

    for i in range(len(q_bin)):
        dq = q_bin[i][1] - q_bin[i][0]
        q = q_bin[i][0] + dq/2.0

        P = p(q, m_P, m_D)**3

        Return.append(
            1.0 * unp.sqrt( dG[i] / dq * (24.0*np.pi**3) / (fermi**2 * P))
        )

    return Q, unp.nominal_values(Return), unp.std_devs(Return)

def crit_gni(Re,Pr):
    """
    Calculeaza criteriul Gnielinski
    
    Intrare
    -------
    Re: criteriul Reynolds
    Pr: criteriul Prandtl
    
    Iesire
    ------
    Nu: criteriul Nusselt
    """
    
    fric=(0.79*unp.log(Re)-1.64)**(-2)
    return (fric/8.)*(Re-1000.)*Pr/(1+12.7*unp.sqrt(fric/8.)*(Pr**(2./3)-1))

def sqrt(number):
    number = Mixed(number)
    if isinstance(number.value, fr.Fraction):
        p = np.sqrt(number.value.numerator)
        q = np.sqrt(number.value.denominator)
        if (p % 1 == 0) and (q % 1 == 0):
            return Mixed(fr.Fraction(int(p), int(q)))
        else:
            return Mixed(np.sqrt(float(number.value)))
    elif isinstance(number.value, uc.UFloat):
        return Mixed(unumpy.sqrt(number.value).item())
    elif isinstance(number.value, int):
        s = np.sqrt(number.value)
        if s % 1 == 0:
            return Mixed(int(s))
        else:
            return Mixed(s)
    elif isinstance(number.value, float):
        return Mixed(float(np.sqrt(number.value)))

def Radius(auf, ab):
    return unp.sqrt((9*eta_l*((1/ab)*1e-3 - (1/auf)*1e-3))/(4*g*(rho_oel - rho_l)))

def Ladung(auf, ab, U):
    return 3*np.pi*unp.sqrt((9*eta_l*((1/ab)*1e-3 - (1/auf)*1e-3))/(4*g*(rho_oel - rho_l)))* ((1/ab)*1e-3 + (1/auf)*1e-3)/(U/d)

p, cov = np.polyfit(i_helm, b_helm, 1,  cov=True) #fit a line to calibration measurement so that we get a calibration

print(cov)
B_earth = np.polyval(p,0) #We get the Earths magnetic field from i=0 of the Helmholtz calibration
B_fm_no_i = unumpy.uarray(np.mean(O[:,2]), np.std(O[:,2])) #Get average and error for initial magnetization
print(B_earth)
print(B_fm_no_i)

mag = B_fm_no_i - B_earth #initial magnetization is the field inside of the ferromagnet before any field is applied minus the earths magnetic field 

print(mag)
Bin = b_fm - mag #internal magnetization is the measured internal field minus the initial magnetization. This correction might not be necessary for a soft ferromagnet
Bext = unumpy.uarray(np.polyval(p,i_fm), 0.0005) #external field
Bext_nom = unumpy.nominal_values(Bext)
Bext_err = unumpy.std_devs(Bext)

#calculate u_r
u=(-2*Bext*b**2 + Bin*a**2 + Bin*b**2 - 2*unumpy.sqrt(b**2*(Bext**2*b**2 - Bext*Bin*a**2 - Bext*Bin*b**2 + Bin**2*a**2)))/(Bin*(a**2 - b**2))

#obtain nominal value and uncertainty from u
u_nom = unumpy.nominal_values(u)
u_err = unumpy.std_devs(u)

#plt.errorbar(unumpy.nominal_values(Bext), u_nom, u_err)
plt.errorbar(Bext_nom, u_nom*Bext_nom, u_err, fmt = '.') #plot u_r vs B_ext
plt.errorbar(Bext_nom, u_nom*Bext_nom, u_err, color = 'b') #plot a line as well so that it doesn't look like the values are just jumping around.
plt.xlabel('$B_{ext}$ [mT]')
plt.ylabel('Relative permeability')
plt.title('Relative permeability measurement for kapton/steel powder')
plt.savefig('permeability_scan.png')

 def norm(v):
     return unumpy.sqrt(np.dot(v, v))

#R2(509.5,0.5)
#Rap(3.3e3)
#1.22627763e-01   1.63147814e+02   3.47381148e+02   5.80136682e-02
#1.15649731e-02

L=ufloat(10.11e-3,0.03e-3)
C=ufloat(2.098e-9,0.006e-9)
Reff=2*L*pi2ny

print ('pi2ny')
print (pi2ny)
#pi2ny(-5.53+/-0.16)e+03

print ('Reff')
print(Reff)
#Reff-111.9+/-3.3
Rap=unp.sqrt(4*L/C)
print('Rap')
print (Rap)
#Rap4390.4+/-9.0



Tex=1/pi2ny
print('Tex')
print(Tex)
#Tex-0.000181+/-0.000005
#print ('Refft')

#print ('Text')

Csp=ufloat(37,1)*10**-12
L=ufloat(31.90,0.05)*10**-3
R= ufloat(50,0)

uCkopplung=unp.uarray([9.99,8.00,6.47,5.02,4.00,3.00,2.03,1.01],[9.99/100,8.00/100,6.47/100,5.02/100,4.00/100,3.00/100,2.03/100,1.01/100])
Ckopplung=uCkopplung*10**(-9)

Cgesamt=(C*Ckopplung)/(2*C+Ckopplung)+Csp

vau1m=np.array([30.56,30.57,30.57,30.57,30.57,30.57,30.58,30.58])*1000
vau2m=np.array([32.85,33.38,33.98,34.88,35.86,37.40,40.12,47.23])*1000

alpham=np.array([14,11,9,7,6,5,4,2])

# Resonanzfrequenz
fres=unp.sqrt(1/(L*(C+Csp))-R**2/(4*L**2))/(2*np.pi)

# Fundamentalfrequenz
vau1 =1/(2*np.pi*unp.sqrt(L*(C+Csp)))
vau2 = 1/(2*np.pi*unp.sqrt(L*Cgesamt))

#Abweichung der Fundamentalfrequenzen
vau1dev=(vau1-vau1m)*100/vau1
vau2dev=(vau2-vau2m)*100/vau2

#Frequenzverhältnisse
alpha=-(vau1+vau2)/(2*(vau1-vau2))
alphav=(alpha-alpham)*100/alpham

# Ausgabe
print("""

print(mag)
Bin = (
    b_fm - mag
)  # internal magnetization is the measured internal field minus the initial magnetization. This correction might not be necessary for a soft ferromagnet
Bext = unumpy.uarray(np.polyval(p, i_fm), 0.0005)  # external field
Bext_nom = unumpy.nominal_values(Bext)
Bext_err = unumpy.std_devs(Bext)

print(Bext - Bin)
# calculate u_r
u = (
    -2 * Bext * b ** 2
    + Bin * a ** 2
    + Bin * b ** 2
    - 2 * unumpy.sqrt(b ** 2 * (Bext ** 2 * b ** 2 - Bext * Bin * a ** 2 - Bext * Bin * b ** 2 + Bin ** 2 * a ** 2))
) / (Bin * (a ** 2 - b ** 2))

# obtain nominal value and uncertainty from u
u_nom = unumpy.nominal_values(u)
u_err = unumpy.std_devs(u)

# plt.errorbar(unumpy.nominal_values(Bext), u_nom, u_err)
plt.errorbar(Bext_nom, u_nom, u_err, fmt=".")  # plot u_r vs B_ext
plt.errorbar(
    Bext_nom, u_nom, u_err, color="b"
)  # plot a line as well so that it doesn't look like the values are just jumping around.
plt.xlabel("$B_{ext}$ [mT]")
plt.ylabel("Relative permeability")
plt.title("Relative permeability measurement for kapton/steel powder")
plt.savefig("permeability_scan.png")

def j(dt, ft):
    return ((-1/(dt*ft))+unp.sqrt((1/(dt*ft))**2+1))

plt.plot(f, uc , 'b.', label="Kondensatorspannung")

plt.xlabel(r'$Frequenz / Hz$')
plt.ylabel(r'$Uc/U$')
plt.xscale('log')
plt.tight_layout()
plt.legend(loc='best')
plt.savefig('Kondensatorspannungen.pdf')

L=ufloat(10.11e-3,0.03e-3)
C=ufloat(2.098e-9,0.006e-9)
R=ufloat(509.5,0.5)
Re=R+50

rdiffomega=Re/L
romegares=unp.sqrt(1/(C*L)-(Re/L)**2/2)
print ('rdiffomega')
print (rdiffomega)
print ('romegares')
print (romegares)
rguete=romegares/rdiffomega
print ('rguete')
print (rguete)
eomegadiff=ufloat(10000,0.1)
eomegares=ufloat(33500,0.1)
eguete=1/eomegadiff*eomegares
print (eguete)
#(rdiffomega)(5.040+/-0.016)e+04
#(romegares)(2.171+/-0.004)e+05
#(rguete)4.309+/-0.010