import math


def point_in_ellipse(origin, point, a, b, pa_rad, verbose=False):
    """
    Identify if the point is inside the ellipse.

    :param origin A SkyCoord defining the centre of the ellipse.
    :param point A SkyCoord defining the point to be checked.
    :param a The semi-major axis in arcsec of the ellipse
    :param b The semi-minor axis in arcsec of the ellipse
    :param pa_rad The position angle of the ellipse. This is the angle of the major axis measured in radians East of 
                   North (or CCW from the y axis).
    """
    # Convert point to be in plane of the ellipse, accounting for distortions at high declinations
    p_ra_dist = (point.icrs.ra.degree - origin.icrs.ra.degree)* math.cos(origin.icrs.dec.rad)
    p_dec_dist = point.icrs.dec.degree - origin.icrs.dec.degree

    # Calculate the angle and radius of the test opoint relative to the centre of the ellipse
    # Note that we reverse the ra direction to reflect the CCW direction
    radius = math.sqrt(p_ra_dist**2 + p_dec_dist**2)
    diff_angle = (math.pi/2 + pa_rad) if p_dec_dist == 0 else math.atan(p_ra_dist / p_dec_dist) - pa_rad

    # Obtain the point position in terms of the ellipse major and minor axes
    minor = radius * math.sin(diff_angle)
    major = radius * math.cos(diff_angle)
    if verbose:
        print ('point relative to ellipse centre angle:{} deg radius:{:.4f}" maj:{:.2f}" min:{:.2f}"'.format(math.degrees(diff_angle), radius*3600, 
                major*3600, minor*3600))
    
    a_deg = a / 3600.0
    b_deg = b / 3600.0

    # Calc distance from origin relative to a and b
    dist = math.sqrt((major / a_deg) ** 2 + (minor / b_deg) ** 2)
    if verbose:
        print("Point %s is %f from ellipse %f, %f, %f at %s." % (point, dist, a, b, math.degrees(pa_rad), origin))
    return round(dist,3) <= 1.0