generate uniform grid of points on Ra-Dec plane around the given target (Ra,Dec)

Question

want to form a uniform square grid in the RA-Dec plane around the given target position. Target may be anywhere i.e. may be near the pol or may near equator but the grid should be uniform.

inputs are :

1. target position: Ra and Dec. ( hours and Deg)
2. grid size: N X M ( N is cols, M is rows)
3. grid spacing: d (in deg)

output:

list of Ra-Dec points

Please provide formulas or python script to generate these points.

here is the python code but it is not giving expected answer

import numpy as np
from astropy.coordinates import SkyCoord
import astropy.units as u

def generate_uniform_grid(target_ra, target_dec, n_cols, n_rows, spacing_deg):
    """
    Generates a uniform RA-Dec grid around a target RA-Dec position.
    This algorithm works for all positions, including poles, without special cases.

    Parameters:
        target_ra (float): Target RA in hours.
        target_dec (float): Target Dec in degrees.
        n_cols (int): Number of columns in the grid (RA offsets).
        n_rows (int): Number of rows in the grid (Dec offsets).
        spacing_deg (float): Spacing between adjacent grid points in degrees.

    Returns:
        list: List of (RA, Dec) tuples for grid points in degrees.
    """
    # Convert target position to Cartesian coordinates
    target_coord = SkyCoord(ra=target_ra * u.hourangle, dec=target_dec * u.deg)
    target_vector = target_coord.cartesian.xyz.value  # 3D Cartesian vector

    # Grid parameters
    half_width = (n_cols - 1) / 2
    half_height = (n_rows - 1) / 2
    dec_offsets = np.linspace(-half_height * spacing_deg, half_height * spacing_deg, n_rows)
    ra_offsets = np.linspace(-half_width * spacing_deg, half_width * spacing_deg, n_cols)

    # Create local grid offsets in Cartesian coordinates
    grid_vectors = []
    for dec_offset in dec_offsets:
        for ra_offset in ra_offsets:
            # Offset spherical coordinates
            offset_coord = SkyCoord(
                ra=ra_offset * u.deg,
                dec=dec_offset * u.deg,
                frame="icrs",
            )
            offset_vector = offset_coord.cartesian.xyz.value

            # Rotate the offset vector to align with the target position
            rotated_vector = _rotate_to_target(offset_vector, target_vector)
            grid_vectors.append(rotated_vector)

    # Convert all grid vectors back to RA-Dec
    grid_points = []
    for vector in grid_vectors:
        point = SkyCoord(
            x=vector[0],
            y=vector[1],
            z=vector[2],
            representation_type="cartesian"
        ).spherical
        grid_points.append((point.lon.deg % 360, point.lat.deg))

    return grid_points


def _rotate_to_target(offset_vector, target_vector):
    """
    Rotates the offset vector to align with the target vector in 3D space.

    Parameters:
        offset_vector (np.ndarray): 3D Cartesian unit vector for the offset.
        target_vector (np.ndarray): 3D Cartesian unit vector for the target.

    Returns:
        np.ndarray: Rotated 3D Cartesian vector.
    """
    # Cross product for the rotation axis
    rotation_axis = np.cross([0, 0, 1], target_vector)
    if np.linalg.norm(rotation_axis) < 1e-10:  # Handle pole case
        return offset_vector

    rotation_axis /= np.linalg.norm(rotation_axis)  # Normalize

    # Angle between vectors
    angle = np.arccos(np.dot([0, 0, 1], target_vector))

    # Rotation matrix
    rotation_matrix = _rotation_matrix(rotation_axis, angle)

    # Apply the rotation
    rotated_vector = np.dot(rotation_matrix, offset_vector)
    return rotated_vector


def _rotation_matrix(axis, angle):
    """
    Constructs a 3D rotation matrix given an axis and angle.

    Parameters:
        axis (np.ndarray): Rotation axis (3D unit vector).
        angle (float): Rotation angle in radians.

    Returns:
        np.ndarray: 3x3 rotation matrix.
    """
    x, y, z = axis
    c = np.cos(angle)
    s = np.sin(angle)
    t = 1 - c

    return np.array([
        [t * x * x + c, t * x * y - s * z, t * x * z + s * y],
        [t * x * y + s * z, t * y * y + c, t * y * z - s * x],
        [t * x * z - s * y, t * y * z + s * x, t * z * z + c],
    ])


# Example Usage
target_ra = 0.0  # Target RA in hours
target_dec = 90.0  # Target Dec in degrees (north pole)
n_cols = 3  # Grid columns
n_rows = 3  # Grid rows
spacing_deg = 45.0  # Spacing between grid points in degrees

grid = generate_uniform_grid(target_ra, target_dec, n_cols, n_rows, spacing_deg)
for ra, dec in grid:
    print(f"RA: {ra:.6f} deg, Dec: {dec:.6f} deg")

The expected answer is

RA: 315.000000 deg, Dec: 45.000000 deg

RA: 0.000000 deg, Dec: 45.000000 deg

RA: 45.000000 deg, Dec: 45.000000 deg

RA: 315.000000 deg, Dec: 90.000000 deg

RA: 0.000000 deg, Dec: 90.000000 deg

RA: 45.000000 deg, Dec: 90.000000 deg

RA: 315.000000 deg, Dec: 45.000000 deg

RA: 0.000000 deg, Dec: 45.000000 deg

RA: 45.000000 deg, Dec: 45.000000 deg