def other(g=True):
    freq = unumpy.uarray([100, 500, 1000, 5000, 10000, 50000], np.array([100, 500, 1000, 5000, 10000, 50000]) * 0.01)
    vin = ufloat(1.01, 0.01)
    vout = unumpy.uarray([0.640, 3.02, 5.27, 9.2, 9.6, 6.4], [0.01, 0.01, 0.01, 0.1, 0.1, 0.1])
    fase = unumpy.uarray([92, 108, 123, 166, 178, -125], [1, 2, 1, 1, 1, 1])

    Gs = vout / vin
    dB = 10 * unumpy.log10(Gs)

    if not g:
        return None

    f = plt.figure(figsize=(8, 8))
    f.suptitle("Differenziatore", fontsize=15, y=0.98)
    
    ax = f.add_subplot(211)

    ax.errorbar(x=unumpy.nominal_values(freq),
        y=unumpy.nominal_values(Gs),
        c='black', fmt='o-')

    ax.set_xlabel('Frequenza', fontsize=14)
    ax.set_ylabel('Guadagno G', fontsize=14)

    ax.set_xscale('log')
    #ax.set_ylim((-13, 1))
    #ax.set_yticklabels(('', 2, 4, 6, 8, 10, 12))
    ax.set_xticklabels(('', '100 Hz', u"1 kHz", u"10 kHz", u"100 kHz", u"1 MHz"))
    
    ax.minorticks_on()
    ax.grid(b=True, which='major', color='0.7', linestyle='-', zorder=-5)
    ax.grid(b=True, which='minor', color='0.9', linestyle='-', zorder=-9)
    ax.set_axisbelow(True)
    
    ax2 = f.add_subplot(212)

    ax2.errorbar(x=unumpy.nominal_values(freq),
        y=unumpy.nominal_values(fase),
        c='black', fmt='o-')

    ax2.set_ylabel('Sfasamento [Gradi]', fontsize=14)

    ax2.set_xscale('log')
    #ax2.set_yticklabels(('', 25, 50, 75, 100))
    ax2.set_xticklabels(('', '100 Hz', u"1 kHz", u"10 kHz", u"100 kHz"))
    
    ax2.minorticks_on()
    ax2.grid(b=True, which='major', color='0.7', linestyle='-', zorder=-5)
    ax2.grid(b=True, which='minor', color='0.9', linestyle='-', zorder=-9)
    ax2.set_axisbelow(True)

    ax3 = ax2.twiny()
    ax3.set_xticks((0, 0.333, 0.666, 1))
    ax3.set_xticklabels(('', "1 kHz", u"10 kHz", u"100 kHz"))

    f.subplots_adjust(top=0.93, hspace=0.25, bottom=0.07, right=0.95)

    plt.savefig("../latex/diff.pdf")

    plt.show()

def test_x2_in_x1_2():
    """
    x2 has a couple of bins, each of which span more than one original bin
    """
    # old size
    m = 10

    # bin edges
    x_old = np.linspace(0., 1., m+1)
    x_new = np.array([0.25, 0.55, 0.75])

    # some arbitrary distribution
    y_old = 1. + np.sin(x_old[:-1]*np.pi) / np.ediff1d(x_old)

    y_old = unp.uarray(y_old, 0.1*y_old*uniform((m,)))

    # rebin
    y_new = rebin.rebin(x_old, y_old, x_new, interp_kind='piecewise_constant')

    # compute answer here to check rebin
    y_new_here = unp.uarray(np.zeros(2), np.zeros(2))
    y_new_here[0] = 0.5 * y_old[2] + y_old[3] + y_old[4] + 0.5 * y_old[5]
    y_new_here[1] = 0.5 * y_old[5] + y_old[6] + 0.5 * y_old[7]

    assert_allclose(unp.nominal_values(y_new),
                   unp.nominal_values(y_new_here))

    # mean or nominal value comparison
    assert_allclose(unp.std_devs(y_new),
                       unp.std_devs(y_new_here))

