view release on metacpan or  search on metacpan

palsrc/palUe2pv.c  view on Meta::CPAN

*       multiple times, which is faster than multiple calls to palPlanel.
*     - It is not obligatory to use palEl2ue to obtain the parameters.
*       However, it should be noted that because palEl2ue performs its
*       own validation, no checks on the contents of the array U are made
*       by the present routine.
*     - DATE is the instant for which the prediction is required.  It is
*       in the TT timescale (formerly Ephemeris Time, ET) and is a
*       Modified Julian Date (JD-2400000.5).
*     - The universal elements supplied in the array U are in canonical
*       units (solar masses, AU and canonical days).  The position and
*       velocity are not sensitive to the choice of reference frame.  The
*       palEl2ue routine in fact produces coordinates with respect to the
*       J2000 equator and equinox.
*     - The algorithm was originally adapted from the EPHSLA program of
*       D.H.P.Jones (private communication, 1996).  The method is based
*       on Stumpff's Universal Variables.
*     - Reference:  Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.

*  History:
*     2012-03-09 (TIMJ):
*        Initial version cloned from SLA/F.
*        Adapted with permission from the Fortran SLALIB library.
*     {enter_further_changes_here}

*  Copyright:
*     Copyright (C) 2005 Rutherford Appleton Laboratory
*     Copyright (C) 2012 Science and Technology Facilities Council.
*     All Rights Reserved.

*  Licence:
*     This program is free software; you can redistribute it and/or
*     modify it under the terms of the GNU General Public License as
*     published by the Free Software Foundation; either version 3 of
*     the License, or (at your option) any later version.
*
*     This program is distributed in the hope that it will be
*     useful, but WITHOUT ANY WARRANTY; without even the implied
*     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
*     PURPOSE. See the GNU General Public License for more details.
*
*     You should have received a copy of the GNU General Public License
*     along with this program; if not, write to the Free Software
*     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
*     MA 02110-1301, USA.

*  Bugs:
*     {note_any_bugs_here}
*-
*/

#include <math.h>

#include "pal.h"
#include "palmac.h"

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

  /*  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)