def METHOD_NAME(omega, shape, npml, dlf, dlb, dmin_pml=(True, True)):
    """Makes the 'S-matrices'. When dotted with derivative matrices, they add
    PML. If dmin_pml is set to False, PML will not be applied on the "bottom"
    side of the domain."""
    # strip out some information needed
    Nx, Ny = shape
    N = Nx * Ny
    nx_pml, ny_pml = npml
    # Create the sfactor in each direction and for 'f' and 'b'
    s_vector_x_f = create_sfactor("f", omega, dlf[0], Nx, nx_pml, dmin_pml[0])
    s_vector_x_b = create_sfactor("b", omega, dlb[0], Nx, nx_pml, dmin_pml[0])
    s_vector_y_f = create_sfactor("f", omega, dlf[1], Ny, ny_pml, dmin_pml[1])
    s_vector_y_b = create_sfactor("b", omega, dlb[1], Ny, ny_pml, dmin_pml[1])
    # Fill the 2d space with layers of appropriate s-factors
    sx_f_2d = np.zeros(shape, dtype=np.complex128)
    sx_b_2d = np.zeros(shape, dtype=np.complex128)
    sy_f_2d = np.zeros(shape, dtype=np.complex128)
    sy_b_2d = np.zeros(shape, dtype=np.complex128)
    # Insert the cross sections into the S-grids (could be done more elegantly)
    for i in range(Ny):
        sx_f_2d[:, i] = 1 / s_vector_x_f
        sx_b_2d[:, i] = 1 / s_vector_x_b
    for i in range(Nx):
        sy_f_2d[i, :] = 1 / s_vector_y_f
        sy_b_2d[i, :] = 1 / s_vector_y_b
    # Reshape the 2d s-factors into a 1D s-vecay
    sx_f_vec = sx_f_2d.flatten()
    sx_b_vec = sx_b_2d.flatten()
    sy_f_vec = sy_f_2d.flatten()
    sy_b_vec = sy_b_2d.flatten()
    # Construct the 1D total s-vector into a diagonal matrix
    sx_f = sp.spdiags(sx_f_vec, 0, N, N)
    sx_b = sp.spdiags(sx_b_vec, 0, N, N)
    sy_f = sp.spdiags(sy_f_vec, 0, N, N)
    sy_b = sp.spdiags(sy_b_vec, 0, N, N)
    return sx_f, sx_b, sy_f, sy_b