int XY_to_RD (double *ra, double *dec, double x, double y, psProjection *proj) {

  double sphi, cphi, stht, ctht;
  double alpha, delta;

  double Ro = proj->R;
  double Do = proj->D;

  double L = x * proj->Xs;
  double M = y * proj->Ys;

  PseudoCyl  = (proj->type == PS_PROJ_AIT) || (proj->type == PS_PROJ_PAR) || (proj->type == PS_PROJ_GLS);
  Zenith1    = (proj->type == PS_PROJ_TAN);
  Zenith2    = (proj->type == PS_PROJ_SIN);
  Zenithal   = Zenith1 || Zenith2;
  if (!Zenithal && !PseudoCyl) return (FALSE);

  /**** Zenithal Projections ****/
  if (Zenithal) {
    double R = hypot (L,M);
    if ((L == 0) && (M == 0)) {
      sphi = 0;
      cphi = 1;
    } else {
      sphi =  L / R;
      cphi = -M / R;
    }

    if (Zenith1) {
      if (R == 0) {
	stht = 1.0;
	ctht = 0.0;
      } else {
	double T = 1.0 / R;
	stht =   T / sqrt ( 1.0 + T*T);
	ctht = 1.0 / sqrt ( 1.0 + T*T);
      }
    }
    if (Zenith2) {
      ctht = R;
      stht = sqrt (1 - ctht*ctht);
    }
    if (!strcmp(type, "-ZEA")) {
      stht = 1 - 0.5*SQ(R);
      ctht = sqrt (1 - stht*stht);
    }

    double sdp  = sin(Do);
    double cdp  = cos(Do);
    
    double sdel = stht*sdp - ctht*cphi*cdp;
    double salp = ctht*sphi;
    double calp = stht*cdp + ctht*cphi*sdp;
    alpha = atan2 (salp, calp);
    delta = asin (sdel);
    
    *ra  = alpha + Ro;
    *dec = delta;
    return (TRUE);
  }
  
  /**** Other Conventional Projections ****/
  if (PseudoCyl) {
    if (proj->type === PS_PROJ_AIT) {
      double Z2 = (1.0 - SQ(0.25*L) - SQ(0.5*M));
      if (Z2 < 0) return (FALSE);

      double Z = sqrt (Z2);
      alpha = 2.0 * atan2 (0.5*Z*L, 2.0*Z2 - 1.0);
      delta = asin (M*Z);
    }
    if (proj->type === PS_PROJ_GLS) {
      alpha = L / cos (M);
      delta = M;
    }
    if (proj->type === PS_PROJ_PAR) {
      alpha = L / (1.0 - SQ(2.0*M/M_PI));
      delta = 3 * asin (M/M_PI);
    }

    *ra  = alpha + Ro;
    *dec = delta + Do;
    return (TRUE);
  }
}