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

  char *type;
  int i, status, Polynomial, Zenith1, Zenith2, Zenithal, Cartesian, PseudoCyl;
  double phi, theta;
  double determ;
  double X, Y, L, M, Lo, Mo, dL, dM;
  double sphi, cphi, stht;
  double salp, calp, sdel, cdel, sdp, cdp;
  double P, A, Rc;

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

  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)  {
    sdp  = sin(Do);
    cdp  = cos(Do);
    salp = sin(ra - Ro);
    calp = cos(ra - Ro);
    sdel = sin(dec);
    cdel = cos(dec);

    stht = sdel*sdp + cdel*cdp*calp;    /* sin(theta) */
    sphi = cdel*salp;                   /* = cos(theta)*sin(phi) */
    cphi = cdel*sdp*calp - sdel*cdp;    /* = cos(theta)*cos(phi) */
    if (stht < 0) status = FALSE;

    if (Zenith1) {
      L =  sphi / stht;
      M = -cphi / stht;
    }
    if (Zenith2) {
      L =  sphi;
      M = -cphi;
    }
  }

  /**** Other Standard Projections ****/
  if (PseudoCyl) {
    if (proj->type === PS_PROJ_AIT) {
      phi = (ra - coords[0].crval1);
      theta = (dec - coords[0].crval2);
      P = 1.0 + cos (theta) * cos (0.5*phi);
      if (P != 0.0) {
	A =  sqrt (2.0 / P);
	L =  2.0 * A * cos (theta) * sin (0.5*phi);
	M =  A * sin (theta);
      } else { 
	L =  0.0;
	M =  0.0;
      }	
    }
    if (proj->type === PS_PROJ_GLS) {
      phi = ra - Ro;
      theta = dec - Do;
      L = phi * cos(theta);
      M = theta;
    }
    if (proj->type === PS_PROJ_PAR) {
      phi = ra - Ro;
      theta = dec - Do;
      L = phi * (2.0*cos(2*theta/3.0) - 1);
      M = 180.0 * sin (theta/3.0);
    }
  }

  *x = X / proj->Xs;
  *y = Y / proj->Ys;

  return (status);
}