def detcoefs(f,nu,re,fh,fp,pr=None):
    """
    functia derermina coeficientii fit-ului functiei
    
    Nu=a*Re^b*Fh^c sau a functiei
    cf=a*Re^b*Fh^c
    
    Parametrii
    ----------
    f  : functie cu proptotipul :
        f((x,y),*coefs)
    nu : numpy.array of ufloats
        valorile nusselt pentru toate masuratorile
    re : numpy.array of ufloats
        valorile reynolds pentru toate masuratorile
    fh : numpy.array of floats
        valorile raportului h/Dh pentru toate masuratorile
    fp : numpy.array of floats
        valorile raportului p/Dh pentru toate masuratorile
        
    Intoarce
    --------
    coefs : array of ufloat 
        coeficientii functiei fitate cu erorile lor
    rchi2 : float 
        chi redus pentru analiza corectitudinii fit-ului
    dof : integer
        gradele de libertate

    """
    nu_n=uns.nominal_values(nu)
    nu_s=uns.std_devs(nu)
    w_nu=nu_s/nu_n
    re_n=uns.nominal_values(re)
    pr_n=uns.nominal_values(pr)
    if pr!=None:
        popt,pcov=curve_fit(f,(re_n,pr_n,fh,fp),nu_n,sigma=w_nu,
                        maxfev=1500)
        chi2=sum(((f((re_n,pr_n,fh,fp),*popt)-nu_n)/nu_s)**2)
    else:
        popt,pcov=curve_fit(f,(re_n,fh,fp),nu_n,sigma=w_nu,
                            maxfev=1500)
        chi2=sum(((f((re_n,fh,fp),*popt)-nu_n)/nu_s)**2)
    
    dof=len(nu_n)-len(popt)
    rchi2=chi2/dof
    
    coefs=[]
    for i in range(len(popt)):
        coefs.append(un.ufloat(popt[i],np.sqrt(pcov[i,i])))
    if pr!=None:
        func=lambda x,y,z,k:f((x,y,z,k),*popt)
    else:
        func=lambda x,y,z:f((x,y,z),*popt)
        
    return {'coefs':np.array(coefs),
            'rchi2':rchi2,
            'DOF':dof,
            'f':func
            }

def auswerten(name, d, n, t, z, V_mol, eps, raw):
    d *= 1e-3
    N = unp.uarray(n/t, np.sqrt(n)/t) - N_u

    if name=="Cu":
        tools.table((raw[0], raw[1], N), ("D/mm", "n", "(N-N_U)/\per\second"), "build/{}.tex".format(name), "Messdaten von {}.".format(name), "tab:daten{}".format(name), split=2, footer=r"$\Delta t = \SI{60}{s}$")#"(N-N_U)/\per\second"
    else:
        tools.table((raw[0], raw[1], raw[2], N), ("D/mm", "n", "\Delta t/s", "(N-N_U)/\per\second"), "build/{}.tex".format(name), "Messdaten von {}.".format(name), "tab:daten{}".format(name), split=2)
    mu = z * const.N_A / V_mol * 2 * np.pi * (const.e**2 / (4 * np.pi * const.epsilon_0 * const.m_e * const.c**2))**2 * ((1+eps)/eps**2 * ((2 * (1+eps))/(1+2*eps) - 1/eps * np.log(1+2*eps)) + 1/(2*eps) * np.log(1+ 2*eps) - (1+ 3*eps)/(1+2*eps)**2)

    params, pcov = curve_fit(fit, d, unp.nominal_values(N), sigma=unp.std_devs(N))
    params_ = unc.correlated_values(params, pcov)
    print("{}: N(0) = {}, µ = {}, µ_com = {}".format(name, params_[0], -params_[1], mu))

    sd = np.linspace(0, .07, 1000)

    valuesp = (fit(sd, *(unp.nominal_values(params_) + 10*unp.std_devs(params_)))).astype(float)
    valuesm = (fit(sd, *(unp.nominal_values(params_) - 10*unp.std_devs(params_)))).astype(float)

    #plt.xlim(0,7)
    plt.xlabel(r"$D/\si{mm}$")
    plt.ylabel(r"$(N-N_U)/\si{\per\second}$")
    plt.plot(1e3*sd, fit(sd, *params), 'b-', label="Fit")
    plt.fill_between(1e3*sd, valuesm, valuesp, facecolor='blue', alpha=0.125, edgecolor='none', label=r'$1\sigma$-Umgebung ($\times 10$)')
    plt.errorbar(1e3*d, unp.nominal_values(N), yerr=unp.std_devs(N), fmt='rx', label="Messdaten")
    plt.legend(loc='best')
    plt.yscale('linear')
    plt.tight_layout(pad=0)
    plt.savefig("build/{}.pdf".format(name))
    plt.yscale('log')
    plt.savefig("build/{}_log.pdf".format(name))
    plt.clf()

