def fit_Zr_L(sZr, ll=None):
    "Fits the Zr L lines, return a hyperspy model. Input is a Zr spectrum image, and optionally, the accompanying low loss SI."
    print("Will produce 10 progress bars")

    mZr = sZr.create_model(ll = ll, GOS="Hartree-Slater", auto_add_edges=False)
    # Fit Background
    mZr.fit_component(mZr["PowerLaw"], bounded=True, signal_range=[2150.,2210.], fit_independent=True, only_current=True)
    mZr.assign_current_values_to_all()
    mZr.fit_component(mZr["PowerLaw"], bounded=True, signal_range=[2150.,2210.], fit_independent=True, only_current=False)

    Zr_L3 = hs.model.components1D.EELSCLEdge("Zr_L3", GOS="Hartree-Slater")
    Zr_L2 = hs.model.components1D.EELSCLEdge("Zr_L2", GOS="Hartree-Slater")
 # Fit Zr L3 white line first, use centre position of gaussian to choose onset position of core loss edge, then repeat
    Zr_L3_white = hs.model.components1D.Gaussian()
    Zr_L3_white.name = "Zr_L3 White Line"
    Zr_L3_white.centre.bmin = 2218.0
    Zr_L3_white.centre.bmax = 2240.0
    Zr_L3_white.A.bmin = 0
    Zr_L3_white.sigma.bmax=3
    mZr.append(Zr_L3_white)

    mZr.fit_component(Zr_L3_white,fitter="leastsq", signal_range=[Zr_L3_white.centre.bmin-3.0, Zr_L3_white.centre.bmax+3.0], bounded=True,)
    mZr.assign_current_values_to_all([Zr_L3_white])
    mZr.fit_component(Zr_L3_white,fitter="leastsq", signal_range=[Zr_L3_white.centre.bmin-3.0, Zr_L3_white.centre.bmax+3.0], bounded=True, only_current=False)

    mZr.append(Zr_L3)
    Zr_L3.onset_energy.twin = Zr_L3_white.centre
    mZr.fit_component(Zr_L3, bounded=True, signal_range=[2240., 2290.])
    mZr.assign_current_values_to_all([Zr_L3])
    mZr.fit_component(Zr_L3, bounded=True, signal_range=[2240., 2290.], only_current=False)

# Fit Zr L2 white line, use centre position of gaussian to choose onset position of core loss edge, then repeat

    Zr_L2_white = hs.model.components1D.Gaussian()
    Zr_L2_white.name = "Zr_L2 White Line"
    Zr_L2_white.centre.bmin = 2309.0
    Zr_L2_white.centre.bmax = 2317.0
    Zr_L2_white.A.bmin = 0
    Zr_L2_white.sigma.bmax=3

    mZr.append(Zr_L2_white)
    mZr.fit_component(Zr_L2_white, signal_range=[Zr_L2_white.centre.bmin-3.0, Zr_L2_white.centre.bmax+3.0], bounded=True,)

    mZr.assign_current_values_to_all([Zr_L2_white])
    mZr.fit_component(Zr_L2_white, signal_range=[Zr_L2_white.centre.bmin-3.0, Zr_L2_white.centre.bmax+3.0], bounded=True, only_current=False)

    mZr.append(Zr_L2)
    Zr_L2.onset_energy.twin = Zr_L2_white.centre
    mZr.fit_component(Zr_L2, bounded=True, signal_range=[2335., 2392.])
    mZr.assign_current_values_to_all([Zr_L2])
    mZr.fit_component(Zr_L2, bounded=True, signal_range=[2335., 2392.], only_current=False)
    mZr.fit_component(Zr_L2_white, signal_range=[Zr_L2_white.centre.bmin-3.0, Zr_L2_white.centre.bmax+3.0], bounded=True, only_current=False)

# Finally cycle over components one more time to ensure that fit is good
    mZr.fit_component(Zr_L3, bounded=True, signal_range=[2240., 2290.], only_current=False)
    mZr.fit_component(Zr_L3_white, signal_range=[Zr_L3_white.centre.bmin-3.0, Zr_L3_white.centre.bmax+3.0], bounded=True, only_current=False)
    mZr.fit_component(Zr_L2, bounded=True, signal_range=[2335., 2392.], only_current=False)
    mZr.fit_component(Zr_L2_white, signal_range=[Zr_L2_white.centre.bmin-3.0, Zr_L2_white.centre.bmax+3.0], bounded=True, only_current=False)
    print("Finished fitting Zirconium")
    return mZr