void palUe2pv	(	double	date,
		double	u[13],
		double	pv[6],
		int *	jstat
	)

Definition at line 124 of file palUe2pv.c.

References PAL__GCON, and PAL__SPD.

                                                                     {

  /*  Canonical days to seconds */
  const double CD2S = PAL__GCON / PAL__SPD;

  /*  Test value for solution and maximum number of iterations */
  const double TEST = 1e-13;
  const int NITMAX = 25;

  int I, NIT, N;
  double CM,ALPHA,T0,P0[3],V0[3],R0,SIGMA0,T,PSI,DT,W,
    TOL,PSJ,PSJ2,BETA,S0,S1,S2,S3,
    FF,R,F,G,FD,GD;

  double PLAST = 0.0;
  double FLAST = 0.0;

  /*  Unpack the parameters. */
  CM = u[0];
  ALPHA = u[1];
  T0 = u[2];
  for (I=0; I<3; I++) {
    P0[I] = u[I+3];
    V0[I] = u[I+6];
  }
  R0 = u[9];
  SIGMA0 = u[10];
  T = u[11];
  PSI = u[12];

  /*  Approximately update the universal eccentric anomaly. */
  PSI = PSI+(date-T)*PAL__GCON/R0;

  /*  Time from reference epoch to date (in Canonical Days: a canonical
   *  day is 58.1324409... days, defined as 1/PAL__GCON). */
  DT = (date-T0)*PAL__GCON;

  /*  Refine the universal eccentric anomaly, psi. */
  NIT = 1;
  W = 1.0;
  TOL = 0.0;
  while (fabs(W) >= TOL) {

    /*     Form half angles until BETA small enough. */
    N = 0;
    PSJ = PSI;
    PSJ2 = PSJ*PSJ;
    BETA = ALPHA*PSJ2;
    while (fabs(BETA) > 0.7) {
      N = N+1;
      BETA = BETA/4.0;
      PSJ = PSJ/2.0;
      PSJ2 = PSJ2/4.0;
    }

    /*     Calculate Universal Variables S0,S1,S2,S3 by nested series. */
    S3 = PSJ*PSJ2*((((((BETA/210.0+1.0)
                       *BETA/156.0+1.0)
                      *BETA/110.0+1.0)
                     *BETA/72.0+1.0)
                    *BETA/42.0+1.0)
                   *BETA/20.0+1.0)/6.0;
    S2 = PSJ2*((((((BETA/182.0+1.0)
                   *BETA/132.0+1.0)
                  *BETA/90.0+1.0)
                 *BETA/56.0+1.0)
                *BETA/30.0+1.0)
               *BETA/12.0+1.0)/2.0;
    S1 = PSJ+ALPHA*S3;
    S0 = 1.0+ALPHA*S2;

    /*     Undo the angle-halving. */
    TOL = TEST;
    while (N > 0) {
      S3 = 2.0*(S0*S3+PSJ*S2);
      S2 = 2.0*S1*S1;
      S1 = 2.0*S0*S1;
      S0 = 2.0*S0*S0-1.0;
      PSJ = PSJ+PSJ;
      TOL += TOL;
      N--;
    }

    /*     Values of F and F' corresponding to the current value of psi. */
    FF = R0*S1+SIGMA0*S2+CM*S3-DT;
    R = R0*S0+SIGMA0*S1+CM*S2;

    /*     If first iteration, create dummy "last F". */
    if ( NIT == 1) FLAST = FF;

    /*     Check for sign change. */
    if ( FF*FLAST < 0.0 ) {

      /*        Sign change:  get psi adjustment using secant method. */
      W = FF*(PLAST-PSI)/(FLAST-FF);
    } else {

      /*        No sign change:  use Newton-Raphson method instead. */
      if (R == 0.0) {
        /* Null radius vector */
        *jstat = -1;
        return;
      }
      W = FF/R;
    }

    /*     Save the last psi and F values. */
    PLAST = PSI;
    FLAST = FF;

    /*     Apply the Newton-Raphson or secant adjustment to psi. */
    PSI = PSI-W;

    /*     Next iteration, unless too many already. */
    if (NIT > NITMAX) {
      *jstat = -2; /* Failed to converge */
      return;
    }
    NIT++;
  }

  /*  Project the position and velocity vectors (scaling velocity to AU/s). */
  W = CM*S2;
  F = 1.0-W/R0;
  G = DT-CM*S3;
  FD = -CM*S1/(R0*R);
  GD = 1.0-W/R;
  for (I=0; I<3; I++) {
    pv[I] = P0[I]*F+V0[I]*G;
    pv[I+3] = CD2S*(P0[I]*FD+V0[I]*GD);
  }

  /*  Update the parameters to allow speedy prediction of PSI next time. */
  u[11] = date;
  u[12] = PSI;

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