def METHOD_NAME(label="Volume", nphi=None, ntheta=None):
    """
    Returns a boozer surface that will be used in unit tests.
    """
    assert label == "Volume" or label == "ToroidalFlux" or label == "Area"
    base_curves, base_currents, ma = get_ncsx_data()
    coils = coils_via_symmetries(base_curves, base_currents, 3, True)
    bs = BiotSavart(coils)
    current_sum = sum(abs(c.current.get_value()) for c in coils)
    G0 = 2. * np.pi * current_sum * (4 * np.pi * 10**(-7) / (2 * np.pi))
    ## RESOLUTION DETAILS OF SURFACE ON WHICH WE OPTIMIZE FOR QA
    mpol = 6
    ntor = 6
    stellsym = True
    nfp = 3
    phis = np.linspace(0, 1/nfp, 2*ntor+1, endpoint=False)
    thetas = np.linspace(0, 1, 2*mpol+1, endpoint=False)
    s = SurfaceXYZTensorFourier(
        mpol=mpol, ntor=ntor, stellsym=stellsym, nfp=nfp, quadpoints_phi=phis, quadpoints_theta=thetas)
    s.fit_to_curve(ma, 0.1, flip_theta=True)
    iota = -0.406
    if label == "Volume":
        lab = Volume(s, nphi=nphi, ntheta=ntheta)
        lab_target = lab.J()
    elif label == "ToroidalFlux":
        bs_tf = BiotSavart(coils)
        lab = ToroidalFlux(s, bs_tf, nphi=nphi, ntheta=ntheta)
        lab_target = lab.J()
    elif label == "Area":
        lab = Area(s, nphi=nphi, ntheta=ntheta)
        lab_target = lab.J()
    ## COMPUTE THE SURFACE
    boozer_surface = BoozerSurface(bs, s, lab, lab_target)
    res = boozer_surface.solve_residual_equation_exactly_newton(tol=1e-13, maxiter=20, iota=iota, G=G0)
    print(f"NEWTON {res['success']}: iter={res['iter']}, iota={res['iota']:.3f}, vol={s.volume():.3f}")
    return bs, boozer_surface