n the following code a projection of a cone onto a plane is created. The problem is that actual position of the cone base is affected by inclination. Fix it.

 # Define cone parameters
cone_top = SkyCoord('7h28m43.5s 69d11m24s', frame='icrs')  # Cone top (apex) as SkyCoord
cone_base_center = SkyCoord('7h28m38s 69d12m28s', frame='icrs')  # Cone base center as SkyCoord
cone_opening_angle = 20 * u.deg  # Opening angle of the cone
cone_orientation = 10 * u.deg  # In-plane orientation angle
cone_inclination = 50 * u.deg  # Tilt relative to the plane of the sky

# Load WCS for the image (assumed to be already defined)
wcs = WCS(HD(h_ha))

def radec_to_cartesian(ra_deg, dec_deg):
    """Convert RA and Dec in degrees to Cartesian coordinates."""
    ra_rad = np.deg2rad(ra_deg)
    dec_rad = np.deg2rad(dec_deg)
    x = np.cos(dec_rad) * np.cos(ra_rad)
    y = np.cos(dec_rad) * np.sin(ra_rad)
    z = np.sin(dec_rad)
    return np.array([x, y, z])

def cartesian_to_radec(x, y, z):
    """Convert Cartesian coordinates to RA and Dec in degrees."""
    hyp = np.hypot(x, y)
    ra_rad = np.arctan2(y, x)
    dec_rad = np.arctan2(z, hyp)
    ra_deg = np.rad2deg(ra_rad) % 360.0
    dec_deg = np.rad2deg(dec_rad)
    return ra_deg, dec_deg

def cone_base_vertices(cone_top, cone_base_center, cone_opening_angle, cone_inclination, num_points=100):
    """
    Generate vertices for the cone's base, projected onto the sky,
    taking into account the cone_inclination.
    """
    # Convert apex and base center to Cartesian coordinates
    apex_vector = radec_to_cartesian(cone_top.ra.deg, cone_top.dec.deg)
    base_center_vector = radec_to_cartesian(cone_base_center.ra.deg, cone_base_center.dec.deg)
    
    # Compute the axis vector from apex to base center
    axis_vector = base_center_vector - apex_vector
    distance = np.linalg.norm(axis_vector)
    axis_vector = axis_vector / distance
    
    # Compute current inclination angle between apex_vector and axis_vector
    current_inclination = np.arccos(np.dot(apex_vector, axis_vector))
    
    # Convert cone_inclination to radians
    cone_inclination_rad = cone_inclination.to(u.radian).value
    
    # Compute the angle difference
    delta_inclination = cone_inclination_rad - current_inclination
    
    # If there is a difference, rotate axis_vector around rotation_axis by delta_inclination
    if not np.isclose(delta_inclination, 0):
        # Compute rotation axis (perpendicular to both apex_vector and axis_vector)
        rotation_axis = np.cross(apex_vector, axis_vector)
        if np.linalg.norm(rotation_axis) == 0:
            # apex_vector and axis_vector are parallel, choose an arbitrary perpendicular axis
            rotation_axis = np.cross(apex_vector, np.array([1, 0, 0]))
        rotation_axis /= np.linalg.norm(rotation_axis)
        
        # Rotate axis_vector by delta_inclination using Rodrigues' rotation formula
        axis_vector_rotated = (axis_vector * np.cos(delta_inclination) +
                               np.sin(delta_inclination) * np.cross(rotation_axis, axis_vector) +
                               (1 - np.cos(delta_inclination)) * np.dot(rotation_axis, axis_vector) * rotation_axis)
        axis_vector = axis_vector_rotated / np.linalg.norm(axis_vector_rotated)
    
    # Now axis_vector has the correct inclination
    # Recompute base_center_vector using adjusted axis_vector
    base_center_vector = apex_vector + axis_vector * distance
    base_center_vector /= np.linalg.norm(base_center_vector)
    
    # Compute base radius
    cone_opening_angle_rad = cone_opening_angle.to(u.radian).value
    base_radius = np.tan(cone_opening_angle_rad) * distance
    
    # At base center, compute orthonormal vectors perpendicular to axis_vector
    # Find a vector not parallel to axis_vector
    if np.allclose(axis_vector, [0, 0, 1]) or np.allclose(axis_vector, [0, 0, -1]):
        v1 = np.array([1, 0, 0])
    else:
        v1 = np.array([0, 0, 1])
    # Compute perpendicular vectors
    v_east = np.cross(axis_vector, v1)
    v_east /= np.linalg.norm(v_east)
    v_north = np.cross(axis_vector, v_east)
    v_north /= np.linalg.norm(v_north)
    
    # Generate circle points around the base
    theta = np.linspace(0, 2 * np.pi, num_points)
    circle_points = base_center_vector[:, np.newaxis] + base_radius * (
        v_east[:, np.newaxis] * np.cos(theta) + v_north[:, np.newaxis] * np.sin(theta)
    )
    circle_points /= np.linalg.norm(circle_points, axis=0)  # Normalize to unit sphere
    
    # Convert circle points to RA and Dec
    ra_base_vertices, dec_base_vertices = cartesian_to_radec(
        circle_points[0], circle_points[1], circle_points[2]
    )
    
    # Return base vertices as SkyCoord
    base_vertices = SkyCoord(ra=ra_base_vertices * u.deg, dec=dec_base_vertices * u.deg, frame='icrs')
    return base_vertices

# Generate the base vertices
cone_base = cone_base_vertices(cone_top, cone_base_center, cone_opening_angle, cone_inclination)