code: def dispos(dra0, decd0, dra, decd):
      
      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)
docs: 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.