void palPv2el	(	const double	pv[6],
		double	date,
		double	pmass,
		int	jformr,
		int *	jform,
		double *	epoch,
		double *	orbinc,
		double *	anode,
		double *	perih,
		double *	aorq,
		double *	e,
		double *	aorl,
		double *	dm,
		int *	jstat
	)

Definition at line 201 of file palPv2el.c.

References C, DMAX, eraAnp(), PAL__GCON, and PAL__SPD.

Referenced by palUe2el(), and t_planet().

                                                       {

  /*  Sin and cos of J2000 mean obliquity (IAU 1976) */
  const double SE = 0.3977771559319137;
  const double CE = 0.9174820620691818;

  /*  Minimum allowed distance (AU) and speed (AU/day) */
  const double RMIN = 1e-3;
  const double VMIN = 1e-8;

  /*  How close to unity the eccentricity has to be to call it a parabola */
  const double PARAB = 1.0e-8;

  double X,Y,Z,XD,YD,ZD,R,V2,V,RDV,GMU,HX,HY,HZ,
    HX2PY2,H2,H,OI,BIGOM,AR,ECC,S,C,AT,U,OM,
    GAR3,EM1,EP1,HAT,SHAT,CHAT,AE,AM,DN,PL,
    EL,Q,TP,THAT,THHF,F;

  int JF;

  /*  Validate arguments PMASS and JFORMR.*/
  if (pmass < 0.0) {
    *jstat = -1;
    return;
  }
  if (jformr < 1 || jformr > 3) {
    *jstat = -2;
    return;
  }

  /*  Provisionally assume the elements will be in the chosen form. */
  JF = jformr;

  /*  Rotate the position from equatorial to ecliptic coordinates. */
  X = pv[0];
  Y = pv[1]*CE+pv[2]*SE;
  Z = -pv[1]*SE+pv[2]*CE;

  /*  Rotate the velocity similarly, scaling to AU/day. */
  XD = PAL__SPD*pv[3];
  YD = PAL__SPD*(pv[4]*CE+pv[5]*SE);
  ZD = PAL__SPD*(-pv[4]*SE+pv[5]*CE);

  /*  Distance and speed. */
  R = sqrt(X*X+Y*Y+Z*Z);
  V2 = XD*XD+YD*YD+ZD*ZD;
  V = sqrt(V2);

  /*  Reject unreasonably small values. */
  if (R < RMIN || V < VMIN) {
    *jstat = -3;
    return;
  }

  /*  R dot V. */
  RDV = X*XD+Y*YD+Z*ZD;

  /*  Mu. */
  GMU = (1.0+pmass)*PAL__GCON*PAL__GCON;

  /*  Vector angular momentum per unit reduced mass. */
  HX = Y*ZD-Z*YD;
  HY = Z*XD-X*ZD;
  HZ = X*YD-Y*XD;

  /*  Areal constant. */
  HX2PY2 = HX*HX+HY*HY;
  H2 = HX2PY2+HZ*HZ;
  H = sqrt(H2);

  /*  Inclination. */
  OI = atan2(sqrt(HX2PY2),HZ);

  /*  Longitude of ascending node. */
  if (HX != 0.0 || HY != 0.0) {
    BIGOM = atan2(HX,-HY);
  } else {
    BIGOM=0.0;
  }

  /*  Reciprocal of mean distance etc. */
  AR = 2.0/R-V2/GMU;

  /*  Eccentricity. */
  ECC = sqrt(DMAX(1.0-AR*H2/GMU,0.0));

  /*  True anomaly. */
  S = H*RDV;
  C = H2-R*GMU;
  if (S != 0.0 || C != 0.0) {
    AT = atan2(S,C);
  } else {
    AT = 0.0;
  }

  /*  Argument of the latitude. */
  S = sin(BIGOM);
  C = cos(BIGOM);
  U = atan2((-X*S+Y*C)*cos(OI)+Z*sin(OI),X*C+Y*S);

  /*  Argument of perihelion. */
  OM = U-AT;

  /*  Capture near-parabolic cases. */
  if (fabs(ECC-1.0) < PARAB) ECC=1.0;

  /*  Comply with JFORMR = 1 or 2 only if orbit is elliptical. */
  if (ECC > 1.0) JF=3;

  /*  Functions. */
  GAR3 = GMU*AR*AR*AR;
  EM1 = ECC-1.0;
  EP1 = ECC+1.0;
  HAT = AT/2.0;
  SHAT = sin(HAT);
  CHAT = cos(HAT);

  /*  Variable initializations to avoid compiler warnings. */
  AM = 0.0;
  DN = 0.0;
  PL = 0.0;
  EL = 0.0;
  Q = 0.0;
  TP = 0.0;

  /*  Ellipse? */
  if (ECC < 1.0 ) {

    /*     Eccentric anomaly. */
    AE = 2.0*atan2(sqrt(-EM1)*SHAT,sqrt(EP1)*CHAT);

    /*     Mean anomaly. */
    AM = AE-ECC*sin(AE);

    /*     Daily motion. */
    DN = sqrt(GAR3);
  }

  /*  "Major planet" element set? */
  if (JF == 1) {

    /*     Longitude of perihelion. */
    PL = BIGOM+OM;

    /*     Longitude at epoch. */
    EL = PL+AM;
  }

  /*  "Comet" element set? */
  if (JF == 3) {

    /*     Perihelion distance. */
    Q = H2/(GMU*EP1);

    /*     Ellipse, parabola, hyperbola? */
    if (ECC < 1.0) {

      /*        Ellipse: epoch of perihelion. */
      TP = date-AM/DN;

    } else {

      /*        Parabola or hyperbola: evaluate tan ( ( true anomaly ) / 2 ) */
      THAT = SHAT/CHAT;
      if (ECC == 1.0) {

        /*           Parabola: epoch of perihelion. */
        TP = date-THAT*(1.0+THAT*THAT/3.0)*H*H2/(2.0*GMU*GMU);

      } else {

        /*           Hyperbola: epoch of perihelion. */
        THHF = sqrt(EM1/EP1)*THAT;
        F = log(1.0+THHF)-log(1.0-THHF);
        TP = date-(ECC*sinh(F)-F)/sqrt(-GAR3);
      }
    }
  }

  /*  Return the appropriate set of elements. */
  *jform = JF;
  *orbinc = OI;
  *anode = eraAnp(BIGOM);
  *e = ECC;
  if (JF == 1) {
    *perih = eraAnp(PL);
    *aorl = eraAnp(EL);
    *dm = DN;
  } else {
    *perih = eraAnp(OM);
    if (JF == 2) *aorl = eraAnp(AM);
  }
  if (JF != 3) {
    *epoch = date;
    *aorq = 1.0/AR;
  } else {
    *epoch = TP;
    *aorq = Q;
  }
  *jstat = 0;

}

Here is the call graph for this function:

Here is the caller graph for this function: