1 /* kepleq.f -- translated by f2c (version 19980913).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* $Procedure KEPLEQ ( Kepler's Equation - Equinoctial Version ) */
kepleq_(doublereal * ml,doublereal * h__,doublereal * k)9 doublereal kepleq_(doublereal *ml, doublereal *h__, doublereal *k)
10 {
11 /* System generated locals */
12 doublereal ret_val;
13
14 /* Builtin functions */
15 double cos(doublereal), sin(doublereal);
16
17 /* Local variables */
18 doublereal evec[2];
19 extern /* Subroutine */ int chkin_(char *, ftnlen), errdp_(char *,
20 doublereal *, ftnlen);
21 doublereal e2;
22 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
23 ftnlen), setmsg_(char *, ftnlen);
24 extern doublereal kpsolv_(doublereal *);
25
26 /* $ Abstract */
27
28 /* This function solves the equinoctial version of Kepler's */
29 /* equation. */
30
31 /* $ Disclaimer */
32
33 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
34 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
35 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
36 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
37 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
38 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
39 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
40 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
41 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
42 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
43
44 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
45 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
46 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
47 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
48 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
49 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
50
51 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
52 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
53 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
54 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
55
56 /* $ Required_Reading */
57
58 /* None. */
59
60 /* $ Keywords */
61
62 /* SPK */
63
64 /* $ Declarations */
65 /* $ Brief_I/O */
66
67 /* VARIABLE I/O DESCRIPTION */
68 /* -------- --- -------------------------------------------------- */
69 /* ML I Mean longitude */
70 /* H I h component of equinoctial elements */
71 /* K I k component of equinoctial elements */
72
73 /* $ Detailed_Input */
74
75 /* ML mean longitude of some body following two body */
76 /* motion. (Mean longitude = Mean anomaly + argument */
77 /* of periapse + longitude of ascending node.) */
78
79 /* H The h component of the equinoctial element set */
80 /* ( h = ECC*SIN( arg of periapse + long ascending node) ) */
81
82 /* K The k component of the equinoctial element set */
83 /* ( k = ECC*COS( arg of periapse + long ascending node) ) */
84
85 /* Note that ECC = DSQRT ( K*K + H*H ) */
86
87 /* $ Detailed_Output */
88
89 /* The function returns the value of F such that */
90 /* ML = F + h*COS(F) - k*SIN(F) */
91
92 /* $ Parameters */
93
94 /* None. */
95
96 /* $ Files */
97
98 /* None. */
99
100 /* $ Exceptions */
101
102 /* 1) If the sum of the squares of F and K is not less than .9 */
103 /* the error 'SPICE(ECCOUTOFBOUNDS)' will be signalled. */
104
105 /* 2) If the iteration for a solution to the equinoctial Kepler's */
106 /* equation does not converge in 10 or fewer steps, the error */
107 /* 'SPICE(NOCONVERGENCE)' is signalled. */
108
109 /* $ Particulars */
110
111 /* This routine solves the equinoctial element version of */
112 /* Kepler's equation. */
113
114 /* ML = F + h*COS(F) - k*SIN(F) */
115
116 /* Here F is an offset from the eccentric anomaly E. */
117
118 /* F = E - argument of periapse - longitude of ascending node. */
119
120 /* where E is eccentric anomaly. */
121
122 /* $ Examples */
123
124 /* None. */
125
126 /* $ Restrictions */
127
128 /* None. */
129
130 /* $ Author_and_Institution */
131
132 /* W.L. Taber (JPL) */
133
134 /* $ Literature_References */
135
136 /* "Optical Navigation Program Mathematical Models" JPL */
137 /* Engineering Memorandum 314-513. By William M. Owen */
138 /* August 9, 1991. */
139
140 /* $ Version */
141
142 /* - SPICELIB Version 1.0.0, 11-DEC-1996 (WLT) */
143
144
145 /* -& */
146 /* $ Index_Entries */
147
148 /* Solve the equinoctial version of Kepler's equation */
149
150 /* -& */
151
152 /* SPICELIB Functions */
153
154
155 /* Local variables */
156
157
158 /* Make sure that H and K are in the expected range. */
159
160 e2 = *h__ * *h__ + *k * *k;
161 if (e2 >= .81) {
162 ret_val = 0.;
163 chkin_("KEPLEQ", (ftnlen)6);
164 setmsg_("The values of H and K supplied to KEPLEQ must satisfy the i"
165 "nequality H*H + K*K < ECC**2 where ECC is the eccentricity t"
166 "hreshold of 0.9. The values of H and K are: # and # respect"
167 "ively. H*H + K*K = #. ", (ftnlen)201);
168 errdp_("#", h__, (ftnlen)1);
169 errdp_("#", k, (ftnlen)1);
170 errdp_("#", &e2, (ftnlen)1);
171 sigerr_("SPICE(ECCOUTOFBOUNDS)", (ftnlen)21);
172 chkout_("KEPLEQ", (ftnlen)6);
173 return ret_val;
174 }
175
176 /* Instead of solving the equation */
177
178 /* ML = F + H*DCOS(F) - K*DSIN(F) */
179
180 /* We set X equal to F - ML and solve the equivalent equation */
181
182 /* 0 = X + H*DCOS(ML+X) - K*DSIN(ML+X) */
183
184 /* = X + H*{DCOS(ML)*DCOS(X) - DSIN(ML)*DSIN(X)} */
185 /* - K*{DSIN(ML)*DCOS(X) + DCOS(ML)*DSIN(X)} */
186
187 /* = X + { H*DCOS(ML) - K*DSIN(ML) }*DCOS(X) */
188 /* - { H*DSIN(ML) + K*DCOS(ML) }*DSIN(X) */
189
190
191 /* We can rearrange this to: */
192
193 /* - - - - */
194 /* | DCOS(ML) -DSIN(ML) | | DCOS(X) | */
195 /* 0 = X + [ H -K ] * | DSIN(ML) DCOS(ML) | * | DSIN(X) | */
196 /* - - - - */
197
198 /* Finally if we let */
199
200 /* C - - */
201 /* | DCOS(ML) -DSIN(ML) | */
202 /* EVEC = [ EX EY ] = [ -H K ] * | DSIN(ML) DCOS(ML) | */
203 /* - - */
204
205 /* and */
206
207 /* DCOS(X) */
208 /* U(X) = DSIN(X) */
209
210 /* Then we can rewrite the equation as: */
211
212 /* 0 = X - < EVEC, U(X) > */
213
214 /* where <,> denotes the dot product operation. Note that X */
215 /* is necessarily in the range from -ECC to ECC where ECC = | EVEC | */
216
217 /* Once we've computed X, F is just ML + X. */
218
219 /* For those of you who are fans of the classical keplerian */
220 /* elements: */
221
222 /* x = F - ML = E - M */
223
224 /* where E denotes eccentric anomaly and M denotes mean anomaly. */
225
226 /* The routine KPEVEC returns the value of X that solves */
227 /* the equation X - < EVEC, UVEC(X) > */
228
229 evec[0] = -(*h__) * cos(*ml) + *k * sin(*ml);
230 evec[1] = *h__ * sin(*ml) + *k * cos(*ml);
231 ret_val = *ml + kpsolv_(evec);
232 return ret_val;
233 } /* kepleq_ */
234
235