Генерировать равномерную сетку точек на плоскости Ra-Dec вокруг заданной цели (Ra,Dec) ⇐ Python

[Anonymous]

хотим сформировать равномерную квадратную сетку в плоскости RA-Dec вокруг заданной целевой позиции. Цель может находиться где угодно, например, вблизи полюса или вблизи экватора, но сетка должна быть однородной.
Входные данные:

• целевая позиция: Ra и Dec. (часы и градусы).
• размер сетки: N X M ( N — столбцы, M — строки)
• шаг сетки: d (в град)

выход:
список точек Ra-Dec
Пожалуйста, предоставьте формулы или скрипт Python для генерации этих точек.
вот код Python, но он не дает ожидаемого ответа

Код: Выделить всё

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")

Ожидаемый ответ:
RA: 315,000000 градусов, декабрь: 45,000000 градусов
RA: 0,000000 градусов , декабрь: 45.000000 град.
RA: 45.000000 град., декабрь: 45,000000 град
RA: 315,000000 град, декабрь: 90,000000 град
RA: 0,000000 град, декабрь: 90,000000 град
RA: 45.000000 град, декабрь: 90,000000 град
RA: 315,000000 град, декабрь: 45,000000 град
RA: 0,000000 град, декабрь: 45,000000 град
RA: 45.000000 град, декабрь: 45,000000 град

Подробнее здесь: https://stackoverflow.com/questions/792 ... get-ra-dec