def generate_em_constraint_data(mNS_min, mNS_max, delta_mNS, sBH_min, sBH_max, delta_sBH, eos_name, threshold, eta_default):
    """
    Wrapper that calls find_em_constraint_data_point over a grid
    of points to generate the bh_spin_z x ns_g_mass x eta surface
    above which NS-BH binaries yield a remnant disk mass that
    exceeds the threshold required by the user.  The user must also
    specify the default symmetric mass ratio value to be assigned
    to points for which the NS mass exceeds the maximum NS mass
    allowed by the chosend NS equation of state.
    The 2D surface that is generated is saved to file in two formats:
    constraint_em_bright.npz and constraint_em_bright.npz.

    Parameters
    -----------
    mNS_min: float
        lower boundary of the grid in the NS mass direction
    mNS_max: float
        upper boundary of the grid in the NS mass direction
    delta_mNS: float
        grid spacing in the NS mass direction
    sBH_min: float
        lower boundary of the grid in the direction of the
        BH dimensionless spin component along the orbital
        angular momentum
    sBH_max: float
        upper boundary of the grid in the direction of the
        BH dimensionless spin component along the orbital
        angular momentum
    delta_sBH: float
        grid spacing in the direction of the BH dimensionless
        spin component along the orbital angular momentum
    eos_name: string
        NS equation of state label ('2H' is the only supported
        choice at the moment)
    threshold: float
        an amount to be subtracted to the remnant mass upper limit
        predicted by the model (in solar masses)
    eta_default: float
        the value to be returned for points in the grids in which
        the NS mass is too high
    """
    # Build a grid of points in the mNS x sBHz space,
    # making sure maxima and minima are included
    mNS_nsamples = complex(0,int(np.ceil((mNS_max-mNS_min)/delta_mNS)+1))
    sBH_nsamples = complex(0,int(np.ceil((sBH_max-sBH_min)/delta_sBH)+1))
    mNS_vec, sBH_vec = np.mgrid[mNS_min:mNS_max:mNS_nsamples, sBH_min:sBH_max:sBH_nsamples] # pylint:disable=invalid-slice-index
    mNS_locations = np.array(mNS_vec[:,0])
    sBH_locations = np.array(sBH_vec[0])
    mNS_sBH_grid = zip(mNS_vec.ravel(), sBH_vec.ravel())
    mNS_sBH_grid = np.array(mNS_sBH_grid)
    mNS_vec = np.array(mNS_sBH_grid[:,0])
    sBH_vec = np.array(mNS_sBH_grid[:,1])

    # Until a numpy v>=1.7 is available everywhere, we have to use a silly
    # vectorization of find_em_constraint_data_point and pass to it a bunch of
    # constant arguments as vectors with one entry repeated several times
    eos_name_vec=[eos_name for _ in range(len(mNS_vec))]
    eos_name_vec=np.array(eos_name_vec)
    threshold_vec=np.empty(len(mNS_vec))
    threshold_vec.fill(threshold)
    eta_default_vec=np.empty(len(mNS_vec))
    eta_default_vec.fill(eta_default)

    # Compute the minimum etas at all point in the mNS x sBHz grid
    eta_sol = find_em_constraint_data_points(mNS_vec, sBH_vec, eos_name_vec, threshold_vec, eta_default_vec)
    eta_sol = eta_sol.reshape(-1,len(sBH_locations))
    # Save the results
    np.savez('constraint_em_bright', mNS_pts=mNS_locations, sBH_pts=sBH_locations, eta_mins=eta_sol)
    # Cast the results in a format that is quick to plot from textfile
    constraint_data = zip(mNS_vec.ravel(), sBH_vec.ravel(), eta_sol.ravel())
    np.savetxt('constraint_em_bright.dat', constraint_data)