def stability(ux, ug_crit, omega, label=None):
    formatter = FuncFormatter(lambda x, pos: "%g\u00B2" % np.sqrt(x))

    ux, ug_crit = _double_sort(ux, ug_crit)

    x = np.power(ux, 2)
    y = ug_crit

    line = plt.errorbar(
        unp.nominal_values(x), unp.nominal_values(y), xerr=unp.std_devs(x), yerr=unp.std_devs(y), fmt="o", label=label
    )
    color = line.lines[0].get_color()

    plt.gca().xaxis.set_major_formatter(formatter)
    plt.gca().xaxis.set_ticks(np.power(np.arange(4, 10) * 100, 2))
    plt.xlabel(r"$U_x^2$")
    plt.ylabel(r"$U_g$")

    linear = lambda x, a: a * x
    a = _uarray_fit(linear, x, y, x0=(0.0001,), epsfcn=1e-7)[0]

    x = np.linspace(0, unp.nominal_values(x).max() * 1.1, 20)
    y = linear(x, a.n)

    plt.plot(x, linear(x, a.n), color=color)
    plt.fill_between(x, linear(x, a.n + a.s), linear(x, a.n - a.s), color=color, alpha=0.1)

    a = _normalize_ufloat(a)
    print("Slope:", a * 1000, "1/kV")

    qm = -2 / 3 * a * r0 ** 2 * omega ** 2 / K

    print(label, "q/m = {:.4P} uC/kg".format(qm * 1e6))

def import_data_from_objLog(FilesList, Objects_Include, pv):
    
    
    List_Abundances     = ['OI_HI', 'NI_HI', 'SI_HI', 'SI_HI_ArCorr', 'Y_Mass_O', 'Y_Mass_S', 'Y_Inference_O', 'Y_Inference_S']
    #List_Abundances    = ['OI_HI', 'NI_HI', 'SI_HI', 'SI_HI_ArCorr', 'Y_Mass_O', 'Y_Mass_S', 'Y_inf_O', 'Y_inf_S']

    #Dictionary of dictionaries to store object abundances
    Abund_dict = OrderedDict()
    for abund in List_Abundances:
        Abund_dict[abund] = OrderedDict()

    #Loop through files
    for i in range(len(FilesList)):
        #Analyze file address
        CodeName, FileName, FileFolder  = pv.Analyze_Address(FilesList[i])  
      
        if CodeName in Objects_Include:
            #Loop through abundances in the log
            for abund in List_Abundances:
                Abund_Mag = pv.GetParameter_ObjLog(CodeName, FileFolder, Parameter = abund, Assumption = 'float')
                #If the abundance was measure store it 
                if Abund_Mag != None:
                    Abund_dict[abund][CodeName] = Abund_Mag
        
    #Dictionary to store objects with abundances pairs for regressions. 
    #As an initial value for the keys we define the abundances we want to use for the regression
    Abundances_Pairs_dict = OrderedDict()
    Abundances_Pairs_dict['O_Regression']                   = ('OI_HI','Y_Mass_O')      
    Abundances_Pairs_dict['N_Regression']                   = ('NI_HI','Y_Mass_O')      
    Abundances_Pairs_dict['S_Regression']                   = ('SI_HI','Y_Mass_S')      
    Abundances_Pairs_dict['S_ArCorr_Regression']            = ('SI_HI_ArCorr','Y_Mass_S')
    Abundances_Pairs_dict['O_Regression_Inference']         = ('OI_HI','Y_Inference_O')      
    Abundances_Pairs_dict['N_Regression_Inference']         = ('NI_HI','Y_Inference_O')      
    Abundances_Pairs_dict['S_Regression_Inference']         = ('SI_HI','Y_Inference_S')      
    Abundances_Pairs_dict['S_ArCorr_Regression_Inference']  = ('SI_HI_ArCorr','Y_Inference_S') 
        
    #Loop through the regression lists and get objects with both abundances observed
    for j in range(len(Abundances_Pairs_dict)):
        
        #Get the elements keys for the regression
        Vector, Elem_X, Elem_Y = Abundances_Pairs_dict.keys()[j], Abundances_Pairs_dict.values()[j][0], Abundances_Pairs_dict.values()[j][1]
        
        #Determine objects with both abundances observed
        Obj_vector  = intersect1d(Abund_dict[Elem_X].keys(), Abund_dict[Elem_Y].keys(), assume_unique = True)
        X_vector    = zeros(len(Obj_vector))
        Y_vector    = zeros(len(Obj_vector))
        X_vector_E  = zeros(len(Obj_vector))
        Y_vector_E  = zeros(len(Obj_vector))
                        
        #Generate abundances vectors
        for z in range(len(Obj_vector)):  
            X_vector[z] = nominal_values(Abund_dict[Elem_X][Obj_vector[z]])
            X_vector_E[z] = std_devs(Abund_dict[Elem_X][Obj_vector[z]])            
            Y_vector[z] = nominal_values(Abund_dict[Elem_Y][Obj_vector[z]])
            Y_vector_E[z] = std_devs(Abund_dict[Elem_Y][Obj_vector[z]])
    
        Abundances_Pairs_dict[Vector] = (list(Obj_vector), uarray(X_vector, X_vector_E), uarray(Y_vector, Y_vector_E))
        
    return Abundances_Pairs_dict

