def METHOD_NAME(ps_full, N, L, S, check=False):
    # generate rs_index
    rs_index = []
    for alpha in range(0, S):
        for a in range(0, N):
            rs_index.append([a, alpha])
    # initiialise power spec
    d = np.zeros(shape=(S, S, L + 1, N, N))
    # assign entries
    ipow = 0
    for ia in range(0, N * S):
        a, i_species = rs_index[ia]
        for jb in range(0, ia + 1):
            b, j_species = rs_index[jb]
            for ll in range(0, L + 1):
                val = ps_full[ipow]
                if ia != jb:
                    val /= 2**0.5
                    d[i_species, j_species, ll, a, b] = val
                    d[j_species, i_species, ll, b, a] = val
                else:
                    d[i_species, j_species, ll, a, b] = val
                ipow += 1
    if check:
        d_sort = sorted(d.flatten())
        fp_sort = sorted(ps_full)
        assert np.allclose(d_sort, fp_sort)
    return d