| 1 | /* | 
|---|
| 2 | *+ | 
|---|
| 3 | *  Name: | 
|---|
| 4 | *     palPlanel | 
|---|
| 5 |  | 
|---|
| 6 | *  Purpose: | 
|---|
| 7 | *     Transform conventional elements into position and velocity | 
|---|
| 8 |  | 
|---|
| 9 | *  Language: | 
|---|
| 10 | *     Starlink ANSI C | 
|---|
| 11 |  | 
|---|
| 12 | *  Type of Module: | 
|---|
| 13 | *     Library routine | 
|---|
| 14 |  | 
|---|
| 15 | *  Invocation: | 
|---|
| 16 | *     void palPlanel ( double date, int jform, double epoch, double orbinc, | 
|---|
| 17 | *                      double anode, double perih, double aorq, double e, | 
|---|
| 18 | *                      double aorl, double dm, double pv[6], int *jstat ); | 
|---|
| 19 |  | 
|---|
| 20 | *  Arguments: | 
|---|
| 21 | *     date = double (Given) | 
|---|
| 22 | *        Epoch (TT MJD) of osculation (Note 1) | 
|---|
| 23 | *     jform = int (Given) | 
|---|
| 24 | *        Element set actually returned (1-3; Note 3) | 
|---|
| 25 | *     epoch = double (Given) | 
|---|
| 26 | *        Epoch of elements (TT MJD) (Note 4) | 
|---|
| 27 | *     orbinc = double (Given) | 
|---|
| 28 | *        inclination (radians) | 
|---|
| 29 | *     anode = double (Given) | 
|---|
| 30 | *        longitude of the ascending node (radians) | 
|---|
| 31 | *     perih = double (Given) | 
|---|
| 32 | *        longitude or argument of perihelion (radians) | 
|---|
| 33 | *     aorq = double (Given) | 
|---|
| 34 | *        mean distance or perihelion distance (AU) | 
|---|
| 35 | *     e = double (Given) | 
|---|
| 36 | *        eccentricity | 
|---|
| 37 | *     aorl = double (Given) | 
|---|
| 38 | *        mean anomaly or longitude (radians, JFORM=1,2 only) | 
|---|
| 39 | *     dm = double (Given) | 
|---|
| 40 | *        daily motion (radians, JFORM=1 only) | 
|---|
| 41 | *     u = double [13] (Returned) | 
|---|
| 42 | *        Universal orbital elements (Note 1) | 
|---|
| 43 | *            (0)  combined mass (M+m) | 
|---|
| 44 | *            (1)  total energy of the orbit (alpha) | 
|---|
| 45 | *            (2)  reference (osculating) epoch (t0) | 
|---|
| 46 | *          (3-5)  position at reference epoch (r0) | 
|---|
| 47 | *          (6-8)  velocity at reference epoch (v0) | 
|---|
| 48 | *            (9)  heliocentric distance at reference epoch | 
|---|
| 49 | *           (10)  r0.v0 | 
|---|
| 50 | *           (11)  date (t) | 
|---|
| 51 | *           (12)  universal eccentric anomaly (psi) of date, approx | 
|---|
| 52 | *     jstat = int * (Returned) | 
|---|
| 53 | *        status:  0 = OK | 
|---|
| 54 | *              - -1 = illegal JFORM | 
|---|
| 55 | *              - -2 = illegal E | 
|---|
| 56 | *              - -3 = illegal AORQ | 
|---|
| 57 | *              - -4 = illegal DM | 
|---|
| 58 | *              - -5 = numerical error | 
|---|
| 59 |  | 
|---|
| 60 | *  Description: | 
|---|
| 61 | *     Heliocentric position and velocity of a planet, asteroid or comet, | 
|---|
| 62 | *     starting from orbital elements. | 
|---|
| 63 |  | 
|---|
| 64 | *  Authors: | 
|---|
| 65 | *     PTW: Pat Wallace (STFC) | 
|---|
| 66 | *     TIMJ: Tim Jenness (JAC, Hawaii) | 
|---|
| 67 | *     {enter_new_authors_here} | 
|---|
| 68 |  | 
|---|
| 69 | *  Notes: | 
|---|
| 70 | *     - DATE is the instant for which the prediction is required.  It is | 
|---|
| 71 | *       in the TT timescale (formerly Ephemeris Time, ET) and is a | 
|---|
| 72 | *       Modified Julian Date (JD-2400000.5). | 
|---|
| 73 | *     - The elements are with respect to the J2000 ecliptic and equinox. | 
|---|
| 74 | *     - A choice of three different element-set options is available: | 
|---|
| 75 | * | 
|---|
| 76 | *       Option JFORM = 1, suitable for the major planets: | 
|---|
| 77 | * | 
|---|
| 78 | *         EPOCH  = epoch of elements (TT MJD) | 
|---|
| 79 | *         ORBINC = inclination i (radians) | 
|---|
| 80 | *         ANODE  = longitude of the ascending node, big omega (radians) | 
|---|
| 81 | *         PERIH  = longitude of perihelion, curly pi (radians) | 
|---|
| 82 | *         AORQ   = mean distance, a (AU) | 
|---|
| 83 | *         E      = eccentricity, e (range 0 to <1) | 
|---|
| 84 | *         AORL   = mean longitude L (radians) | 
|---|
| 85 | *         DM     = daily motion (radians) | 
|---|
| 86 | * | 
|---|
| 87 | *       Option JFORM = 2, suitable for minor planets: | 
|---|
| 88 | * | 
|---|
| 89 | *         EPOCH  = epoch of elements (TT MJD) | 
|---|
| 90 | *         ORBINC = inclination i (radians) | 
|---|
| 91 | *         ANODE  = longitude of the ascending node, big omega (radians) | 
|---|
| 92 | *         PERIH  = argument of perihelion, little omega (radians) | 
|---|
| 93 | *         AORQ   = mean distance, a (AU) | 
|---|
| 94 | *         E      = eccentricity, e (range 0 to <1) | 
|---|
| 95 | *         AORL   = mean anomaly M (radians) | 
|---|
| 96 | * | 
|---|
| 97 | *       Option JFORM = 3, suitable for comets: | 
|---|
| 98 | * | 
|---|
| 99 | *         EPOCH  = epoch of elements and perihelion (TT MJD) | 
|---|
| 100 | *         ORBINC = inclination i (radians) | 
|---|
| 101 | *         ANODE  = longitude of the ascending node, big omega (radians) | 
|---|
| 102 | *         PERIH  = argument of perihelion, little omega (radians) | 
|---|
| 103 | *         AORQ   = perihelion distance, q (AU) | 
|---|
| 104 | *         E      = eccentricity, e (range 0 to 10) | 
|---|
| 105 | * | 
|---|
| 106 | *       Unused arguments (DM for JFORM=2, AORL and DM for JFORM=3) are not | 
|---|
| 107 | *       accessed. | 
|---|
| 108 | *     - Each of the three element sets defines an unperturbed heliocentric | 
|---|
| 109 | *       orbit.  For a given epoch of observation, the position of the body | 
|---|
| 110 | *       in its orbit can be predicted from these elements, which are | 
|---|
| 111 | *       called "osculating elements", using standard two-body analytical | 
|---|
| 112 | *       solutions.  However, due to planetary perturbations, a given set | 
|---|
| 113 | *       of osculating elements remains usable for only as long as the | 
|---|
| 114 | *       unperturbed orbit that it describes is an adequate approximation | 
|---|
| 115 | *       to reality.  Attached to such a set of elements is a date called | 
|---|
| 116 | *       the "osculating epoch", at which the elements are, momentarily, | 
|---|
| 117 | *       a perfect representation of the instantaneous position and | 
|---|
| 118 | *       velocity of the body. | 
|---|
| 119 | * | 
|---|
| 120 | *       Therefore, for any given problem there are up to three different | 
|---|
| 121 | *       epochs in play, and it is vital to distinguish clearly between | 
|---|
| 122 | *       them: | 
|---|
| 123 | * | 
|---|
| 124 | *       . The epoch of observation:  the moment in time for which the | 
|---|
| 125 | *         position of the body is to be predicted. | 
|---|
| 126 | * | 
|---|
| 127 | *       . The epoch defining the position of the body:  the moment in time | 
|---|
| 128 | *         at which, in the absence of purturbations, the specified | 
|---|
| 129 | *         position (mean longitude, mean anomaly, or perihelion) is | 
|---|
| 130 | *         reached. | 
|---|
| 131 | * | 
|---|
| 132 | *       . The osculating epoch:  the moment in time at which the given | 
|---|
| 133 | *         elements are correct. | 
|---|
| 134 | * | 
|---|
| 135 | *       For the major-planet and minor-planet cases it is usual to make | 
|---|
| 136 | *       the epoch that defines the position of the body the same as the | 
|---|
| 137 | *       epoch of osculation.  Thus, only two different epochs are | 
|---|
| 138 | *       involved:  the epoch of the elements and the epoch of observation. | 
|---|
| 139 | * | 
|---|
| 140 | *       For comets, the epoch of perihelion fixes the position in the | 
|---|
| 141 | *       orbit and in general a different epoch of osculation will be | 
|---|
| 142 | *       chosen.  Thus, all three types of epoch are involved. | 
|---|
| 143 | * | 
|---|
| 144 | *       For the present routine: | 
|---|
| 145 | * | 
|---|
| 146 | *       . The epoch of observation is the argument DATE. | 
|---|
| 147 | * | 
|---|
| 148 | *       . The epoch defining the position of the body is the argument | 
|---|
| 149 | *         EPOCH. | 
|---|
| 150 | * | 
|---|
| 151 | *       . The osculating epoch is not used and is assumed to be close | 
|---|
| 152 | *         enough to the epoch of observation to deliver adequate accuracy. | 
|---|
| 153 | *         If not, a preliminary call to palPertel may be used to update | 
|---|
| 154 | *         the element-set (and its associated osculating epoch) by | 
|---|
| 155 | *         applying planetary perturbations. | 
|---|
| 156 | *     - The reference frame for the result is with respect to the mean | 
|---|
| 157 | *       equator and equinox of epoch J2000. | 
|---|
| 158 | *     - The algorithm was originally adapted from the EPHSLA program of | 
|---|
| 159 | *       D.H.P.Jones (private communication, 1996).  The method is based | 
|---|
| 160 | *       on Stumpff's Universal Variables. | 
|---|
| 161 |  | 
|---|
| 162 | *  See Also: | 
|---|
| 163 | *     Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983. | 
|---|
| 164 |  | 
|---|
| 165 | *  History: | 
|---|
| 166 | *     2012-03-12 (TIMJ): | 
|---|
| 167 | *        Initial version taken directly from SLA/F. | 
|---|
| 168 | *        Adapted with permission from the Fortran SLALIB library. | 
|---|
| 169 | *     {enter_further_changes_here} | 
|---|
| 170 |  | 
|---|
| 171 | *  Copyright: | 
|---|
| 172 | *     Copyright (C) 2002 Rutherford Appleton Laboratory | 
|---|
| 173 | *     Copyright (C) 2012 Science and Technology Facilities Council. | 
|---|
| 174 | *     All Rights Reserved. | 
|---|
| 175 |  | 
|---|
| 176 | *  Licence: | 
|---|
| 177 | *     This program is free software; you can redistribute it and/or | 
|---|
| 178 | *     modify it under the terms of the GNU General Public License as | 
|---|
| 179 | *     published by the Free Software Foundation; either version 3 of | 
|---|
| 180 | *     the License, or (at your option) any later version. | 
|---|
| 181 | * | 
|---|
| 182 | *     This program is distributed in the hope that it will be | 
|---|
| 183 | *     useful, but WITHOUT ANY WARRANTY; without even the implied | 
|---|
| 184 | *     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR | 
|---|
| 185 | *     PURPOSE. See the GNU General Public License for more details. | 
|---|
| 186 | * | 
|---|
| 187 | *     You should have received a copy of the GNU General Public License | 
|---|
| 188 | *     along with this program; if not, write to the Free Software | 
|---|
| 189 | *     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, | 
|---|
| 190 | *     MA 02110-1301, USA. | 
|---|
| 191 |  | 
|---|
| 192 | *  Bugs: | 
|---|
| 193 | *     {note_any_bugs_here} | 
|---|
| 194 | *- | 
|---|
| 195 | */ | 
|---|
| 196 |  | 
|---|
| 197 | #include "pal.h" | 
|---|
| 198 |  | 
|---|
| 199 |  | 
|---|
| 200 | void palPlanel ( double date, int jform, double epoch, double orbinc, | 
|---|
| 201 | double anode, double perih, double aorq, double e, | 
|---|
| 202 | double aorl, double dm, double pv[6], int *jstat ) { | 
|---|
| 203 |  | 
|---|
| 204 | int j; | 
|---|
| 205 | double u[13]; | 
|---|
| 206 |  | 
|---|
| 207 | /*  Validate elements and convert to "universal variables" parameters. */ | 
|---|
| 208 | palEl2ue( date, jform, epoch, orbinc, anode, perih, aorq, e, aorl, | 
|---|
| 209 | dm, u, &j ); | 
|---|
| 210 |  | 
|---|
| 211 | /* Determine the position and velocity */ | 
|---|
| 212 | if (j == 0) { | 
|---|
| 213 | palUe2pv( date, u, pv, &j); | 
|---|
| 214 | if (j != 0) j = -5; | 
|---|
| 215 | } | 
|---|
| 216 |  | 
|---|
| 217 | /* Wrap up */ | 
|---|
| 218 | *jstat = j; | 
|---|
| 219 |  | 
|---|
| 220 | } | 
|---|