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 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 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 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 ucurve_fit(f, x, y, **kwargs):
    if np.any(unp.std_devs(y) == 0):
        sigma = None
    else:
        sigma = unp.std_devs(y)

    popt, pcov = scipy.optimize.curve_fit(f, x, unp.nominal_values(y), sigma=sigma, **kwargs)

    return unc.uarray(popt, np.sqrt(np.diag(pcov)))

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 ucurve_fit(f, x, y, **kwargs):
    if np.any(unp.std_devs(y) == 0):
        sigma = None
    else:
        sigma = unp.std_devs(y)

    popt, pcov = scipy.optimize.curve_fit(f, x, unp.nominal_values(y), sigma=sigma, **kwargs)

    return unc.correlated_values(popt, pcov)

def CdTe_Am():
    detector_name = "CdTe"
    sample_name = "Am"
    suff = "_" + sample_name + "_" + detector_name 
    npy_file = npy_dir + "spectra" + suff + ".npy"
    n = np.load(npy_file)[:500]
    s_n = np.sqrt(n)
    x = np.arange(len(n))
    
    # only fit
    p0 = [250000, 250, 25]
    red = (x > 235) * (x < 270)
    n_red = n[red]
    x_red = x[red]
    s_n_red = s_n[red]
    fit, cov_fit = curve_fit(gauss, x_red, n_red, p0=p0, sigma=s_n_red)
    fit  = np.abs(fit)
    fit_corr = uc.correlated_values(fit, cov_fit)
    s_fit = un.std_devs(fit_corr)
    # chi square
    n_d = 4
    chi2 = np.sum(((gauss(x_red, *fit) - n_red) / s_n_red) ** 2 )
    chi2_test = chi2/n_d
    fit_r, s_fit_r = err_round(fit, cov_fit)
    fit_both = np.array([fit_r, s_fit_r])
    fit_both = np.reshape(fit_both.T, np.size(fit_both))
    As[2], mus[2], sigmas[2] = fit
    s_As[2], s_mus[2], s_sigmas[2] = un.std_devs(fit_corr)

    all = np.concatenate((fit_both, [chi2_test]), axis=0)
    def plot_two_gauss():
        fig1, ax1 = plt.subplots(1, 1)
        if not save_fig:
            fig1.suptitle("Detector: " + detector_name + "; sample: " + sample_name)
        plot1, = ax1.plot(x, n, '.', alpha=0.3)   # histo plot
        ax1.errorbar(x, n, yerr=s_n, fmt=',', alpha=0.99, c=plot1.get_color(), errorevery=10) # errors of t are not changed!
        ax1.plot(x_red, gauss(x_red, *fit))
        ax1.set_xlabel("channel")
        ax1.set_ylabel("counts")
        textstr = 'Results of fit:\n\
                \\begin{eqnarray*}\
                A     &=& (%.0f \pm %.0f) \\\\ \
                \mu   &=& (%.1f \pm %.1f) \\\\ \
                \sigma&=& (%.1f \pm %.1f) \\\\ \
                \chi^2 / n_d &=& %.1f\
                \end{eqnarray*}'%tuple(all)
        ax1.text(0.65, 0.95, textstr, transform=ax1.transAxes, va='top', bbox=props)
        if show_fig:
            fig1.show()
        if save_fig:
            file_name = "detector" + suff
            fig1.savefig(fig_dir + file_name + ".pdf")
            fig1.savefig(fig_dir + file_name + ".png")
        return 0
    plot_two_gauss()
    return fit, s_fit

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 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_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)

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 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 PlotBasins(Sc, BasinData, ErrorBars):
    """
    Plot basin average E*R* data.
    """
    e_star = E_Star(Sc, BasinData[2], BasinData[0])
    r_star = R_Star(Sc, BasinData[1], BasinData[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='go', label='Basin Data')
    else:
        plt.errorbar(unp.nominal_values(e_star), unp.nominal_values(r_star),
                     fmt='go', label='Basin Data')

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 test_component_extraction():
    "Extracting the nominal values and standard deviations from an array"

    arr = unumpy.uarray(([1, 2], [0.1, 0.2]))

    assert numpy.all(unumpy.nominal_values(arr) == [1, 2])
    assert numpy.all(unumpy.std_devs(arr) == [0.1, 0.2])

    # unumpy matrices, in addition, should have nominal_values that
    # are simply numpy matrices (not unumpy ones, because they have no
    # uncertainties):
    mat = unumpy.matrix(arr)
    assert numpy.all(unumpy.nominal_values(mat) == [1, 2])
    assert numpy.all(unumpy.std_devs(mat) == [0.1, 0.2])
    assert type(unumpy.nominal_values(mat)) == numpy.matrix

def _uarray_fit(func, x, y, x0, B=1000, **kwargs):
    def _derivative(func, x, *params, h=1e-6):
        return (
            -func(x + 2 * h, *params) + 8 * func(x + h, *params) - 8 * func(x - h, *params) + func(x - 2 * h, *params)
        ) / (12 * h)

    def _chi(params, func, xn, yn, xs, ys):
        difference = yn - func(xn, *params)
        error = np.sqrt(np.power(_derivative(func, xn, *params) * xs, 2) + np.power(ys, 2))
        chi = difference / error
        return chi

    xs = unp.std_devs(x)
    ys = unp.std_devs(y)

    means = []
    for i in range(B):
        indices = np.random.randint(0, len(x), size=len(x))

        shifts1 = np.random.normal(loc=0, scale=1, size=len(x))
        x_simulated = unp.nominal_values(x) + xs * shifts1

        shifts2 = np.random.normal(loc=0, scale=1, size=len(y))
        y_simulated = unp.nominal_values(y) + ys * shifts2

        popt, pcov, infodict, mesg, ier = leastsq(
            _chi, x0=tuple(x0), args=(func, x_simulated, y_simulated, xs, ys), full_output=True, **kwargs
        )
        if ier in (1, 2, 3, 4):
            means.append(popt)

    popt, pcov, infodict, mesg, ier = leastsq(
        _chi,
        x0=tuple(x0),
        args=(func, unp.nominal_values(x), unp.nominal_values(y), xs, ys),
        full_output=True,
        **kwargs
    )

    errors = np.std(means, axis=0)
    results = tuple(ufloat(a, b) for a, b in zip(popt, errors))

    chisqndof = np.power(
        _chi(popt, func, unp.nominal_values(x), unp.nominal_values(y), unp.std_devs(x), unp.std_devs(y)), 2
    ).sum() / (len(x) - len(x0))
    print("Chi^2/ndof =", chisqndof)

    return results