def test_nom_ueig():
    sA = array([[1, 2], [3, 4]])
    A = array([[0.1, 0.2], [0.1, 0.3]])
    w, v = eig(A)
    uA = uarray((A, sA))
    uw, uv = ueig(A)
    assert nominal_values(uw) == w
    assert nominal_values(uv) == v

def fitting_powerlaw_LP(lx_min, lx_max,
        const=[ufloat(0.083,0.058), ufloat(2.11,0.21)]):
    x = lx_range(lx_min, lx_max)
    y = (10**24.5) * ((x/1e45)**const[1]) * (10**const[0])
    y_nom = unumpy.nominal_values(y)
    y_min = unumpy.nominal_values(y) - unumpy.std_devs(y)
    y_max = unumpy.nominal_values(y) + unumpy.std_devs(y)
    return y_nom, y_min, y_max

def fitting_powerlaw_YP(sz_min, sz_max,
        const=[ufloat(-0.133,0.069), ufloat(2.03,0.30)]):
    x = lx_range(sz_min, sz_max)
    y = (10**24.5) * ((x/1e-4)**const[1]) * (10**const[0])
    y_nom = unumpy.nominal_values(y)
    y_min = unumpy.nominal_values(y) - unumpy.std_devs(y)
    y_max = unumpy.nominal_values(y) + unumpy.std_devs(y)
    return y_nom, y_min, y_max

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 ulinReg(xdata,ydata):
    """
    Führt über linReg_iter eine lineare Regression durch. Parameter sind hierbei
    allerdings uncertainties.ufloat
    """
    popt, perr = linReg_iter(unumpy.nominal_values(xdata),
                 unumpy.nominal_values(ydata), unumpy.std_devs(ydata),
                 unumpy.std_devs(xdata))
    return uc.ufloat(popt[0],perr[0]), uc.ufloat(popt[1],perr[1])

    def state_space(self, speed, nominal=False):
        """
        Returns the A and B matrices for the Whipple model linearized about
        the upright constant velocity configuration.


        Parameters
        ----------
        speed : float
            The speed of the bicycle.
        nominal : boolean, optional
            The default is false and uarrays are returned with the calculated
            uncertainties. If true ndarrays are returned without uncertainties.

        Returns
        -------

        A : ndarray, shape(4,4)
            The state matrix.
        B : ndarray, shape(4,2)
            The input matrix.

        Notes
        -----
        ``A`` and ``B`` describe the Whipple model in state space form:

            x' = A * x + B * u

        where

        The states are [roll angle,
                        steer angle,
                        roll rate,
                        steer rate]

        The inputs are [roll torque,
                        steer torque]

        If you have a flywheel defined, body D, it will completely be ignored
        in these results. These results are strictly for the Whipple bicycle
        model.

        """

        M, C1, K0, K2 = self.canonical()

        g = self.parameters['Benchmark']['g']

        A, B = bicycle.ab_matrix(M, C1, K0, K2, speed, g)

        if nominal is True:
            return (unumpy.nominal_values(A), unumpy.nominal_values(B))
        elif nominal is False:
            return A, B
        else:
            raise ValueError('nominal must be True or False')

