1 /* edlimb.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 /* Table of constant values */
9
10 static doublereal c_b18 = 1.;
11 static integer c__9 = 9;
12
13 /* $Procedure EDLIMB ( Ellipsoid Limb ) */
edlimb_(doublereal * a,doublereal * b,doublereal * c__,doublereal * viewpt,doublereal * limb)14 /* Subroutine */ int edlimb_(doublereal *a, doublereal *b, doublereal *c__,
15 doublereal *viewpt, doublereal *limb)
16 {
17 /* System generated locals */
18 doublereal d__1, d__2, d__3;
19
20 /* Local variables */
21 doublereal scla, sclb, sclc;
22 extern /* Subroutine */ int vscl_(doublereal *, doublereal *, doublereal *
23 );
24 doublereal scla2, sclb2, sclc2, v[3], scale;
25 extern /* Subroutine */ int chkin_(char *, ftnlen);
26 doublereal level;
27 extern /* Subroutine */ int moved_(doublereal *, integer *, doublereal *);
28 logical found;
29 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen), vsclg_(
30 doublereal *, doublereal *, integer *, doublereal *);
31 doublereal tmpel[9];
32 extern /* Subroutine */ int nvc2pl_(doublereal *, doublereal *,
33 doublereal *);
34 doublereal lplane[4];
35 extern /* Subroutine */ int inedpl_(doublereal *, doublereal *,
36 doublereal *, doublereal *, doublereal *, logical *);
37 doublereal normal[3];
38 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
39 ftnlen), setmsg_(char *, ftnlen);
40 extern logical return_(void);
41
42 /* $ Abstract */
43
44 /* Find the limb of a triaxial ellipsoid, viewed from a specified */
45 /* point. */
46
47 /* $ Disclaimer */
48
49 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
50 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
51 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
52 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
53 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
54 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
55 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
56 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
57 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
58 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
59
60 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
61 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
62 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
63 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
64 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
65 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
66
67 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
68 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
69 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
70 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
71
72 /* $ Required_Reading */
73
74 /* ELLIPSES */
75
76 /* $ Keywords */
77
78 /* ELLIPSE */
79 /* ELLIPSOID */
80 /* GEOMETRY */
81 /* MATH */
82
83 /* $ Declarations */
84 /* $ Brief_I/O */
85
86 /* Variable I/O Description */
87 /* -------- --- -------------------------------------------------- */
88 /* A I Length of ellipsoid semi-axis lying on the x-axis. */
89 /* B I Length of ellipsoid semi-axis lying on the y-axis. */
90 /* C I Length of ellipsoid semi-axis lying on the z-axis. */
91 /* VIEWPT I Location of viewing point. */
92 /* LIMB O Limb of ellipsoid as seen from viewing point. */
93
94 /* $ Detailed_Input */
95
96 /* A, */
97 /* B, */
98 /* C are the lengths of the semi-axes of a triaxial */
99 /* ellipsoid. The ellipsoid is centered at the */
100 /* origin and oriented so that its axes lie on the */
101 /* x, y and z axes. A, B, and C are the lengths of */
102 /* the semi-axes that point in the x, y, and z */
103 /* directions respectively. */
104
105 /* VIEWPT is a point from which the ellipsoid is viewed. */
106 /* VIEWPT must be outside of the ellipsoid. */
107
108 /* $ Detailed_Output */
109
110 /* LIMB is a SPICELIB ellipse that represents the limb of */
111 /* the ellipsoid. */
112
113 /* $ Parameters */
114
115 /* None. */
116
117 /* $ Exceptions */
118
119 /* 1) If the length of any semi-axis of the ellipsoid is */
120 /* non-positive, the error SPICE(INVALIDAXISLENGTH) is signalled. */
121 /* LIMB is not modified. */
122
123 /* 2) If the length of any semi-axis of the ellipsoid is zero after */
124 /* the semi-axis lengths are scaled by the reciprocal of the */
125 /* magnitude of the longest semi-axis and then squared, the error */
126 /* SPICE(DEGENERATECASE) is signalled. LIMB is not modified. */
127
128 /* 3) If the viewing point VIEWPT is inside the ellipse, the error */
129 /* SPICE(INVALIDPOINT) is signalled. LIMB is not modified. */
130
131 /* 4) If the geometry defined by the input ellipsoid and viewing */
132 /* point is so extreme that the limb cannot be found, the error */
133 /* SPICE(DEGENERATECASE) is signalled. */
134
135 /* 5) If the shape of the ellipsoid and the viewing geometry are */
136 /* such that the limb is an excessively flat ellipsoid, the */
137 /* limb may be a degenerate ellipse. You must determine whether */
138 /* this possibility poses a problem for your application. */
139
140 /* $ Files */
141
142 /* None. */
143
144 /* $ Particulars */
145
146 /* The limb of a body, as seen from a viewing point, is the boundary */
147 /* of the portion of the body's surface that is visible from that */
148 /* viewing point. In this definition, we consider a surface point */
149 /* to be `visible' if it can be connected to the viewing point by a */
150 /* line segment that doen't pass through the body. This is a purely */
151 /* geometrical definition that ignores the matter of which portions */
152 /* of the surface are illuminated, or whether the view is obscured by */
153 /* any additional objects. */
154
155 /* If a body is modelled as a triaxial ellipsoid, the limb is always */
156 /* an ellipse. The limb is determined by its center, a semi-major */
157 /* axis vector, and a semi-minor axis vector. */
158
159 /* We note that the problem of finding the limb of a triaxial */
160 /* ellipsoid is mathematically identical to that of finding its */
161 /* terminator, if one makes the simplifying assumption that the */
162 /* terminator is the limb of the body as seen from the vertex of the */
163 /* umbra. So, this routine can be used to solve this simplified */
164 /* version of the problem of finding the terminator. */
165
166 /* $ Examples */
167
168 /* 1) We'd like to find the apparent limb of Jupiter, corrected for */
169 /* light time, as seen from a spacecraft's position at time ET. */
170
171 /* C */
172 /* C Find the viewing point in Jupiter-centered */
173 /* C coordinates. To do this, find the apparent position */
174 /* C of Jupiter as seen from the spacecraft and negate */
175 /* C this vector. In this case we'll use light time */
176 /* C correction to arrive at the apparent limb. JSTAT is */
177 /* C the Jupiter's state (position and velocity) as seen */
178 /* C from the spacecraft. SCPOS is the spacecraft's */
179 /* C position relative to Jupiter. */
180 /* C */
181 /* CALL SPKEZ ( JUPID, ET, 'J2000', 'LT', SCID, */
182 /* . SCSTAT, LT ) */
183
184 /* CALL VMINUS ( SCSTAT, SCPOS ) */
185
186 /* C */
187 /* C Get Jupiter's semi-axis lengths... */
188 /* C */
189 /* CALL BODVCD ( JUPID, 'RADII', 3, N, RAD ) */
190
191 /* C */
192 /* C ...and the transformation from J2000 to Jupiter */
193 /* C equator and prime meridian coordinates. Note that we */
194 /* C use the orientation of Jupiter at the time of */
195 /* C emission of the light that arrived at the */
196 /* C spacecraft at time ET. */
197 /* C */
198 /* CALL BODMAT ( JUPID, ET-LT, TIPM ) */
199
200 /* C */
201 /* C Transform the spacecraft's position into Jupiter- */
202 /* C fixed coordinates. */
203 /* C */
204 /* CALL MXV ( TIPM, SCPOS, SCPOS ) */
205
206 /* C */
207 /* C Find the apparent limb. LIMB is a SPICELIB ellipse */
208 /* C representing the limb. */
209 /* C */
210 /* CALL EDLIMB ( RAD(1), RAD(2), RAD(3), SCPOS, LIMB ) */
211
212 /* C */
213 /* C LCENTR, SMAJOR, and SMINOR are the limb's center, */
214 /* C semi-major axis of the limb, and a semi-minor axis */
215 /* C of the limb. We obtain these from LIMB using the */
216 /* C SPICELIB routine EL2CGV ( Ellipse to center and */
217 /* C generating vectors ). */
218 /* C */
219 /* CALL EL2CGV ( LIMB, CENTER, SMAJOR, SMINOR ) */
220
221 /* $ Restrictions */
222
223 /* None. */
224
225 /* $ Literature_References */
226
227 /* None. */
228
229 /* $ Author_and_Institution */
230
231 /* N.J. Bachman (JPL) */
232
233 /* $ Version */
234
235 /* - SPICELIB Version 1.3.0, 23-OCT-2005 (NJB) */
236
237 /* Updated to remove non-standard use of duplicate arguments */
238 /* in VSCLG call. Updated header to refer to BODVCD instead */
239 /* of BODVAR. */
240
241 /* - SPICELIB Version 1.2.0, 06-OCT-1993 (NJB) */
242
243 /* Declaration of unused local variable NEAR was removed. */
244
245 /* - SPICELIB Version 1.1.1, 10-MAR-1992 (WLT) */
246
247 /* Comment section for permuted index source lines was added */
248 /* following the header. */
249
250 /* - SPICELIB Version 1.1.0, 04-DEC-1990 (NJB) */
251
252 /* Error message and description changed for non-positive */
253 /* axis length error. */
254
255 /* - SPICELIB Version 1.0.0, 02-NOV-1990 (NJB) */
256
257 /* -& */
258 /* $ Index_Entries */
259
260 /* ellipsoid limb */
261
262 /* -& */
263 /* $ Revisions */
264
265 /* - SPICELIB Version 1.3.0, 23-OCT-2005 (NJB) */
266
267 /* Updated to remove non-standard use of duplicate arguments */
268 /* in VSCLG call. Updated header to refer to BODVCD instead */
269 /* of BODVAR. */
270
271 /* - SPICELIB Version 1.2.0, 06-OCT-1993 (NJB) */
272
273 /* Declaration of unused local variable NEAR was removed. */
274
275 /* - SPICELIB Version 1.1.0, 04-DEC-1990 (NJB) */
276
277 /* Error message and description changed for non-positive */
278 /* axis length error. The former message and description did */
279 /* not match, and the description was incorrect: it described */
280 /* `zero-length', rather than `non-positive' axes as invalid. */
281
282 /* -& */
283
284 /* SPICELIB functions */
285
286
287 /* Local parameters */
288
289
290 /* Local variables */
291
292
293 /* Standard SPICE error handling. */
294
295 if (return_()) {
296 return 0;
297 } else {
298 chkin_("EDLIMB", (ftnlen)6);
299 }
300
301 /* The semi-axes must have positive length. */
302
303 if (*a <= 0. || *b <= 0. || *c__ <= 0.) {
304 setmsg_("Semi-axis lengths: A = #, B = #, C = #. ", (ftnlen)41);
305 errdp_("#", a, (ftnlen)1);
306 errdp_("#", b, (ftnlen)1);
307 errdp_("#", c__, (ftnlen)1);
308 sigerr_("SPICE(INVALIDAXISLENGTH)", (ftnlen)24);
309 chkout_("EDLIMB", (ftnlen)6);
310 return 0;
311 }
312
313 /* Scale the semi-axes lengths for better numerical behavior. */
314 /* If squaring any one of the scaled lengths causes it to */
315 /* underflow to zero, we cannot continue the computation. Otherwise, */
316 /* scale the viewing point too. */
317
318 /* Computing MAX */
319 d__1 = abs(*a), d__2 = abs(*b), d__1 = max(d__1,d__2), d__2 = abs(*c__);
320 scale = max(d__1,d__2);
321 scla = *a / scale;
322 sclb = *b / scale;
323 sclc = *c__ / scale;
324 /* Computing 2nd power */
325 d__1 = scla;
326 scla2 = d__1 * d__1;
327 /* Computing 2nd power */
328 d__1 = sclb;
329 sclb2 = d__1 * d__1;
330 /* Computing 2nd power */
331 d__1 = sclc;
332 sclc2 = d__1 * d__1;
333 if (scla2 == 0. || sclb2 == 0. || sclc2 == 0.) {
334 setmsg_("Semi-axis too small: A = #, B = #, C = #. ", (ftnlen)43);
335 errdp_("#", a, (ftnlen)1);
336 errdp_("#", b, (ftnlen)1);
337 errdp_("#", c__, (ftnlen)1);
338 sigerr_("SPICE(DEGENERATECASE)", (ftnlen)21);
339 chkout_("EDLIMB", (ftnlen)6);
340 return 0;
341 }
342 d__1 = 1. / scale;
343 vscl_(&d__1, viewpt, v);
344
345 /* The viewing point must be outside of the ellipsoid. LEVEL is the */
346 /* constant of the level surface that V lies on. The ellipsoid */
347 /* itself is the level surface corresponding to LEVEL = 1. */
348
349 /* Computing 2nd power */
350 d__1 = v[0];
351 /* Computing 2nd power */
352 d__2 = v[1];
353 /* Computing 2nd power */
354 d__3 = v[2];
355 level = d__1 * d__1 / scla2 + d__2 * d__2 / sclb2 + d__3 * d__3 / sclc2;
356 if (level < 1.) {
357 setmsg_("Viewing point is inside the ellipsoid.", (ftnlen)38);
358 sigerr_("SPICE(DEGENERATECASE)", (ftnlen)21);
359 chkout_("EDLIMB", (ftnlen)6);
360 return 0;
361 }
362
363 /* Find a normal vector for the limb plane. */
364
365 /* To compute this vector, we use the fact that the surface normal at */
366 /* each limb point is orthogonal to the line segment connecting the */
367 /* viewing point and the limb point. Let the notation */
368
369 /* < a, b > */
370
371 /* indicate the dot product of the vectors a and b. If we call the */
372 /* viewing point V and the limb point X, then */
373
374
375
376 /* X(1) X(2) X(3) */
377 /* 0 = < V - X, ( -------- , -------- , -------- ) > */
378 /* 2 2 2 */
379 /* SCLA SCLB SCLC */
380
381
382 /* X(1) X(2) X(3) */
383 /* = < V, ( -------- , -------- , -------- ) > */
384 /* 2 2 2 */
385 /* SCLA SCLB SCLC */
386
387
388 /* X(1) X(2) X(3) */
389 /* - < X, ( -------- , -------- , -------- ) > */
390 /* 2 2 2 */
391 /* SCLA SCLB SCLC */
392
393 /* 2 2 2 */
394 /* X(1) X(2) X(3) */
395 /* = -------- + -------- + -------- */
396 /* 2 2 2 */
397 /* SCLA SCLB SCLC */
398
399
400 /* = 1 */
401
402
403 /* This last equation is just the equation of the scaled ellipsoid. */
404 /* We can combine the last two equalities and interchange the */
405 /* positions of X and V to obtain */
406
407
408 /* V(1) V(2) V(3) */
409 /* < X, ( -------- , -------- , -------- ) > = 1 */
410 /* 2 2 2 */
411 /* SCLA SCLB SCLC */
412
413
414 /* This is the equation of the limb plane. */
415
416
417 /* Put together a SPICELIB plane, LPLANE, that represents the limb */
418 /* plane. */
419
420 normal[0] = v[0] / scla2;
421 normal[1] = v[1] / sclb2;
422 normal[2] = v[2] / sclc2;
423 nvc2pl_(normal, &c_b18, lplane);
424
425 /* Find the limb by intersecting the limb plane with the ellipsoid. */
426
427 inedpl_(&scla, &sclb, &sclc, lplane, limb, &found);
428
429 /* FOUND should be true unless we've encountered numerical problems. */
430
431 if (! found) {
432 setmsg_("Ellipsoid shape and viewing geometry are too extreme; the l"
433 "imb was not found. ", (ftnlen)78);
434 sigerr_("SPICE(DEGENERATECASE)", (ftnlen)21);
435 chkout_("EDLIMB", (ftnlen)6);
436 return 0;
437 }
438
439 /* Undo the scaling before returning the limb. */
440
441 vsclg_(&scale, limb, &c__9, tmpel);
442 moved_(tmpel, &c__9, limb);
443 chkout_("EDLIMB", (ftnlen)6);
444 return 0;
445 } /* edlimb_ */
446
447