def dispos(dra0, decd0, dra, decd):
    """Compute distance and position angle solving a spherical
    triangle (no approximations).

    Source/credit: Skycat
    Author: A.P. Martinez

    Parameters
    ----------
    dra0 : float
        Center RA in decimal degrees.

    decd0 : float
        Center DEC in decimal degrees.

    dra : float
        Point RA in decimal degrees.

    decd : float
        Point DEC in decimal degrees.

    Returns
    -------
    phi : float
        Phi in degrees (East of North).

    dist : float
        Distance in arcmin.

    """
    radian = 180.0 / math.pi

    # coord transformed in radians
    alf = dra / radian
    alf0 = dra0 / radian
    del_ = decd / radian
    del0 = decd0 / radian

    sd0 = math.sin(del0)
    sd = math.sin(del_)
    cd0 = math.cos(del0)
    cd = math.cos(del_)
    cosda = math.cos(alf - alf0)
    cosd = sd0 * sd + cd0 * cd * cosda
    dist = math.acos(cosd)
    phi = 0.0
    if dist > 0.0000004:
        sind = math.sin(dist)
        cospa = (sd * cd0 - cd * sd0 * cosda) / sind
        #if cospa > 1.0:
        #    cospa=1.0
        if math.fabs(cospa) > 1.0:
            # 2005-06-02: fix from awicenec@eso.org
            cospa = cospa / math.fabs(cospa)
        sinpa = cd * math.sin(alf - alf0) / sind
        phi = math.acos(cospa) * radian
        if sinpa < 0.0:
            phi = 360.0 - phi
    dist *= radian
    dist *= 60.0
    if decd0 == 90.0:
        phi = 180.0
    if decd0 == -90.0:
        phi = 0.0
    return (phi, dist)
Compute distance and position angle solving a spherical
    triangle (no approximations).

    Source/credit: Skycat
    Author: A.P. Martinez

    Parameters
    ----------
    dra0 : float
        Center RA in decimal degrees.

    decd0 : float
        Center DEC in decimal degrees.

    dra : float
        Point RA in decimal degrees.

    decd : float
        Point DEC in decimal degrees.

    Returns
    -------
    phi : float
        Phi in degrees (East of North).

    dist : float
        Distance in arcmin.