def plot(x,y,*args,**kwargs):
    nominal_curve = pyplot.plot(x, unumpy.nominal_values(y), *args, **kwargs)
    pyplot.fill_between(x, 
                        unumpy.nominal_values(y)-unumpy.std_devs(y), 
                        unumpy.nominal_values(y)+unumpy.std_devs(y),
                        facecolor=nominal_curve[0].get_color(),
                        edgecolor='face',
                        alpha=0.1,
                        linewidth=0)
    return nominal_curve

def GenExpErrorPlot(ras,figoffset,save=True,close=True,**figaspect): 
    colors=['b','g','r','c','m','y','k']
    markers=['.','*','o','v','^','<','>']
    index=0  
    Nufig=figoffset
    Ffig=Nufig+1
    for ra in ras:
        fh=Fh(ra)
        plt.figure(figoffset)
        plt.errorbar(uns.nominal_values(ra.Re.magnitude),
                 uns.nominal_values(ra.Nu.magnitude),
                 uns.std_devs(ra.Nu.magnitude),
                 uns.std_devs(ra.Re.magnitude),
                 '{:s}{:s}'.format(colors[index%len(colors)],
                                markers[index%len(markers)]),
                 label='{:.3f}'.format(fh))
        plt.figure(Ffig)
        plt.errorbar(uns.nominal_values(ra.Re.magnitude),
                 uns.nominal_values(ra.fr.magnitude),
                 uns.std_devs(ra.Re.magnitude),
                 uns.std_devs(ra.fr.magnitude),
                 '{:s}{:s}'.format(colors[index%len(colors)],
                                markers[index%len(markers)]),
                 label='{:.3f}'.format(fh))
    
    plt.figure(Nufig)
    plt.legend(loc='upper left')
    plt.xlabel('Re')
    plt.ylabel('Nu')
    
    plt.figure(Ffig)
    plt.legend(loc='upper right')
    plt.xlabel('Re')
    plt.ylabel('$c_f$')

    if save:
        plt.figure(Nufig)
        plt.savefig('Nu_exp_all.png',dpi=300,
                    figsize=(166.54/2.54,81/2.54),
                    orientation='landscape',
                    facecolor='w',
                    edgecolor='k')
        plt.figure(Ffig)
        plt.savefig('cf_exp_all.png',dpi=300,
                    figsize=(166.54/2.54,81/2.54),
                    orientation='landscape',
                    facecolor='w',
                    edgecolor='k')
    if close:
        plt.close(Nufig)
        plt.close(Ffig)
        return figoffset,-1,-1
    
    return Ffig+1,Nufig,Ffig

def compute_hwz(N_list, ttor, fit, plotname, title, sl=slice(None,None), Uscale=1, p0=None, eq=None, plabels=None, punits=None, Th_erw=None):
    
    N = np.sum(unp.uarray(N_list,np.sqrt(N_list)), axis=0)
    t = np.arange(len(N))*ttor+ttor/2.

    table = pt.PrettyTable()
    table.add_column('t [s]', t.astype(int), align='r')
    if len(N_list) > 1:
        for i in range(len(N_list)):
            table.add_column('N'+str(i+1), N_list[i].astype(int), align='r')
        table.add_column('Summe', N, align='r')
    else:
        table.add_column('N', N, align='r')
    with open("Resources/table_"+plotname+".txt", "w") as text_file:
        text_file.write(table.get_string())


    global N_U
    N_U = N_U0*Uscale*ttor
    popt, pstats = papstats.curve_fit(fit, t[sl], N[sl], p0=p0)

    # Untergrundfehler
    N_U = (N_U0-N_U0.s)*Uscale*ttor
    popt_min, pstats_min = papstats.curve_fit(fit, t[sl], N[sl], p0=p0)
    N_U = (N_U0+N_U0.s)*Uscale*ttor
    popt_max, pstats_max = papstats.curve_fit(fit, t[sl], N[sl], p0=p0)
    N_U = N_U0*Uscale*ttor
    s_U = unp.nominal_values(((np.abs(popt-popt_min)+np.abs(popt-popt_max))/2.))
    s_corrected = np.sqrt(unp.std_devs(popt)**2 + s_U**2)
    popt_corrected = unp.uarray(unp.nominal_values(popt),s_corrected)
    
    # Halbwertszeit
    Th = popt_corrected[::2]*unc.umath.log(2)
    for i in range(len(Th)):
        papstats.print_rdiff(Th[i]/60, Th_erw[i]/60)

    # Plot
    plt.clf()
    plt.title('Diagramm '+plotname+': '+title)
    plt.xlabel('Messzeit $t \, [s]$')
    plt.ylabel('Ereigniszahl $N$')
    xspace = np.linspace(0, t[-1])
    papstats.plot_data(t, N, label='Messpunkte')
    papstats.plot_fit(fit, popt, pstats, xspace, eq=eq, plabels=plabels, punits=punits)
    plt.fill_between(xspace, fit(xspace, *unp.nominal_values(popt_min)), fit(xspace, *unp.nominal_values(popt_max)), color='g', alpha=0.2)
    Nmin = np.amin(unp.nominal_values(N))
    for i in range(len(Th)):
        plt.hlines(popt[1::2][i].n/2.+N_U.n, 0, Th[i].n, lw=2, label='Halbwertszeit $'+papstats.pformat(Th[i], label=r'T_{\frac{1}{2}}'+('^'+str(i+1) if len(Th)>1 else ''), unit='s')+'$')
    handles, labels = plt.gca().get_legend_handles_labels()
    p = plt.Rectangle((0, 0), 1, 1, color='g', alpha=0.2)
    handles.append(p)
    labels.append('Fit im '+r'$1 \sigma$'+'-Bereich von $N_U$:'+''.join(['\n$'+papstats.pformat(s_U[i], label='\Delta '+plabels[i]+'^{U}', unit=punits[i])+'$' for i in range(len(plabels))]))
    plt.legend(handles, labels)
    papstats.savefig_a4(plotname+'.png')

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 residuals(self):
		'''
		usually called by draw()
		calculates residuals
		'''
		from uncertainties import unumpy
		from numpy import array
		from ROOT import TGraphErrors
		residuals = array( range( len( self.__x ) ) ,'float')
		for i in range( len( unumpy.nominal_values(self.__x) ) ):
			residuals[i] = unumpy.nominal_values(self.__y)[i] - self.func.Eval( unumpy.nominal_values(self.__x)[i] )
		self.resgraph = TGraphErrors( len( self.__x ), unumpy.nominal_values(self.__x), residuals, unumpy.std_devs(self.__x), unumpy.std_devs(self.__y) )

def PlotPatches(Sc,PatchData,ErrorBars):
    """
    Plot E*R* data binned from hilltop patches.
    """
    e_star = E_Star(Sc,PatchData[2],PatchData[0])
    r_star = R_Star(Sc,PatchData[1],PatchData[0])
    if ErrorBars:
        plt.errorbar(unp.nominal_values(e_star),unp.nominal_values(r_star),yerr=unp.std_devs(r_star),xerr=unp.std_devs(e_star),
                     fmt='ro',label='Hilltop Patch Data')
    else:
        plt.errorbar(unp.nominal_values(e_star),unp.nominal_values(r_star),
                     fmt='ro',label='Hilltop Patch Data')

def curve_fit(fit, xdata, ydata, sigma=None, p0=None):
    if sigma is None:
        # TODO: Nur y-Fehler? chisquared vergleicht nur y-Differenzen, x-Fehler relevant für Fit?
        sigma = unp.std_devs(ydata)
        if np.sum(sigma) == 0:
            sigma = None
    xdata = unp.nominal_values(xdata)
    ydata = unp.nominal_values(ydata)
    popt, pcov = opt.curve_fit(fit, xdata, ydata, sigma=sigma, p0=p0)
    popt = unp.uarray(popt, np.sqrt(np.diagonal(pcov)))
    pstats = PAPStats(ydata, fit(xdata, *unp.nominal_values(popt)), sigma=sigma, ddof=len(popt))
    return popt, pstats

def plot_data(xdata, ydata, ax=plt, **kwargs):
    xerr = unp.std_devs(xdata)
    if np.sum(xerr)==0:
        xerr = None
    yerr = unp.std_devs(ydata)
    if np.sum(yerr)==0:
        yerr = None
    if not (kwargs.has_key('ls') or kwargs.has_key('linestyle')):
        kwargs['ls'] = 'none'
    if not kwargs.has_key('marker'):
        kwargs['marker'] = '.'
    return ax.errorbar(unp.nominal_values(xdata), unp.nominal_values(ydata), xerr=xerr, yerr=yerr, **kwargs)