1 /* nearpt.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 integer c__3 = 3;
11 static doublereal c_b36 = 2.;
12 static integer c__2048 = 2048;
13 static doublereal c_b108 = 1e-16;
14
15 /* $Procedure NEARPT ( Nearest point on an ellipsoid ) */
nearpt_(doublereal * positn,doublereal * a,doublereal * b,doublereal * c__,doublereal * npoint,doublereal * alt)16 /* Subroutine */ int nearpt_(doublereal *positn, doublereal *a, doublereal *b,
17 doublereal *c__, doublereal *npoint, doublereal *alt)
18 {
19 /* Initialized data */
20
21 static char mssg[80*7] = "Axis A was nonpositive. ? "
22 " " "Axis B was nonpositive. ? "
23 " " "Axes A a"
24 "nd B were nonpositive. ? "
25 " " "Axis C was nonpositive. ? "
26 " " "Axes A and C were nonpositive. ? "
27 " " "Axes B and C we"
28 "re nonpositive. ? "
29 "All three axes were nonpositive. ? "
30 " ";
31
32 /* System generated locals */
33 integer i__1, i__2, i__3, i__4, i__5, i__6;
34 doublereal d__1, d__2, d__3, d__4, d__5;
35
36 /* Builtin functions */
37 integer s_rnge(char *, integer, char *, integer);
38 double sqrt(doublereal), pow_dd(doublereal *, doublereal *);
39
40 /* Local variables */
41 extern /* Subroutine */ int vadd_(doublereal *, doublereal *, doublereal *
42 );
43 doublereal sign, axis[3], temp, term[3], errp[3], copy[3];
44 logical trim;
45 extern /* Subroutine */ int vequ_(doublereal *, doublereal *);
46 integer i__;
47 doublereal q, scale;
48 extern /* Subroutine */ int chkin_(char *, ftnlen);
49 doublereal denom;
50 extern /* Subroutine */ int errch_(char *, char *, ftnlen, ftnlen);
51 extern doublereal dpmax_(void);
52 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen);
53 logical extra;
54 doublereal lower;
55 extern doublereal vdist_(doublereal *, doublereal *);
56 doublereal point[3], pnorm, upper;
57 extern /* Subroutine */ int vperp_(doublereal *, doublereal *, doublereal
58 *);
59 extern doublereal vnorm_(doublereal *);
60 doublereal denom2, denom3, lambda, tlambd[3], height;
61 extern doublereal brcktd_(doublereal *, doublereal *, doublereal *);
62 logical inside;
63 doublereal factor;
64 extern /* Subroutine */ int orderd_(doublereal *, integer *, integer *),
65 reordd_(integer *, integer *, doublereal *);
66 doublereal toobig;
67 integer iorder[3];
68 extern doublereal touchd_(doublereal *);
69 doublereal olderr, normal[3], bestht, orignl[3], prodct;
70 extern /* Subroutine */ int sigerr_(char *, ftnlen);
71 doublereal epoint[3];
72 extern /* Subroutine */ int chkout_(char *, ftnlen), vsclip_(doublereal *,
73 doublereal *);
74 doublereal bestpt[3], newerr;
75 extern /* Subroutine */ int setmsg_(char *, ftnlen);
76 doublereal axisqr[3];
77 extern logical approx_(doublereal *, doublereal *, doublereal *);
78 doublereal qlower;
79 integer snglpt;
80 doublereal qupper, spoint[3];
81 extern logical return_(void);
82 logical solvng;
83 extern /* Subroutine */ int errint_(char *, integer *, ftnlen), surfnm_(
84 doublereal *, doublereal *, doublereal *, doublereal *,
85 doublereal *);
86 integer solutn, bad;
87 doublereal err[3];
88 integer itr;
89
90 /* $ Abstract */
91
92 /* This routine locates the point on the surface of an ellipsoid */
93 /* that is nearest to a specified position. It also returns the */
94 /* altitude of the position above the ellipsoid. */
95
96 /* $ Disclaimer */
97
98 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
99 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
100 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
101 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
102 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
103 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
104 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
105 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
106 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
107 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
108
109 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
110 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
111 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
112 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
113 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
114 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
115
116 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
117 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
118 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
119 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
120
121 /* $ Required_Reading */
122
123 /* None. */
124
125 /* $ Keywords */
126
127 /* ALTITUDE */
128 /* ELLIPSOID */
129 /* GEOMETRY */
130
131 /* $ Declarations */
132 /* $ Brief_I/O */
133
134 /* VARIABLE I/O DESCRIPTION */
135 /* -------- --- -------------------------------------------------- */
136 /* POSITN I Position of a point in body-fixed frame. */
137 /* A I Length of semi-axis parallel to x-axis. */
138 /* B I Length of semi-axis parallel to y-axis. */
139 /* C I Length on semi-axis parallel to z-axis. */
140 /* NPOINT O Point on the ellipsoid closest to POSITN. */
141 /* ALT O Altitude of POSITN above the ellipsoid. */
142
143 /* $ Detailed_Input */
144
145 /* POSITN 3-vector giving the position of a point with respect */
146 /* to the center of an ellipsoid. The vector is expressed */
147 /* in a body-fixed reference frame. The semi-axes of the */
148 /* ellipsoid are aligned with the x, y, and z-axes of the */
149 /* body-fixed frame. */
150
151 /* A Length of the semi-axis of the ellipsoid that is */
152 /* parallel to the x-axis of the body-fixed reference */
153 /* frame. */
154
155 /* B Length of the semi-axis of the ellipsoid that is */
156 /* parallel to the y-axis of the body-fixed reference */
157 /* frame. */
158
159 /* C Length of the semi-axis of the ellipsoid that is */
160 /* parallel to the z-axis of the body-fixed reference */
161 /* frame. */
162
163 /* $ Detailed_Output */
164
165 /* NPOINT is the nearest point on the ellipsoid to POSITN. */
166 /* NPOINT is a 3-vector expressed in the body-fixed */
167 /* reference frame. */
168
169 /* ALT is the altitude of POSITN above the ellipsoid. If */
170 /* POSITN is inside the ellipsoid, ALT will be negative */
171 /* and have magnitude equal to the distance between */
172 /* NPOINT and POSITN. */
173
174 /* $ Parameters */
175
176 /* None. */
177
178 /* $ Exceptions */
179
180 /* 1) If any of the axis lengths A, B or C are non-positive, the */
181 /* error SPICE(BADAXISLENGTH) will be signaled. */
182
183 /* 2) If the ratio of the longest to the shortest ellipsoid axis */
184 /* is large enough so that arithmetic expressions involving its */
185 /* squared value may overflow, the error SPICE(BADAXISLENGTH) */
186 /* will be signaled. */
187
188 /* 3) If any of the expressions */
189
190 /* A * ABS( POSITN(1) ) / m**2 */
191 /* B * ABS( POSITN(2) ) / m**2 */
192 /* C * ABS( POSITN(3) ) / m**2 */
193
194 /* where m is the minimum of { A, B, C }, is large enough so */
195 /* that arithmetic expressions involving these sub-expressions */
196 /* may overflow, the error SPICE(INPUTSTOOLARGE) is signaled. */
197
198 /* 4) If the axes of the ellipsoid have radically different */
199 /* magnitudes, for example if the ratios of the axis lengths vary */
200 /* by 10 orders of magnitude, the results may have poor */
201 /* precision. No error checks are done to identify this problem. */
202
203 /* 5) If the axes of the ellipsoid and the input point POSITN have */
204 /* radically different magnitudes, for example if the ratio of */
205 /* the magnitude of POSITN to the length of the shortest axis is */
206 /* 1.E25, the results may have poor precision. No error checks */
207 /* are done to identify this problem. */
208
209 /* $ Files */
210
211 /* None. */
212
213 /* $ Particulars */
214
215 /* Many applications of this routine are more easily performed */
216 /* using the higher-level SPICELIB routine SUBPNT. This routine */
217 /* is the mathematical workhorse on which SUBPNT relies. */
218
219 /* $ Examples */
220
221 /* Example 1. */
222
223 /* The code fragment below illustrates how you can use SPICELIB to */
224 /* compute the apparent sub-earth point on the moon. */
225
226 /* C */
227 /* C Load the ephemeris, leapseconds and physical constants */
228 /* C files first. We assume the names of these files are */
229 /* C stored in the character variables SPK, LSK and */
230 /* C PCK. */
231 /* C */
232 /* CALL FURNSH ( SPK ) */
233 /* CALL FURNSH ( LSK ) */
234 /* CALL FURNSH ( PCK ) */
235
236 /* C */
237 /* C Get the apparent position of the moon as seen from the */
238 /* C earth. Look up this position vector in the moon */
239 /* C body-fixed frame IAU_MOON. The orientation of the */
240 /* C IAU_MOON frame will be computed at epoch ET-LT. */
241 /* C */
242 /* CALL SPKPOS ( 'moon', ET, 'IAU_MOON', 'LT+S', */
243 /* . 'earth', TRGPOS, LT ) */
244
245 /* C */
246 /* C Negate the moon's apparent position to obtain the */
247 /* C position of the earth in the moon's body-fixed frame. */
248 /* C */
249 /* CALL VMINUS ( TRGPOS, EVEC ) */
250
251 /* C */
252 /* C Get the lengths of the principal axes of the moon. */
253 /* C Transfer the elements of the array RADII to the */
254 /* C variables A, B, C to enhance readability. */
255 /* C */
256 /* CALL BODVRD ( 'MOON', 'RADII', DIM, RADII ) */
257 /* CALL VUPACK ( RADII, A, B, C ) */
258
259 /* C */
260 /* C Finally get the point SUBPNT on the surface of the */
261 /* C moon closest to the earth --- the sub-earth point. */
262 /* C SUBPNT is expressed in the IAU_MOON reference frame. */
263 /* C */
264 /* CALL NEARPT ( EVEC, A, B, C, SUBPNT, ALT ) */
265
266
267 /* Example 2. */
268
269 /* One can use this routine to define a generalization of GEODETIC */
270 /* coordinates called GAUSSIAN coordinates of a triaxial body. (The */
271 /* name is derived from the famous Gauss-map of classical */
272 /* differential geometry). The coordinates are longitude, */
273 /* latitude, and altitude. */
274
275 /* We let the x-axis of the body fixed coordinate system point */
276 /* along the longest axis of the triaxial body. The y-axis points */
277 /* along the middle axis and the z-axis points along the shortest */
278 /* axis. */
279
280 /* Given a point P, there is a point on the ellipsoid that is */
281 /* closest to P, call it Q. The latitude and longitude of P */
282 /* are determined by constructing the outward pointing unit normal */
283 /* to the ellipsoid at Q. Latitude of P is the latitude that the */
284 /* normal points toward in the body-fixed frame. Longitude is */
285 /* the longitude the normal points to in the body-fixed frame. */
286 /* The altitude is the signed distance from P to Q, positive if P */
287 /* is outside the ellipsoid, negative if P is inside. */
288 /* (the mapping of the point Q to the unit normal at Q is the */
289 /* Gauss-map of Q). */
290
291 /* To obtain the Gaussian coordinates of a point whose position */
292 /* in body-fixed rectangular coordinates is given by a vector P, */
293 /* the code fragment below will suffice. */
294
295 /* CALL NEARPT ( P, A, B, C, Q, ALT ) */
296 /* CALL SURFNM ( A, B, C Q, NRML ) */
297
298 /* CALL RECLAT ( NRML, R, LONG, LAT ) */
299
300 /* The Gaussian coordinates are LONG, LAT, and ALT. */
301
302
303 /* $ Restrictions */
304
305 /* See the Exceptions header section above. */
306
307 /* $ Literature_References */
308
309 /* None. */
310
311 /* $ Author_and_Institution */
312
313 /* C.H. Acton (JPL) */
314 /* N.J. Bachman (JPL) */
315 /* W.L. Taber (JPL) */
316 /* E.D. Wright (JPL) */
317
318 /* $ Version */
319
320 /* - SPICELIB Version 1.4.0, 27-JUN-2013 (NJB) */
321
322 /* Updated in-line comments. */
323
324 /* Last update was 04-MAR-2013 (NJB) */
325
326 /* Bug fix: now correctly computes off-axis solution for */
327 /* the case of a prolate ellipsoid and a viewing point */
328 /* on the interior long axis. */
329
330 /* - SPICELIB Version 1.3.1, 07-FEB-2008 (NJB) */
331
332 /* Header update: header now refers to SUBPNT rather */
333 /* than deprecated routine SUBPT. */
334
335 /* - SPICELIB Version 1.3.0, 07-AUG-2006 (NJB) */
336
337 /* Bug fix: added initialization of variable SNGLPT to support */
338 /* operation under the Macintosh Intel Fortran */
339 /* compiler. Note that this bug did not affect */
340 /* operation of this routine on other platforms. */
341
342 /* - SPICELIB Version 1.2.0, 15-NOV-2005 (EDW) (NJB) */
343
344 /* Various changes were made to ensure that all loops terminate. */
345
346 /* Bug fix: scale of transverse component of error vector */
347 /* was corrected for the exterior point case. */
348
349 /* Bug fix: non-standard use of duplicate arguments in VSCL */
350 /* calls was corrected. */
351
352 /* Error checking was added to screen out inputs that might */
353 /* cause numeric overflow. */
354
355 /* Replaced BODVAR call in examples to BODVRD. */
356
357 /* - SPICELIB Version 1.1.1, 28-JUL-2003 (NJB) (CHA) */
358
359 /* Various header corrections were made. */
360
361 /* - SPICELIB Version 1.1.0, 27-NOV-1990 (WLT) */
362
363 /* The routine was substantially rewritten to achieve */
364 /* more robust behavior and document the mathematics */
365 /* of the routine. */
366
367 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
368
369 /* Comment section for permuted index source lines was added */
370 /* following the header. */
371
372 /* - SPICELIB Version 1.0.0, 31-JAN-1990 (WLT) */
373
374 /* -& */
375 /* $ Index_Entries */
376
377 /* distance from point to ellipsoid */
378 /* nearest point on an ellipsoid */
379
380 /* -& */
381 /* $ Revisions */
382
383 /* - SPICELIB Version 1.1.0, 07-AUG-2006 (NJB) */
384
385 /* Bug fix: added initialization of variable SNGLPT to support */
386 /* operation under the Macintosh Intel Fortran */
387 /* compiler. Note that this bug did not affect */
388 /* operation of this routine on other platforms. The */
389 /* statement referencing the uninitialized variable */
390 /* was: */
391
392 /* IF ( INSIDE .AND. ( SNGLPT .EQ. 2 */
393 /* . .OR. SNGLPT .EQ. 3 ) ) THEN */
394
395 /* SNGLPT is uninitialized only if INSIDE is .FALSE., */
396 /* so the value of the logical expression is not affected by */
397 /* the uninitialized value of SNGLPT. */
398
399 /* However, the Intel Fortran compiler for the Mac flags a runtime */
400 /* error when the above code is exercised. So SNGLPT is now */
401 /* initialized prior to the above IF statement. */
402
403
404 /* - SPICELIB Version 1.2.0, 15-NOV-2005 (EDW) (NJB) */
405
406 /* Bug fix: scale of transverse component of error vector */
407 /* was corrected for the exterior point case. */
408 /* Replaced BODVAR call in examples to BODVRD. */
409
410 /* Bug fix: non-standard use of duplicate arguments in VSCL */
411 /* calls was corrected. */
412
413 /* Various changes were made to ensure that all loops terminate. */
414
415 /* Error checking was added to screen out inputs that might */
416 /* cause numeric overflow. */
417
418 /* Removed secant solution branch from root-finding loop. */
419 /* Although the secant solution sped up some root searches, */
420 /* it caused large numbers of unnecessary iterations in others. */
421
422 /* Changed the expression: */
423
424 /* IF ( LAMBDA .EQ. LOWER */
425 /* . .OR. LAMBDA .EQ. UPPER ) THEN */
426
427 /* to */
428
429 /* IF ( APPROX( LAMBDA, LOWER, CNVTOL ) */
430 /* . .OR. APPROX( LAMBDA, UPPER, CNVTOL ) ) THEN */
431
432 /* Use of APPROX eliminates the possibility of an infinite loop */
433 /* when LAMBDA approaches to within epsilon of, but does not */
434 /* equate to UPPER or LOWER. Infinite loops occurred under some */
435 /* compiler's optimizations. */
436
437 /* The loop also includes a check on number of iterations, */
438 /* signaling an error if the bisection loop uses more than */
439 /* MAXITR passes. */
440
441 /* TOUCHD is now used to defeat extended-register usage in */
442 /* cases where such usage may cause logic problems. */
443
444 /* Some minor code changes were made to ensure that various */
445 /* variables remain in their expected ranges. */
446
447 /* A few code changes were made to enhance clarity. */
448
449
450 /* - SPICELIB Version 1.1.0, 27-NOV-1990 */
451
452 /* The routine was nearly rewritten so that points */
453 /* near the coordinate planes in the interior of the ellipsoid */
454 /* could be handled without fear of floating point overflow */
455 /* or divide by zero. */
456
457 /* While the mathematical ideas involved in the original routine */
458 /* are retained, the code is for the most part new. In addition, */
459 /* the new code has been documented far more extensively than was */
460 /* NEARPT 1.0.0. */
461
462
463 /* - Beta Version 2.0.0, 9-JAN-1989 (WLT) */
464
465 /* Error handling added has been added for bad axes values. */
466
467 /* The algorithm did not work correctly for some points inside */
468 /* the ellipsoid lying on the plane orthogonal to the shortest */
469 /* axis of the ellipsoid. The problem was corrected. */
470
471 /* Finally the algorithm was made slightly more robust and clearer */
472 /* by use of SPICELIB routines and by normalizing the inputs. */
473
474 /* Add an example to the header section. */
475 /* -& */
476
477 /* SPICELIB functions */
478
479
480 /* Parameters */
481
482
483 /* The convergence tolerance CNVTOL is used to terminate the */
484 /* bisection loop when the solution interval is very small but */
485 /* hasn't converged to length zero. This situation can occur when */
486 /* the root is extremely close to zero. */
487
488
489 /* Various potentially large numbers we'll compute must not */
490 /* exceed DPMAX()/MARGIN: */
491
492
493 /* The parameter MAXSOL determines the maximum number of */
494 /* iterations that will be performed in locating the */
495 /* near point. This must be at least 3. To get strong */
496 /* robustness in the routine, MAXSOL should be at least 4. */
497
498
499 /* MAXITR defines the maximum number of iterations allowed in */
500 /* the bisection loop used to find LAMBDA. If this loop requires */
501 /* more than MAXITR iterations to achieve convergence, NEARPT */
502 /* will signal an error. */
503
504 /* On a PC/Linux/g77 platform, it has been observed that each */
505 /* bisection loop normally completes in fewer than 70 iterations. */
506 /* MAXITR is used as a "backstop" to prevent infinite looping in */
507 /* case the normal loop termination conditions don't take effect. */
508 /* The value selected is based on the range of exponents for IEEE */
509 /* double precision floating point numbers. */
510
511
512 /* Length of lines in message buffer. */
513
514
515 /* Local Variables */
516
517
518 /* Saved variables */
519
520
521 /* Initial values */
522
523
524 /* Here's what you can expect to find in the routine below. */
525
526 /* Chapter 1. Error and Exception Handling. */
527
528 /* Chapter 2. Mathematical background for the solution---the */
529 /* lambda equation. */
530
531 /* Chapter 3. Initializations for the main processing loop. */
532
533 /* Chapter 4. Mathematical Solution of the lambda equation. */
534
535 /* Section 4.1 Avoiding numerical difficulties. */
536 /* Section 4.2 Bracketing the root of the lambda */
537 /* equation. */
538 /* Section 4.3 Refining the estimate of lambda. */
539 /* Section 4.4 Handling points on the central plane. */
540
541 /* Chapter 5. Decisions and initializations for sharpening */
542 /* the solution. */
543
544 /* Chapter 6. Clean up. */
545
546
547 /* Error and Exception Handling. */
548 /* ================================================================ */
549 /* ---------------------------------------------------------------- */
550
551 if (return_()) {
552 return 0;
553 } else {
554 chkin_("NEARPT", (ftnlen)6);
555 }
556
557 /* Check the axes to make sure that none of them is less than or */
558 /* equal to zero. If one is, signal an error and return. */
559
560 bad = 0;
561 if (*a <= 0.) {
562 ++bad;
563 }
564 if (*b <= 0.) {
565 bad += 2;
566 }
567 if (*c__ <= 0.) {
568 bad += 4;
569 }
570 if (bad > 0) {
571 setmsg_(mssg + ((i__1 = bad - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge(
572 "mssg", i__1, "nearpt_", (ftnlen)591)) * 80, (ftnlen)80);
573 errch_("?", "The A,B, and C axes were #, #, and # respectively.", (
574 ftnlen)1, (ftnlen)50);
575 errdp_("#", a, (ftnlen)1);
576 errdp_("#", b, (ftnlen)1);
577 errdp_("#", c__, (ftnlen)1);
578 sigerr_("SPICE(BADAXISLENGTH)", (ftnlen)20);
579 chkout_("NEARPT", (ftnlen)6);
580 return 0;
581 }
582
583 /* Mathematical background for the solution---the lambda equation. */
584 /* ================================================================ */
585 /* ---------------------------------------------------------------- */
586
587
588 /* Here is the background and general outline of how this problem is */
589 /* going to be solved. */
590
591 /* We want to find a point on the ellipsoid */
592
593
594 /* X**2 Y**2 Z**2 */
595 /* ------ + ------ + ------ = 1 */
596 /* A**2 B**2 C**2 */
597
598 /* that is closest to the input point POSITN. */
599
600 /* If one cares about the gory details, we know that such a */
601 /* point must exist because the ellipsoid is a compact subset of */
602 /* Euclidean 3-space and the distance function between the input */
603 /* point and the ellipsoid is continuous. Since any continuous */
604 /* function on a compact set actually achieves its minimum at */
605 /* some point of the compact set, we are guaranteed that a */
606 /* closest point exists. The closest point may not be unique if */
607 /* the input point is inside the ellipsoid. */
608
609 /* If we let NPOINT be a closest point to POSITN, then the line */
610 /* segment joining POSITN to NPOINT is parallel to the normal to the */
611 /* ellipsoid at NPOINT. Moreover, suppose we let SEGMENT(P) be the */
612 /* line segment that connects an arbitrary point P with POSITN. It */
613 /* can be shown that there is only one point P on the ellipsoid in */
614 /* the same octant at POSITN such that the normal at P is parallel */
615 /* to SEGMENT(P). */
616
617
618 /* More gory details: A normal to a point (X,Y,Z) */
619 /* on the ellipsoid is given by */
620
621 /* (X/A**2, Y/B**2, Z/C**2) */
622
623 /* Given a fixed LAMBDA, and allowing (X,Y,Z) to */
624 /* range over all points on the ellipsoid, the set */
625 /* of points */
626
627
628 /* LAMBDA*X LAMBDA*Y LAMBDA*Z */
629 /* ( X + --------, Y + --------, Z + -------- ) */
630 /* A**2 B**2 C**2 */
631
632 /* describes another ellipsoid with axes having lengths */
633
634 /* LAMBDA LAMBDA LAMBDA */
635 /* A + ------ , B + ------ , C + ------ . */
636 /* A B C */
637
638
639 /* Moreover, as long as LAMBDA > - MIN( A**2, B**2, C**2 ) */
640 /* none of these ellipsoids intersect. Thus, as long as */
641 /* the normal lines are not allowed to cross the coordinate plane */
642 /* orthogonal to the smallest axis (called the central plane) */
643 /* they do not intersect. */
644
645
646 /* Finally every point that does not lie on the central plane */
647 /* lies on one of the "lambda" ellipsoids described above. */
648
649 /* Consequently, for each point, P, not on the central plane */
650 /* there is a unique point, NPOINT, on the ellipsoid, such that */
651 /* the normal line at NPOINT also contains P and does not cross */
652 /* the central plane. */
653
654
655 /* From the above discussion we see that we can mathematically */
656 /* solve the near point problem by finding a point NPOINT */
657 /* on the ellipsoid given by the equation: */
658
659 /* X**2 Y**2 Z**2 */
660 /* ------ + ------ + ------ = 1 */
661 /* A**2 B**2 C**2 */
662
663
664 /* such that for some value of LAMBDA */
665
666 /* POSITN = NPOINT + LAMBDA*NORMAL(NPOINT). */
667
668 /* Moreover, if POSITN = (X_o,Y_o,Z_o) then LAMBDA must satisfy */
669 /* the equation: */
670
671 /* 2 2 2 */
672 /* X_o Y_o Z_o */
673 /* ----------------- + ------------------ + ------------------ = 1 */
674 /* 2 2 2 */
675 /* ( A + LAMBDA/A ) ( B + LAMBDA/B ) ( C + LAMBDA/C ) */
676
677
678 /* and LAMBDA must be greater than -MIN(A**2,B**2,C**2) */
679
680
681 /* Once LAMBDA is known, NPOINT can be computed from the equation: */
682
683 /* POSITN = NPOINT + LAMBDA*NORMAL(NPOINT). */
684
685
686 /* The process of solving for LAMBDA can be viewed as selecting */
687 /* that ellipsoid */
688
689 /* 2 2 2 */
690 /* x y z */
691 /* --------------- + ---------------- + --------------- = 1 */
692 /* 2 2 2 */
693 /* (a + lambda/a) ( b + lambda/b) (c + lambda/c) */
694
695 /* that contains the input point POSITN. For lambda = 0, this */
696 /* ellipsoid is just the input ellipsoid. When we increase */
697 /* lambda we get a larger "inflated" ellipsoid. When we */
698 /* decrease lambda we get a smaller "deflated" ellipsoid. Thus, */
699 /* the search for lambda can be viewed as inflating or deflating */
700 /* the input ellipsoid (in a specially prescribed manner) until */
701 /* the resulting ellipsoid contains the input point POSITN. */
702
703 /* The mathematical solution laid out above, has some numeric */
704 /* flaws. However, it is robust enough so that if it is applied */
705 /* repeatedly, we can converge to a good solution of the near point */
706 /* problem. */
707
708 /* In the code that follows, we will first solve the lambda equation */
709 /* using the original input point. However, for points near the */
710 /* central plane the solution we obtain may not lie on the */
711 /* ellipsoid. But, it should lie along the correct normal line. */
712
713 /* Using this first candidate solution, we find the closest point */
714 /* to it on the ellipsoid. This second iteration always produces */
715 /* a point that is as close as you can get to the ellipsoid. */
716 /* However, the normal at this second solution may not come as close */
717 /* as desired to pointing toward the input position. To overcome */
718 /* this deficiency we sharpen the second solution. */
719
720 /* To sharpen a solution we use the computed near point, the */
721 /* computed altitude of POSITN and the normal at the near point to */
722 /* approximate POSITN. The difference between the approximated */
723 /* position of POSITN and the input value of POSITN is called the */
724 /* error vector. To get a sharpened solution we translate the */
725 /* computed near point in the direction of the component of the */
726 /* error vector orthogonal to the normal. We then find the */
727 /* mathematical near point to our translated solution. */
728
729 /* The sharpening process is repeated until it no longer produces */
730 /* an "improved" near point. */
731
732 /* At each step of this procedure, we must compute a solution to */
733 /* the "lambda" equation in order to produce our next estimate of */
734 /* the near point. If it were possible to create a "private" */
735 /* routine in FORTRAN that only this routine could access, we */
736 /* would do it. However, things being what they are, we have to */
737 /* compute the lambda solution in a loop. We keep track of which */
738 /* refinement we are working on by counting the number of */
739 /* lambda solutions that are computed. */
740
741
742 /* Initializations for the main processing loop */
743 /* ================================================================ */
744 /* ---------------------------------------------------------------- */
745
746
747 /* Let the game begin! */
748
749 /* First order the axes of the ellipsoid and corresponding */
750 /* component of POSITN by the size of lengths of axes. Knowing */
751 /* which axes are smallest will simplify our task of computing */
752 /* lambda when the time comes. */
753
754 axis[0] = *a;
755 axis[1] = *b;
756 axis[2] = *c__;
757 vequ_(positn, point);
758 orderd_(axis, &c__3, iorder);
759 reordd_(iorder, &c__3, axis);
760 reordd_(iorder, &c__3, point);
761
762 /* Rescale everything so as to avoid underflows when squaring */
763 /* quantities and copy the original starting point. */
764
765 /* Be sure that this is UNDONE at the end of the routine. */
766
767 scale = 1. / axis[0];
768 vsclip_(&scale, axis);
769 vsclip_(&scale, point);
770 vequ_(point, orignl);
771
772 /* Save the norm of the scaled input point. */
773
774 pnorm = vnorm_(point);
775
776 /* The scaled axis lengths must be small enough so they can */
777 /* be squared. */
778
779 toobig = sqrt(dpmax_() / 100.);
780
781 /* Note the first axis has length 1.D0, so we don't check it. */
782
783 for (i__ = 2; i__ <= 3; ++i__) {
784 if (axis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("axis",
785 i__1, "nearpt_", (ftnlen)818)] > toobig) {
786 setmsg_("Ratio of length of axis #* to length of axis #* is *; t"
787 "his value may cause numeric overflow.", (ftnlen)92);
788 errint_("*", &iorder[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
789 s_rnge("iorder", i__1, "nearpt_", (ftnlen)823)], (ftnlen)
790 1);
791 errint_("*", iorder, (ftnlen)1);
792 errdp_("*", &axis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
793 s_rnge("axis", i__1, "nearpt_", (ftnlen)825)], (ftnlen)1);
794 sigerr_("SPICE(BADAXISLENGTH)", (ftnlen)20);
795 chkout_("NEARPT", (ftnlen)6);
796 return 0;
797 }
798 }
799
800 /* We also must limit the size of the products */
801
802 /* AXIS(I)*POINT(I), I = 1, 3 */
803
804 /* We can safely check these by comparing the products of */
805 /* the square roots of the factors to TOOBIG. */
806
807 for (i__ = 1; i__ <= 3; ++i__) {
808 prodct = sqrt(axis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
809 "axis", i__1, "nearpt_", (ftnlen)844)]) * sqrt((d__1 = point[(
810 i__2 = i__ - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("point",
811 i__2, "nearpt_", (ftnlen)844)], abs(d__1)));
812 if (prodct > toobig) {
813 setmsg_("Product of length of scaled axis #* and size of corresp"
814 "onding scaled component of POSITN is > *; these values m"
815 "ay cause numeric overflow.", (ftnlen)137);
816 errint_("*", &iorder[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
817 s_rnge("iorder", i__1, "nearpt_", (ftnlen)852)], (ftnlen)
818 1);
819 d__1 = pow_dd(&toobig, &c_b36);
820 errdp_("*", &d__1, (ftnlen)1);
821 sigerr_("SPICE(INPUTSTOOLARGE)", (ftnlen)21);
822 chkout_("NEARPT", (ftnlen)6);
823 return 0;
824 }
825 }
826
827 /* Compute the squared lengths of the scaled axes. */
828
829 axisqr[0] = axis[0] * axis[0];
830 axisqr[1] = axis[1] * axis[1];
831 axisqr[2] = axis[2] * axis[2];
832
833 /* We will need to "solve" for the NEARPT at least 3 times. */
834 /* SOLUTN is the counter that keeps track of how many times */
835 /* we have actually solved for a near point. SOLVNG indicates */
836 /* whether we should continue solving for NEARPT. */
837
838 snglpt = 4;
839 solutn = 1;
840 solvng = TRUE_;
841 while(solvng) {
842
843 /* Mathematical solution of the lambda equation. */
844 /* ================================================================ */
845 /* ---------------------------------------------------------------- */
846
847
848 /* Make a stab at solving the mathematical problem of finding the */
849 /* near point. In other words, solve the lambda equation. */
850
851
852 /* Avoiding Numerical difficulties */
853 /* ------------------------------- */
854
855 /* First make a copy of POINT, then to avoid numerical */
856 /* difficulties later on, we will assume that any component of */
857 /* POINT that is not sufficiently different from zero to */
858 /* contribute to an addition to the corresponding component */
859 /* of AXIS, is in fact zero. */
860
861 vequ_(point, copy);
862 for (i__ = 1; i__ <= 3; ++i__) {
863 if (point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("poi"
864 "nt", i__1, "nearpt_", (ftnlen)903)] * .5 + axis[(i__2 =
865 i__ - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("axis", i__2,
866 "nearpt_", (ftnlen)903)] == axis[(i__3 = i__ - 1) < 3 &&
867 0 <= i__3 ? i__3 : s_rnge("axis", i__3, "nearpt_", (
868 ftnlen)903)] || point[(i__4 = i__ - 1) < 3 && 0 <= i__4 ?
869 i__4 : s_rnge("point", i__4, "nearpt_", (ftnlen)903)] *
870 .5 - axis[(i__5 = i__ - 1) < 3 && 0 <= i__5 ? i__5 :
871 s_rnge("axis", i__5, "nearpt_", (ftnlen)903)] == -axis[(
872 i__6 = i__ - 1) < 3 && 0 <= i__6 ? i__6 : s_rnge("axis",
873 i__6, "nearpt_", (ftnlen)903)]) {
874 point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("poi"
875 "nt", i__1, "nearpt_", (ftnlen)906)] = 0.;
876 }
877 }
878
879 /* OK. Next we set up the logical that indicates whether */
880 /* the current point is inside the ellipsoid. */
881
882 inside = FALSE_;
883
884 /* Bracketing the root of the lambda equation. */
885 /* ------------------------------------------- */
886
887 /* Let (x,y,z) stand for (POINT(1), POINT(2), POINT(3)) and */
888 /* let (a,b,c) stand for (AXIS(1), AXIS(2), AXIS(3)). */
889
890 /* The main step in finding the near point is to find the */
891 /* root of the lambda equation: */
892
893 /* 2 2 2 */
894 /* x y z */
895 /* 0 = --------------- + ---------------- + --------------- - 1 */
896 /* 2 2 2 */
897 /* (a + lambda/a) ( b + lambda/b) (c + lambda/c) */
898
899
900 /* We let Q(lambda) be the right hand side of this equation. */
901 /* To find the roots of the equation we determine */
902 /* values of lambda that make Q greater than 0 and less than 0. */
903 /* An obvious value to check is lambda = 0. */
904
905 /* Computing 2nd power */
906 d__1 = point[0] / axis[0];
907 /* Computing 2nd power */
908 d__2 = point[1] / axis[1];
909 /* Computing 2nd power */
910 d__3 = point[2] / axis[2];
911 q = d__1 * d__1 + d__2 * d__2 + d__3 * d__3 - 1.;
912
913 /* On the first solution pass, we will determine the sign of */
914 /* the altitude of the input POSITN */
915
916 if (solutn == 1) {
917 if (q >= 0.) {
918 sign = 1.;
919 } else {
920 sign = -1.;
921 }
922 }
923
924 /* OK. Now for the stuff we will have to do on every solution */
925 /* pass. */
926
927 /* Below, LOWER and UPPER are the bounds on our independent */
928 /* variable LAMBDA. QLOWER and QUPPER are the values of Q */
929 /* evaluated at LOWER and UPPER, respectively. The root we */
930 /* seek lies in the interval */
931
932 /* [ LOWER, UPPER ] */
933
934 /* At all points in the algorithm, we have, since Q is a */
935 /* decreasing function to the right of the first non-removable */
936 /* singularity, */
937
938 /* QLOWER > 0 */
939 /* - */
940
941 /* QUPPER < 0 */
942 /* - */
943
944 /* We'll use bracketing to ensure that round-off errors don't */
945 /* violate these inequalities. */
946
947 /* The logical flag INSIDE indicates whether the point is */
948 /* strictly inside the interior of the ellipsoid. Points on the */
949 /* surface are not considered to be inside. */
950
951 if (q == 0.) {
952
953 /* In this case the point is already on the ellipsoid */
954 /* (pretty lucky eh?) We simply set our bracketing values, */
955 /* QLOWER and QUPPER, to zero so that that bisection */
956 /* loop won't ever get executed. */
957
958 qlower = 0.;
959 qupper = 0.;
960 lower = 0.;
961 upper = 0.;
962 lambda = 0.;
963 inside = FALSE_;
964 } else if (q > 0.) {
965
966 /* The input point is outside the ellipsoid (we expect that */
967 /* this is the usual case). We want to choose our lower */
968 /* bracketing value so that the bracketing values for lambda */
969 /* aren't too far apart. So we just make sure that the largest */
970 /* term of the expression for Q isn't bigger than 4. */
971
972 for (i__ = 1; i__ <= 3; ++i__) {
973 tlambd[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
974 "tlambd", i__1, "nearpt_", (ftnlen)1011)] = ((d__1 =
975 point[(i__2 = i__ - 1) < 3 && 0 <= i__2 ? i__2 :
976 s_rnge("point", i__2, "nearpt_", (ftnlen)1011)], abs(
977 d__1)) * .5 - axis[(i__3 = i__ - 1) < 3 && 0 <= i__3 ?
978 i__3 : s_rnge("axis", i__3, "nearpt_", (ftnlen)1011)]
979 ) * axis[(i__4 = i__ - 1) < 3 && 0 <= i__4 ? i__4 :
980 s_rnge("axis", i__4, "nearpt_", (ftnlen)1011)];
981 }
982 /* Computing MAX */
983 d__1 = max(0.,tlambd[0]), d__1 = max(d__1,tlambd[1]);
984 lower = max(d__1,tlambd[2]);
985
986 /* Choose the next value of lambda so that the largest term */
987 /* of Q will be no more than 1/4. */
988
989 /* Computing MAX */
990 d__4 = (d__1 = axis[0] * point[0], abs(d__1)), d__5 = (d__2 =
991 axis[1] * point[1], abs(d__2)), d__4 = max(d__4,d__5),
992 d__5 = (d__3 = axis[2] * point[2], abs(d__3));
993 upper = max(d__4,d__5) * 2.;
994 lambda = upper;
995 inside = FALSE_;
996 } else {
997
998 /* In this case the point POSITN is inside the ellipsoid. */
999
1000 inside = TRUE_;
1001
1002 /* This case is a bit of a nuisance. To solve the lambda */
1003 /* equation we have to find upper and lower bounds on */
1004 /* lambda such that one makes Q greater than 0, the other */
1005 /* makes Q less than 0. Once the root has been bracketed */
1006 /* in this way it is a straightforward problem to find */
1007 /* the value of LAMBDA that is closest to the root we */
1008 /* seek. We already know that for LAMBDA = 0, Q is negative. */
1009 /* So we only need to find a value of LAMBDA that makes */
1010 /* Q positive. But... the expression for Q has singularities */
1011 /* at LAMBDA = -a**2, -b**2, and -c**2. */
1012
1013 /* These singularities are not necessarily to be avoided. */
1014 /* If the numerator of one of the terms for Q is zero, we */
1015 /* can simply compute Q ignoring that particular term. We */
1016 /* say that a singularity corresponding to a term whose */
1017 /* numerator is zero is a viable singularity. By being */
1018 /* careful in our computation of Q, we can assign LAMBDA to */
1019 /* the value of the singularity. A singularity that is not */
1020 /* viable is called a true singularity. */
1021
1022 /* By choosing LAMBDA just slightly greater than the largest */
1023 /* true singularity, we can bracket the root we seek. */
1024
1025 /* First we must decide which singularity is the first true */
1026 /* one. */
1027
1028 snglpt = 4;
1029 for (i__ = 3; i__ >= 1; --i__) {
1030 if (point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
1031 "point", i__1, "nearpt_", (ftnlen)1065)] != 0.) {
1032 snglpt = i__;
1033 }
1034 }
1035
1036 /* If there is a singular point, compute LAMBDA so that the */
1037 /* largest term of Q is equal to 4. */
1038
1039 if (snglpt <= 3) {
1040 for (i__ = 1; i__ <= 3; ++i__) {
1041 if (point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1042 s_rnge("point", i__1, "nearpt_", (ftnlen)1079)] ==
1043 0.) {
1044 tlambd[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1045 s_rnge("tlambd", i__1, "nearpt_", (ftnlen)
1046 1080)] = -axisqr[2];
1047 } else {
1048 tlambd[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1049 s_rnge("tlambd", i__1, "nearpt_", (ftnlen)
1050 1082)] = axis[(i__2 = i__ - 1) < 3 && 0 <=
1051 i__2 ? i__2 : s_rnge("axis", i__2, "nearpt_",
1052 (ftnlen)1082)] * ((d__1 = point[(i__3 = i__ -
1053 1) < 3 && 0 <= i__3 ? i__3 : s_rnge("point",
1054 i__3, "nearpt_", (ftnlen)1082)], abs(d__1)) *
1055 .5 - axis[(i__4 = i__ - 1) < 3 && 0 <= i__4 ?
1056 i__4 : s_rnge("axis", i__4, "nearpt_", (
1057 ftnlen)1082)]);
1058 }
1059 }
1060 /* Computing MAX */
1061 d__1 = max(tlambd[0],tlambd[1]);
1062 lambda = max(d__1,tlambd[2]);
1063 lower = lambda;
1064 upper = max(lower,0.);
1065 } else {
1066
1067 /* The point must be at the origin. In this case */
1068 /* we know where the closest point is. WE DON'T have */
1069 /* to compute anything. It's just at the end of the */
1070 /* shortest semi-major axis. However, since we */
1071 /* may have done some rounding off, we will make */
1072 /* sure that we pick the side of the shortest axis */
1073 /* that has the same sign as COPY(1). */
1074
1075 /* We are going to be a bit sneaky here. We know */
1076 /* where the closest point is so we are going to */
1077 /* simply make POINT and COPY equal to that point */
1078 /* and set the upper and lower bracketing bounds */
1079 /* to zero so that we won't have to deal with any */
1080 /* special cases later on. */
1081
1082 if (copy[0] < 0.) {
1083 point[0] = -axis[0];
1084 copy[0] = -axis[0];
1085 } else {
1086 point[0] = axis[0];
1087 copy[0] = axis[0];
1088 }
1089 copy[1] = 0.;
1090 copy[2] = 0.;
1091 upper = 0.;
1092 lower = 0.;
1093 lambda = 0.;
1094 q = 0.;
1095 inside = FALSE_;
1096 }
1097 }
1098
1099 /* OK. Now compute the value of Q at the two bracketing */
1100 /* values of LAMBDA. */
1101
1102 for (i__ = 1; i__ <= 3; ++i__) {
1103 if (point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("poi"
1104 "nt", i__1, "nearpt_", (ftnlen)1139)] == 0.) {
1105 term[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("term",
1106 i__1, "nearpt_", (ftnlen)1141)] = 0.;
1107 } else {
1108
1109 /* We have to be a bit careful for points inside the */
1110 /* ellipsoid. The denominator of the factor we are going to */
1111 /* compute is ( AXIS + LAMBDA/AXIS ). Numerically this may */
1112 /* be too close to zero for us to actually divide POINT by */
1113 /* it. However, since our solution algorithm for lambda */
1114 /* does not depend upon the differentiability of Q---in */
1115 /* fact it depends only on Q having the correct sign---we */
1116 /* can simply truncate its individual terms when we are in */
1117 /* danger of division overflows. */
1118 denom = axis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1119 s_rnge("axis", i__1, "nearpt_", (ftnlen)1155)] +
1120 lambda / axis[(i__2 = i__ - 1) < 3 && 0 <= i__2 ?
1121 i__2 : s_rnge("axis", i__2, "nearpt_", (ftnlen)1155)];
1122 trim = (d__1 = point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1
1123 : s_rnge("point", i__1, "nearpt_", (ftnlen)1157)],
1124 abs(d__1)) * .5 > denom;
1125 if (inside && trim) {
1126 factor = 2.;
1127 } else {
1128
1129 /* We don't expect DENOM to be zero here, but we'll */
1130 /* check anyway. */
1131
1132 if (denom == 0.) {
1133 setmsg_("AXIS(#) + LAMBDA/AXIS(#) is zero.", (ftnlen)
1134 33);
1135 errint_("#", &i__, (ftnlen)1);
1136 errint_("#", &i__, (ftnlen)1);
1137 sigerr_("SPICE(BUG)", (ftnlen)10);
1138 chkout_("NEARPT", (ftnlen)6);
1139 return 0;
1140 }
1141 factor = point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1142 s_rnge("point", i__1, "nearpt_", (ftnlen)1179)] /
1143 denom;
1144 }
1145 term[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("term",
1146 i__1, "nearpt_", (ftnlen)1183)] = factor * factor;
1147 }
1148 }
1149 if (! inside) {
1150 qlower = q;
1151 qupper = term[0] + term[1] + term[2] - 1.;
1152 } else {
1153 qupper = q;
1154 qlower = term[0] + term[1] + term[2] - 1.;
1155 }
1156
1157 /* Bracket QLOWER and QUPPER. */
1158
1159 qlower = max(0.,qlower);
1160 qupper = min(0.,qupper);
1161 lambda = upper;
1162 q = qupper;
1163
1164 /* Refining the estimate of lambda */
1165 /* ------------------------------- */
1166
1167 /* Now find the root of Q by bisection. */
1168
1169 itr = 0;
1170
1171 /* Throughout this loop we'll use TOUCHD to avoid logic problems */
1172 /* that may be caused by extended precision register usage by */
1173 /* some compilers. */
1174
1175 for(;;) { /* while(complicated condition) */
1176 d__1 = upper - lower;
1177 if (!(touchd_(&d__1) > 0.))
1178 break;
1179 ++itr;
1180 if (itr > 2048) {
1181 setmsg_("Iteration limit # exceeded in NEARPT. A, B, C = # #"
1182 " #; POSITN = ( #, #, # ). LOWER = #; UPPER = #; UPPE"
1183 "R-LOWER = #. Solution pass number = #. This event s"
1184 "hould never occur. Contact NAIF.", (ftnlen)187);
1185 errint_("#", &c__2048, (ftnlen)1);
1186 errdp_("#", a, (ftnlen)1);
1187 errdp_("#", b, (ftnlen)1);
1188 errdp_("#", c__, (ftnlen)1);
1189 errdp_("#", positn, (ftnlen)1);
1190 errdp_("#", &positn[1], (ftnlen)1);
1191 errdp_("#", &positn[2], (ftnlen)1);
1192 errdp_("#", &lower, (ftnlen)1);
1193 errdp_("#", &upper, (ftnlen)1);
1194 d__1 = upper - lower;
1195 errdp_("#", &d__1, (ftnlen)1);
1196 sigerr_("SPICE(BUG)", (ftnlen)10);
1197 chkout_("NEARPT", (ftnlen)6);
1198 return 0;
1199 }
1200
1201 /* Bracket LOWER, QLOWER, and QUPPER. */
1202
1203 lower = min(lower,upper);
1204 qlower = max(0.,qlower);
1205 qupper = min(0.,qupper);
1206
1207 /* Depending upon how Q compares with Q at the */
1208 /* bracketing endpoints we adjust the endpoints */
1209 /* of the bracketing interval */
1210
1211 if (q == 0.) {
1212
1213 /* We've found the root. */
1214
1215 lower = lambda;
1216 upper = lambda;
1217 } else {
1218 if (q < 0.) {
1219 upper = lambda;
1220 qupper = q;
1221 } else {
1222
1223 /* We have Q > 0 */
1224
1225 lower = lambda;
1226 qlower = q;
1227 }
1228
1229 /* Update LAMBDA. */
1230
1231 lambda = lower * .5 + upper * .5;
1232
1233 /* It's quite possible as we get close to the root for Q */
1234 /* that round off errors in the computation of the next */
1235 /* value of LAMBDA will push it outside of the current */
1236 /* bracketing interval. Force it back in to the current */
1237 /* interval. */
1238
1239 lambda = brcktd_(&lambda, &lower, &upper);
1240 }
1241
1242 /* At this point, it's guaranteed that */
1243
1244 /* LOWER < LAMBDA < UPPER */
1245 /* - - */
1246
1247 /* If we didn't find a root, we've set LAMBDA to the midpoint */
1248 /* of the previous values of LOWER and UPPER, and we've moved */
1249 /* either LOWER or UPPER to the old value of LAMBDA, thereby */
1250 /* halving the distance between LOWER and UPPER. */
1251
1252 /* If we are still at an endpoint, we might as well cash in */
1253 /* our chips. We aren't going to be able to get away from the */
1254 /* endpoints. Set LOWER and UPPER equal so that the loop will */
1255 /* finally terminate. */
1256
1257 if (approx_(&lambda, &lower, &c_b108) || approx_(&lambda, &upper,
1258 &c_b108)) {
1259
1260 /* Make the decision as to which way to push */
1261 /* the boundaries, by selecting that endpoint */
1262 /* at which Q is closest to zero. */
1263 if (abs(qlower) < abs(qupper)) {
1264 upper = lower;
1265 } else {
1266 lower = upper;
1267 }
1268
1269 /* Since LOWER is equal to UPPER, the loop will terminate. */
1270
1271 }
1272
1273 /* If LOWER and UPPER aren't the same, we compute the */
1274 /* value of Q at our new guess for LAMBDA. */
1275
1276 d__1 = upper - lower;
1277 if (touchd_(&d__1) > 0.) {
1278 for (i__ = 1; i__ <= 3; ++i__) {
1279 if (point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1280 s_rnge("point", i__1, "nearpt_", (ftnlen)1337)] ==
1281 0.) {
1282 term[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1283 s_rnge("term", i__1, "nearpt_", (ftnlen)1339)]
1284 = 0.;
1285 } else {
1286 denom = axis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1
1287 : s_rnge("axis", i__1, "nearpt_", (ftnlen)
1288 1343)] + lambda / axis[(i__2 = i__ - 1) < 3 &&
1289 0 <= i__2 ? i__2 : s_rnge("axis", i__2,
1290 "nearpt_", (ftnlen)1343)];
1291 trim = (d__1 = point[(i__1 = i__ - 1) < 3 && 0 <=
1292 i__1 ? i__1 : s_rnge("point", i__1, "nearpt_",
1293 (ftnlen)1345)], abs(d__1)) * .5 > denom;
1294 if (inside && trim) {
1295 factor = 2.;
1296 } else {
1297
1298 /* We don't expect DENOM to be zero here, but we'll */
1299 /* check anyway. */
1300
1301 if (denom == 0.) {
1302 setmsg_("AXIS(#) + LAMBDA/AXIS(#) is zero.", (
1303 ftnlen)33);
1304 errint_("#", &i__, (ftnlen)1);
1305 errint_("#", &i__, (ftnlen)1);
1306 sigerr_("SPICE(BUG)", (ftnlen)10);
1307 chkout_("NEARPT", (ftnlen)6);
1308 return 0;
1309 }
1310 factor = point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ?
1311 i__1 : s_rnge("point", i__1, "nearpt_", (
1312 ftnlen)1368)] / denom;
1313 }
1314 term[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1315 s_rnge("term", i__1, "nearpt_", (ftnlen)1372)]
1316 = factor * factor;
1317 }
1318 }
1319 d__1 = term[0] + term[1] + term[2] - 1.;
1320 q = touchd_(&d__1);
1321 }
1322
1323 /* Q(LAMBDA) has been set unless we've already found */
1324 /* a solution. */
1325
1326 /* Loop back through the bracketing refinement code. */
1327
1328 }
1329
1330 /* Now we have LAMBDA, compute the nearest point based upon */
1331 /* this value. */
1332
1333 lambda = lower;
1334 for (i__ = 1; i__ <= 3; ++i__) {
1335 if (point[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("poi"
1336 "nt", i__1, "nearpt_", (ftnlen)1398)] == 0.) {
1337 spoint[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
1338 "spoint", i__1, "nearpt_", (ftnlen)1400)] = 0.;
1339 } else {
1340 denom = lambda / axisqr[(i__1 = i__ - 1) < 3 && 0 <= i__1 ?
1341 i__1 : s_rnge("axisqr", i__1, "nearpt_", (ftnlen)1404)
1342 ] + 1.;
1343
1344 /* We don't expect that DENOM will be non-positive, but we */
1345 /* check for this case anyway. */
1346
1347 if (denom <= 0.) {
1348 setmsg_("Denominator in expression for SPOINT(#) is #.", (
1349 ftnlen)45);
1350 errint_("#", &i__, (ftnlen)1);
1351 errdp_("#", &denom, (ftnlen)1);
1352 sigerr_("SPICE(BUG)", (ftnlen)10);
1353 chkout_("NEARPT", (ftnlen)6);
1354 return 0;
1355 }
1356 spoint[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
1357 "spoint", i__1, "nearpt_", (ftnlen)1421)] = copy[(
1358 i__2 = i__ - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge(
1359 "copy", i__2, "nearpt_", (ftnlen)1421)] / denom;
1360 }
1361 }
1362
1363 /* Handling points on the central plane. */
1364 /* ------------------------------------- */
1365
1366 /* I suppose you thought you were done at this point. Not */
1367 /* necessarily. If POINT is INSIDE the ellipsoid and happens to */
1368 /* lie in the Y-Z plane, there is a possibility (perhaps even */
1369 /* likelihood) that the nearest point on the ellipsoid is NOT in */
1370 /* the Y-Z plane. we must consider this possibility next. */
1371
1372 if (inside && (snglpt == 2 || snglpt == 3)) {
1373
1374 /* There are two ways to get here. SNGLPT = 2 or SNGLPT = 3. */
1375 /* Fortunately these two cases can be handled simultaneously by */
1376 /* code. However, they are most easily understood if explained */
1377 /* separately. */
1378
1379 /* Case 1. SNGLPT = 2 */
1380
1381 /* The input to the lambda solution POINT lies in the Y-Z plane. */
1382 /* We have already detected one critical point of the */
1383 /* distance function to POINT restricted to the ellipsoid. */
1384 /* This point also lies in the Y-Z plane. However, when */
1385 /* POINT lies on the Y-Z plane close to the center of the */
1386 /* ellipsoid, there may be a point that is nearest that does */
1387 /* not lie in the Y-Z plane. Assuming the existence of such a */
1388 /* point, (x,y,z) it must satisfy the equations */
1389
1390 /* lambda*x */
1391 /* x + -------- = POINT(1) = 0 */
1392 /* a*a */
1393
1394
1395 /* lambda*y */
1396 /* y + -------- = POINT(2) */
1397 /* b*b */
1398
1399
1400 /* lambda*z */
1401 /* z + -------- = POINT(3) */
1402 /* c*c */
1403
1404
1405 /* Since we are assuming that this undetected solution (x,y,z) */
1406 /* does not have x equal to 0, it must be the case that */
1407
1408 /* lambda = -a*a. */
1409
1410 /* Because of this, we must have */
1411
1412 /* y = POINT(2) / ( 1 - (a**2/b**2) ) */
1413 /* z = POINT(3) / ( 1 - (a**2/c**2) ) */
1414
1415 /* The value of x is obtained by forcing */
1416
1417 /* (x/a)**2 + (y/b)**2 + (z/c)**2 = 1. */
1418
1419 /* This assumes of course that a and b are not equal. If a and */
1420 /* b are the same, then since POINT(2) is not zero, the */
1421 /* solution we found above by deflating the original ellipsoid */
1422 /* will find the near point. */
1423
1424 /* (If a and b are equal, the ellipsoid deflates to a */
1425 /* segment on the z-axis when lambda = -a**2. Since */
1426 /* POINT(2) is not zero, the deflating ellipsoid must pass */
1427 /* through POINT before it collapses to a segment.) */
1428
1429
1430 /* Case 2. SNGLPT = 3 */
1431
1432 /* The input to the lambda solution POINT lies on the Z-axis. */
1433 /* The solution obtained in the generic case above will */
1434 /* locate the critical point of the distance function */
1435 /* that lies on the Z. However, there will also be */
1436 /* critical points in the X-Z plane and Y-Z plane. The point */
1437 /* in the X-Z plane is the one to examine. Why? We are looking */
1438 /* for the point on the ellipsoid closest to POINT. It must */
1439 /* lie in one of these two planes. But the ellipse of */
1440 /* intersection with the X-Z plane fits inside the ellipse */
1441 /* of intersection with the Y-Z plane. Therefore the closest */
1442 /* point on the Y-Z ellipse must be at a greater distance than */
1443 /* the closest point on the X-Z ellipse. Thus, in solving */
1444 /* the equations */
1445
1446
1447 /* lambda*x */
1448 /* x + -------- = POINT(1) = 0 */
1449 /* a*a */
1450
1451
1452 /* lambda*y */
1453 /* y + -------- = POINT(2) = 0 */
1454 /* b*b */
1455
1456
1457 /* lambda*z */
1458 /* z + -------- = POINT(3) */
1459 /* c*c */
1460
1461
1462 /* We have */
1463
1464 /* lambda = -a*a, */
1465
1466 /* y = 0, */
1467
1468 /* z = POINT(3) / ( 1 - (a**2/c**2) ) */
1469
1470 /* The value of x is obtained by forcing */
1471
1472 /* (x/a)**2 + (y/b)**2 + (z/c)**2 = 1. */
1473
1474 /* This assumes that a and c are not equal. If */
1475 /* a and c are the same, then the solution we found above */
1476 /* by deflating the original ellipsoid, will find the */
1477 /* near point. */
1478
1479 /* ( If a = c, then the input ellipsoid is a sphere. */
1480 /* The ellipsoid will deflate to the center of the */
1481 /* sphere. Since our point is NOT at the center, */
1482 /* the deflating sphere will cross through */
1483 /* (x,y,z) before it collapses to a point ) */
1484
1485 /* We begin by assuming this extra point doesn't exist. */
1486
1487 extra = FALSE_;
1488
1489 /* Next let's note a few simple tests we can apply to */
1490 /* eliminate searching for an extra point. */
1491
1492 /* First of all the smallest axis must be different from */
1493 /* the axis associated with the first true singularity. */
1494
1495
1496 /* Next, whatever point we find, it must be true that */
1497
1498 /* |y| < b, |z| < c */
1499
1500 /* because of the condition on the absolute values, we must */
1501 /* have: */
1502
1503 /* | POINT(2) | <= b - a*(a/b) */
1504
1505 /* | POINT(3) | <= c - a*(a/c) */
1506
1507 if (axis[0] != axis[(i__1 = snglpt - 1) < 3 && 0 <= i__1 ? i__1 :
1508 s_rnge("axis", i__1, "nearpt_", (ftnlen)1575)] && abs(
1509 point[1]) <= axis[1] - axisqr[0] / axis[1] && abs(point[2]
1510 ) <= axis[2] - axisqr[0] / axis[2]) {
1511
1512 /* What happens next depends on whether the ellipsoid is */
1513 /* prolate or triaxial. */
1514
1515 if (axis[0] == axis[1]) {
1516
1517 /* This is the prolate case; we need to compute the */
1518 /* z component of the solution. */
1519
1520 denom3 = 1. - axisqr[0] / axisqr[2];
1521 if (denom3 > 0.) {
1522 epoint[1] = 0.;
1523
1524 /* Concerning the safety of the following division: */
1525
1526 /* - A divide-by-zero check has been done above. */
1527
1528 /* - The numerator can be squared without exceeding */
1529 /* DPMAX(). This was checked near the start of the */
1530 /* routine. */
1531
1532 /* - The denominator was computed as the difference */
1533 /* of 1.D0 and another number. The smallest */
1534 /* possible magnitude of a non-zero value of the */
1535 /* denominator is roughly 1.D-16, assuming IEEE */
1536 /* double precision numeric representation. */
1537
1538
1539 epoint[2] = point[2] / denom3;
1540
1541 /* See if these components can even be on the */
1542 /* ellipsoid... */
1543
1544 /* Computing 2nd power */
1545 d__1 = epoint[1] / axis[1];
1546 /* Computing 2nd power */
1547 d__2 = epoint[2] / axis[2];
1548 temp = 1. - d__1 * d__1 - d__2 * d__2;
1549 if (temp > 0.) {
1550
1551 /* ...and compute the x component of the point. */
1552
1553 epoint[0] = axis[0] * sqrt((max(0.,temp)));
1554 extra = TRUE_;
1555 }
1556 }
1557 } else {
1558
1559 /* This is the triaxial case. */
1560
1561 /* Compute the y and z components (2 and 3) of the extra */
1562 /* point. */
1563
1564 denom2 = 1. - axisqr[0] / axisqr[1];
1565 denom3 = 1. - axisqr[0] / axisqr[2];
1566
1567 /* We expect DENOM2 and DENOM3 will always be positive. */
1568 /* Nonetheless, we check to make sure this is the case. */
1569 /* If not, we don't compute the extra point. */
1570
1571 if (denom2 > 0. && denom3 > 0.) {
1572
1573 /* Concerning the safety of the following divisions: */
1574
1575 /* - Divide-by-zero checks have been done above. */
1576
1577 /* - The numerators can be squared without exceeding */
1578 /* DPMAX(). This was checked near the start of the */
1579 /* routine. */
1580
1581 /* - Each denominator was computed as the difference */
1582 /* of 1.D0 and another number. The smallest */
1583 /* possible magnitude of a non-zero value of */
1584 /* either denominator is roughly 1.D-16, assuming */
1585 /* IEEE double precision numeric representation. */
1586
1587 epoint[1] = point[1] / denom2;
1588 epoint[2] = point[2] / denom3;
1589
1590 /* See if these components can even be on the */
1591 /* ellipsoid... */
1592
1593 /* Computing 2nd power */
1594 d__1 = epoint[1] / axis[1];
1595 /* Computing 2nd power */
1596 d__2 = epoint[2] / axis[2];
1597 temp = 1. - d__1 * d__1 - d__2 * d__2;
1598 if (temp > 0.) {
1599
1600 /* ...and compute the x component of the point. */
1601
1602 epoint[0] = axis[0] * sqrt(temp);
1603 extra = TRUE_;
1604 }
1605 }
1606 }
1607 }
1608
1609 /* Ok. If an extra point is possible, check and see if it */
1610 /* is the one we are searching for. */
1611
1612 if (extra) {
1613 if (vdist_(epoint, point) < vdist_(spoint, point)) {
1614 vequ_(epoint, spoint);
1615 }
1616 }
1617 }
1618
1619 /* Decisions and initializations for sharpening the solution. */
1620 /* ================================================================ */
1621 /* ---------------------------------------------------------------- */
1622
1623 if (solutn == 1) {
1624
1625 /* The first solution for the nearest point may not be */
1626 /* very close to being on the ellipsoid. To */
1627 /* take care of this case, we next find the point on the */
1628 /* ellipsoid, closest to our first solution point. */
1629 /* (Ideally the normal line at this second point should */
1630 /* contain both the current solution point and the */
1631 /* original point). */
1632
1633 vequ_(spoint, point);
1634 vequ_(spoint, bestpt);
1635 bestht = vdist_(bestpt, orignl);
1636 } else if (solutn == 2) {
1637
1638 /* The current solution point will be very close to lying */
1639 /* on the ellipsoid. However, the normal at this solution */
1640 /* may not actually point toward the input point. */
1641
1642 /* With the current solution we can predict */
1643 /* the location of the input point. The difference between */
1644 /* this predicted point and the actual point can be used */
1645 /* to sharpen our estimate of the solution. */
1646
1647 /* The sharpening is performed by */
1648
1649 /* 1) Compute the vector ERR that gives the difference */
1650 /* between the input point (POSITN) and the point */
1651 /* computed using the solution point, normal and */
1652 /* altitude. */
1653
1654 /* 2) Find the component of ERR that is orthogonal to the */
1655 /* normal at the current solution point. If the point */
1656 /* is outside the ellipsoid, scale this component to */
1657 /* the approximate scale of the near point. We use */
1658 /* the scale factor */
1659
1660 /* ||near point|| / ||input point|| */
1661
1662 /* Call this scaled component ERRP. */
1663
1664 /* 3) Translate the solution point by ERRP to get POINT. */
1665
1666 /* 4) Find the point on the ellipsoid closest to POINT. */
1667 /* (step 4 is handled by the solution loop above). */
1668
1669
1670 /* First we need to compute the altitude */
1671
1672 height = sign * vdist_(spoint, orignl);
1673
1674 /* Next compute the difference between the input point and */
1675 /* the one we get by moving out along the normal at our */
1676 /* solution point by the computed altitude. */
1677
1678 surfnm_(axis, &axis[1], &axis[2], spoint, normal);
1679 for (i__ = 1; i__ <= 3; ++i__) {
1680 err[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("err",
1681 i__1, "nearpt_", (ftnlen)1774)] = orignl[(i__2 = i__
1682 - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("orignl", i__2,
1683 "nearpt_", (ftnlen)1774)] - spoint[(i__3 = i__ - 1) <
1684 3 && 0 <= i__3 ? i__3 : s_rnge("spoint", i__3, "near"
1685 "pt_", (ftnlen)1774)] - height * normal[(i__4 = i__ -
1686 1) < 3 && 0 <= i__4 ? i__4 : s_rnge("normal", i__4,
1687 "nearpt_", (ftnlen)1774)];
1688 }
1689
1690 /* Find the component of the error vector that is */
1691 /* perpendicular to the normal, and shift our solution */
1692 /* point by this component. */
1693
1694 vperp_(err, normal, errp);
1695
1696 /* The sign of the original point's altitude tells */
1697 /* us whether the point is outside the ellipsoid. */
1698
1699 if (sign >= 0.) {
1700
1701 /* Scale the transverse component down to the local radius */
1702 /* of the surface point. */
1703
1704 if (pnorm == 0.) {
1705
1706 /* Since the point is outside of the scaled ellipsoid, */
1707 /* we don't expect to arrive here. This is a backstop */
1708 /* check. */
1709
1710 setmsg_("Norm of scaled point is 0. POSITN = ( #, #, # )",
1711 (ftnlen)47);
1712 errdp_("#", positn, (ftnlen)1);
1713 errdp_("#", &positn[1], (ftnlen)1);
1714 errdp_("#", &positn[2], (ftnlen)1);
1715 sigerr_("SPICE(BUG)", (ftnlen)10);
1716 chkout_("NEARPT", (ftnlen)6);
1717 return 0;
1718 }
1719 d__1 = vnorm_(spoint) / pnorm;
1720 vsclip_(&d__1, errp);
1721 }
1722 vadd_(spoint, errp, point);
1723 olderr = vnorm_(err);
1724 bestht = height;
1725
1726 /* Finally store the current solution point, so that if */
1727 /* this sharpening doesn't improve our estimate of the */
1728 /* near point, we can just return our current best guess. */
1729
1730 vequ_(spoint, bestpt);
1731 } else if (solutn > 2) {
1732
1733 /* This branch exists for the purpose of testing our */
1734 /* "sharpened" solution and setting up for another sharpening */
1735 /* pass. */
1736
1737 /* We have already stored our best guess so far in BESTPT and */
1738 /* the vector ERR is the difference */
1739
1740 /* ORIGNL - ( BESTPT + BESTHT*NORMAL ) */
1741
1742 /* We have just computed a new candidate "best" near point. */
1743 /* SPOINT. */
1744
1745 /* If the error vector */
1746
1747 /* ORIGNL - ( SPOINT + HEIGHT*NORMAL) */
1748
1749 /* is shorter than our previous error, we will make SPOINT */
1750 /* our best guess and try to sharpen our estimate again. */
1751
1752 /* If our sharpening method hasn't improved things, we just */
1753 /* call it quits and go with our current best guess. */
1754
1755
1756 /* First compute the altitude, */
1757
1758 height = sign * vdist_(spoint, orignl);
1759
1760 /* ... compute the difference */
1761
1762 /* ORIGNL - SPOINT - HEIGHT*NORMAL, */
1763
1764 surfnm_(axis, &axis[1], &axis[2], spoint, normal);
1765 for (i__ = 1; i__ <= 3; ++i__) {
1766 err[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("err",
1767 i__1, "nearpt_", (ftnlen)1867)] = orignl[(i__2 = i__
1768 - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("orignl", i__2,
1769 "nearpt_", (ftnlen)1867)] - spoint[(i__3 = i__ - 1) <
1770 3 && 0 <= i__3 ? i__3 : s_rnge("spoint", i__3, "near"
1771 "pt_", (ftnlen)1867)] - height * normal[(i__4 = i__ -
1772 1) < 3 && 0 <= i__4 ? i__4 : s_rnge("normal", i__4,
1773 "nearpt_", (ftnlen)1867)];
1774 }
1775
1776 /* ...and determine the magnitude of the error due to our */
1777 /* sharpened estimate. */
1778
1779 newerr = vnorm_(err);
1780
1781 /* If the sharpened estimate yields a smaller error ... */
1782
1783 if (newerr < olderr) {
1784
1785 /* ...our current value of SPOINT becomes our new */
1786 /* best point and the current altitude becomes our */
1787 /* new altitude. */
1788
1789 olderr = newerr;
1790 bestht = height;
1791 vequ_(spoint, bestpt);
1792
1793 /* Next, if we haven't passed the limit on the number of */
1794 /* iterations allowed we prepare the initial point for our */
1795 /* "sharpening" estimate. */
1796
1797 if (solutn <= 6) {
1798 vperp_(err, normal, errp);
1799
1800 /* The sign of the original point's altitude tells */
1801 /* us whether the point is outside the ellipsoid. */
1802
1803 if (sign >= 0.) {
1804
1805 /* Scale the transverse component down to the local */
1806 /* radius of the surface point. */
1807
1808 if (pnorm == 0.) {
1809
1810 /* Since the point is outside of the scaled */
1811 /* ellipsoid, we don't expect to arrive here. */
1812 /* This is a backstop check. */
1813
1814 setmsg_("Norm of scaled point is 0. POSITN = ( #"
1815 ", #, # )", (ftnlen)47);
1816 errdp_("#", positn, (ftnlen)1);
1817 errdp_("#", &positn[1], (ftnlen)1);
1818 errdp_("#", &positn[2], (ftnlen)1);
1819 sigerr_("SPICE(BUG)", (ftnlen)10);
1820 chkout_("NEARPT", (ftnlen)6);
1821 return 0;
1822 }
1823 d__1 = vnorm_(spoint) / pnorm;
1824 vsclip_(&d__1, errp);
1825 }
1826 vadd_(spoint, errp, point);
1827 }
1828 } else {
1829
1830 /* If things didn't get better, there is no point in */
1831 /* going on. Just set the SOLVNG flag to .FALSE. to */
1832 /* terminate the outer loop. */
1833
1834 solvng = FALSE_;
1835 }
1836 }
1837
1838 /* Increment the solution counter so that eventually this */
1839 /* loop will terminate. */
1840
1841 ++solutn;
1842 solvng = solvng && solutn <= 6;
1843 }
1844
1845 /* Clean up. */
1846 /* ================================================================== */
1847 /* ------------------------------------------------------------------ */
1848
1849 /* Re-scale and re-order the components of our solution point. Scale */
1850 /* and copy the value of BESTHT into the output argument. */
1851
1852 d__1 = 1. / scale;
1853 vsclip_(&d__1, bestpt);
1854 for (i__ = 1; i__ <= 3; ++i__) {
1855 npoint[(i__2 = iorder[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
1856 s_rnge("iorder", i__1, "nearpt_", (ftnlen)1966)] - 1) < 3 &&
1857 0 <= i__2 ? i__2 : s_rnge("npoint", i__2, "nearpt_", (ftnlen)
1858 1966)] = bestpt[(i__3 = i__ - 1) < 3 && 0 <= i__3 ? i__3 :
1859 s_rnge("bestpt", i__3, "nearpt_", (ftnlen)1966)];
1860 }
1861 *alt = bestht / scale;
1862 chkout_("NEARPT", (ftnlen)6);
1863 return 0;
1864 } /* nearpt_ */
1865
1866