void palEl2ue	(	double	date,
		int	jform,
		double	epoch,
		double	orbinc,
		double	anode,
		double	perih,
		double	aorq,
		double	e,
		double	aorl,
		double	dm,
		double	u[13],
		int *	jstat
	)

Definition at line 178 of file palEl2ue.c.

References PAL__GCON, palPv2ue(), and palUe2pv().

Referenced by palPertel(), palPlanel(), palPlante(), and t_planet().

                                                                   {

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

  int J;

  double PHT,ARGPH,Q,W,CM,ALPHA,PHS,SW,CW,SI,CI,SO,CO,
    X,Y,Z,PX,PY,PZ,VX,VY,VZ,DT,FC,FP,PSI,
    UL[13],PV[6];

  /*  Validate arguments. */
  if (jform < 1 || jform > 3) {
    *jstat = -1;
    return;
  }
  if (e < 0.0 || e > 10.0 || (e >= 1.0 && jform != 3)) {
    *jstat = -2;
    return;
  }
  if (aorq <= 0.0) {
    *jstat = -3;
    return;
  }
  if (jform == 1 && dm <= 0.0) {
    *jstat = -4;
    return;
  }

  /*
   *  Transform elements into standard form:
   *
   *  PHT   = epoch of perihelion passage
   *  ARGPH = argument of perihelion (little omega)
   *  Q     = perihelion distance (q)
   *  CM    = combined mass, M+m (mu)
   */

  if (jform == 1) {

    /*     Major planet. */
    PHT = epoch-(aorl-perih)/dm;
    ARGPH = perih-anode;
    Q = aorq*(1.0-e);
    W = dm/PAL__GCON;
    CM = W*W*aorq*aorq*aorq;

  } else if (jform == 2) {

    /*     Minor planet. */
    PHT = epoch-aorl*sqrt(aorq*aorq*aorq)/PAL__GCON;
    ARGPH = perih;
    Q = aorq*(1.0-e);
    CM = 1.0;

  } else {

    /*     Comet. */
    PHT = epoch;
    ARGPH = perih;
    Q = aorq;
    CM = 1.0;

  }

  /*  The universal variable alpha.  This is proportional to the total
   *  energy of the orbit:  -ve for an ellipse, zero for a parabola,
   *  +ve for a hyperbola. */

  ALPHA = CM*(e-1.0)/Q;

  /*  Speed at perihelion. */

  PHS = sqrt(ALPHA+2.0*CM/Q);

  /*  In a Cartesian coordinate system which has the x-axis pointing
   *  to perihelion and the z-axis normal to the orbit (such that the
   *  object orbits counter-clockwise as seen from +ve z), the
   *  perihelion position and velocity vectors are:
   *
   *    position   [Q,0,0]
   *    velocity   [0,PHS,0]
   *
   *  To express the results in J2000 equatorial coordinates we make a
   *  series of four rotations of the Cartesian axes:
   *
   *           axis      Euler angle
   *
   *     1      z        argument of perihelion (little omega)
   *     2      x        inclination (i)
   *     3      z        longitude of the ascending node (big omega)
   *     4      x        J2000 obliquity (epsilon)
   *
   *  In each case the rotation is clockwise as seen from the +ve end of
   *  the axis concerned.
   */

  /*  Functions of the Euler angles. */
  SW = sin(ARGPH);
  CW = cos(ARGPH);
  SI = sin(orbinc);
  CI = cos(orbinc);
  SO = sin(anode);
  CO = cos(anode);

  /*  Position at perihelion (AU). */
  X = Q*CW;
  Y = Q*SW;
  Z = Y*SI;
  Y = Y*CI;
  PX = X*CO-Y*SO;
  Y = X*SO+Y*CO;
  PY = Y*CE-Z*SE;
  PZ = Y*SE+Z*CE;

  /*  Velocity at perihelion (AU per canonical day). */
  X = -PHS*SW;
  Y = PHS*CW;
  Z = Y*SI;
  Y = Y*CI;
  VX = X*CO-Y*SO;
  Y = X*SO+Y*CO;
  VY = Y*CE-Z*SE;
  VZ = Y*SE+Z*CE;

  /*  Time from perihelion to date (in Canonical Days: a canonical day
   *  is 58.1324409... days, defined as 1/PAL__GCON). */

  DT = (date-PHT)*PAL__GCON;

  /*  First approximation to the Universal Eccentric Anomaly, PSI,
   *  based on the circle (FC) and parabola (FP) values. */

  FC = DT/Q;
  W = pow(3.0*DT+sqrt(9.0*DT*DT+8.0*Q*Q*Q), 1.0/3.0);
  FP = W-2.0*Q/W;
  PSI = (1.0-e)*FC+e*FP;

  /*  Assemble local copy of element set. */
  UL[0] = CM;
  UL[1] = ALPHA;
  UL[2] = PHT;
  UL[3] = PX;
  UL[4] = PY;
  UL[5] = PZ;
  UL[6] = VX;
  UL[7] = VY;
  UL[8] = VZ;
  UL[9] = Q;
  UL[10] = 0.0;
  UL[11] = date;
  UL[12] = PSI;

  /*  Predict position+velocity at epoch of osculation. */
  palUe2pv( date, UL, PV, &J );
  if (J != 0) {
    *jstat = -5;
    return;
  }

  /*  Convert back to universal elements. */
  palPv2ue( PV, date, CM-1.0, u, &J );
  if (J != 0) {
    *jstat = -5;
    return;
  }

  /*  OK exit. */
  *jstat = 0;

}

Here is the call graph for this function:

Here is the caller graph for this function: