1 /* npelpt.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__1 = 1;
11 static integer c__2 = 2;
12 static doublereal c_b10 = 0.;
13 static doublereal c_b11 = 1.;
14 static doublereal c_b12 = 2.;
15
16 /* $Procedure NPELPT ( Nearest point on ellipse to point ) */
npelpt_(doublereal * point,doublereal * ellips,doublereal * pnear,doublereal * dist)17 /* Subroutine */ int npelpt_(doublereal *point, doublereal *ellips,
18 doublereal *pnear, doublereal *dist)
19 {
20 /* System generated locals */
21 doublereal d__1;
22
23 /* Local variables */
24 extern /* Subroutine */ int vadd_(doublereal *, doublereal *, doublereal *
25 ), vsub_(doublereal *, doublereal *, doublereal *), vequ_(
26 doublereal *, doublereal *), mtxv_(doublereal *, doublereal *,
27 doublereal *);
28 doublereal scale;
29 extern /* Subroutine */ int chkin_(char *, ftnlen), vpack_(doublereal *,
30 doublereal *, doublereal *, doublereal *), errdp_(char *,
31 doublereal *, ftnlen);
32 extern doublereal vdist_(doublereal *, doublereal *);
33 doublereal tempv[3];
34 extern doublereal vnorm_(doublereal *);
35 extern /* Subroutine */ int el2cgv_(doublereal *, doublereal *,
36 doublereal *, doublereal *);
37 doublereal majlen, center[3], minlen;
38 extern /* Subroutine */ int nearpt_(doublereal *, doublereal *,
39 doublereal *, doublereal *, doublereal *, doublereal *);
40 doublereal smajor[3];
41 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
42 ftnlen);
43 doublereal rotate[9] /* was [3][3] */;
44 extern /* Subroutine */ int vsclip_(doublereal *, doublereal *), setmsg_(
45 char *, ftnlen);
46 doublereal sminor[3];
47 extern /* Subroutine */ int twovec_(doublereal *, integer *, doublereal *,
48 integer *, doublereal *);
49 doublereal prjpnt[3];
50 extern logical return_(void);
51 doublereal tmppnt[3];
52 extern /* Subroutine */ int mxv_(doublereal *, doublereal *, doublereal *)
53 ;
54
55 /* $ Abstract */
56
57 /* Find the nearest point on an ellipse to a specified point, both */
58 /* in three-dimensional space, and find the distance between the */
59 /* ellipse and the point. */
60
61 /* $ Disclaimer */
62
63 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
64 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
65 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
66 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
67 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
68 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
69 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
70 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
71 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
72 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
73
74 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
75 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
76 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
77 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
78 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
79 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
80
81 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
82 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
83 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
84 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
85
86 /* $ Required_Reading */
87
88 /* ELLIPSES */
89
90 /* $ Keywords */
91
92 /* CONIC */
93 /* ELLIPSE */
94 /* GEOMETRY */
95 /* MATH */
96
97 /* $ Declarations */
98 /* $ Brief_I/O */
99
100 /* Variable I/O Description */
101 /* -------- --- -------------------------------------------------- */
102 /* POINT I Point whose distance to an ellipse is to be found. */
103 /* ELLIPS I A SPICELIB ellipse. */
104 /* PNEAR O Nearest point on ellipse to input point. */
105 /* DIST O Distance of input point to ellipse. */
106
107 /* $ Detailed_Input */
108
109 /* ELLIPS is a SPICELIB ellipse that represents an ellipse */
110 /* in three-dimensional space. */
111
112 /* POINT is a point in 3-dimensional space. */
113
114 /* $ Detailed_Output */
115
116 /* PNEAR is the nearest point on ELLIPS to POINT. */
117
118 /* DIST is the distance between POINT and PNEAR. This is */
119 /* the distance between POINT and the ellipse. */
120
121 /* $ Parameters */
122
123 /* None. */
124
125 /* $ Exceptions */
126
127 /* 1) Invalid ellipses will be diagnosed by routines called by */
128 /* this routine. */
129
130 /* 2) Ellipses having one or both semi-axis lengths equal to zero */
131 /* are turned away at the door; the error SPICE(DEGENERATECASE) */
132 /* is signalled. */
133
134 /* 3) If the geometric ellipse represented by ELLIPS does not */
135 /* have a unique point nearest to the input point, any point */
136 /* at which the minimum distance is attained may be returned */
137 /* in PNEAR. */
138
139 /* $ Files */
140
141 /* None. */
142
143 /* $ Particulars */
144
145 /* Given an ellipse and a point in 3-dimensional space, if the */
146 /* orthogonal projection of the point onto the plane of the ellipse */
147 /* is on or outside of the ellipse, then there is a unique point on */
148 /* the ellipse closest to the original point. This routine finds */
149 /* that nearest point on the ellipse. If the projection falls inside */
150 /* the ellipse, there may be multiple points on the ellipse that are */
151 /* at the minimum distance from the original point. In this case, */
152 /* one such closest point will be returned. */
153
154 /* This routine returns a distance, rather than an altitude, in */
155 /* contrast to the SPICELIB routine NEARPT. Because our ellipse is */
156 /* situated in 3-space and not 2-space, the input point is not */
157 /* `inside' or `outside' the ellipse, so the notion of altitude does */
158 /* not apply to the problem solved by this routine. In the case of */
159 /* NEARPT, the input point is on, inside, or outside the ellipsoid, */
160 /* so it makes sense to speak of its altitude. */
161
162 /* $ Examples */
163
164 /* 1) For planetary rings that can be modelled as flat disks with */
165 /* elliptical outer boundaries, the distance of a point in */
166 /* space from a ring's outer boundary can be computed using this */
167 /* routine. Suppose CENTER, SMAJOR, and SMINOR are the center, */
168 /* semi-major axis, and semi-minor axis of the ring's boundary. */
169 /* Suppose also that SCPOS is the position of a spacecraft. */
170 /* SCPOS, CENTER, SMAJOR, and SMINOR must all be expressed in */
171 /* the same coordinate system. We can find the distance from */
172 /* the spacecraft to the ring using the code fragment */
173
174 /* C */
175 /* C Make a SPICELIB ellipse representing the ring, */
176 /* C then use NPELPT to find the distance between */
177 /* C the spacecraft position and RING. */
178 /* C */
179 /* CALL CGV2EL ( CENTER, SMAJOR, SMINOR, RING ) */
180 /* CALL NPELPT ( SCPOS, RING, PNEAR, DIST ) */
181
182
183 /* 2) The problem of finding the distance of a line from a tri-axial */
184 /* ellipsoid can be reduced to the problem of finding the */
185 /* distance between the same line and an ellipse; this problem in */
186 /* turn can be reduced to the problem of finding the distance */
187 /* between an ellipse and a point. The routine NPEDLN carries */
188 /* out this process and uses NPELPT to find the ellipse-to-point */
189 /* distance. */
190
191
192 /* $ Restrictions */
193
194 /* None. */
195
196 /* $ Literature_References */
197
198 /* None. */
199
200 /* $ Author_and_Institution */
201
202 /* N.J. Bachman (JPL) */
203
204 /* $ Version */
205
206 /* - SPICELIB Version 1.2.0, 02-SEP-2005 (NJB) */
207
208 /* Updated to remove non-standard use of duplicate arguments */
209 /* in VADD, VSCL, MTXV and MXV calls. */
210
211 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
212
213 /* Comment section for permuted index source lines was added */
214 /* following the header. */
215
216 /* - SPICELIB Version 1.0.0, 02-NOV-1990 (NJB) */
217
218 /* -& */
219 /* $ Index_Entries */
220
221 /* nearest point on ellipse to point */
222
223 /* -& */
224 /* $ Revisions */
225
226 /* - SPICELIB Version 1.2.0, 02-SEP-2005 (NJB) */
227
228 /* Updated to remove non-standard use of duplicate arguments */
229 /* in VADD, VSCL, MTXV and MXV calls. */
230
231 /* -& */
232
233 /* SPICELIB functions */
234
235
236 /* Local variables */
237
238
239 /* Standard SPICE error handling. */
240
241 if (return_()) {
242 return 0;
243 } else {
244 chkin_("NPELPT", (ftnlen)6);
245 }
246
247 /* Here's an overview of our solution: */
248
249 /* Let ELPL be the plane containing the ELLIPS, and let PRJ be */
250 /* the orthogonal projection of the POINT onto ELPL. Let X be */
251 /* any point in the plane ELPL. According to the Pythagorean */
252 /* Theorem, */
253
254 /* 2 2 2 */
255 /* || POINT - X || = || POINT - PRJ || + || PRJ - X ||. */
256
257 /* Then if we can find a point X on ELLIPS that minimizes the */
258 /* rightmost term, that point X is the closest point on the */
259 /* ellipse to POINT. */
260
261 /* So, we find the projection PRJ, and then solve the problem of */
262 /* finding the closest point on ELLIPS to PRJ. To solve this */
263 /* problem, we find a triaxial ellipsoid whose intersection with */
264 /* the plane ELPL is precisely ELLIPS, and two of whose axes lie */
265 /* in the plane ELPL. The closest point on ELLIPS to PRJ is also */
266 /* the closest point on the ellipsoid to ELLIPS. But we have the */
267 /* SPICELIB routine NEARPT on hand to find the closest point on an */
268 /* ellipsoid to a specified point, so we've reduced our problem to */
269 /* a solved problem. */
270
271 /* There is a subtle point to worry about here: if PRJ is outside */
272 /* of ELLIPS (PRJ is in the same plane as ELLIPS, so `outside' */
273 /* does make sense here), then the closest point on ELLIPS to PRJ */
274 /* coincides with the closest point on the ellipsoid to PRJ, */
275 /* regardless of how we choose the z-semi-axis length of the */
276 /* ellipsoid. But the correspondence may be lost if PRJ is inside */
277 /* the ellipse, if we don't choose the z-semi-axis length */
278 /* correctly. */
279
280 /* Though it takes some thought to verify this (and we won't prove */
281 /* it here), making the z-semi-axis of the ellipsoid longer than */
282 /* the other two semi-axes is sufficient to maintain the */
283 /* coincidence of the closest point on the ellipsoid to PRJPNT and */
284 /* the closest point on the ellipse to PRJPNT. */
285
286
287 /* Find the ellipse's center and semi-axes. */
288
289 el2cgv_(ellips, center, smajor, sminor);
290
291 /* Find the lengths of the semi-axes, and scale the vectors to try */
292 /* to prevent arithmetic unpleasantness. Degenerate ellipses are */
293 /* turned away at the door. */
294
295 minlen = vnorm_(sminor);
296 majlen = vnorm_(smajor);
297 if (min(majlen,minlen) == 0.) {
298 setmsg_("Semi-axis lengths: # #. ", (ftnlen)24);
299 errdp_("#", &majlen, (ftnlen)1);
300 errdp_("#", &minlen, (ftnlen)1);
301 sigerr_("SPICE(DEGENERATECASE)", (ftnlen)21);
302 chkout_("NPELPT", (ftnlen)6);
303 return 0;
304 }
305 scale = 1. / majlen;
306 vsclip_(&scale, smajor);
307 vsclip_(&scale, sminor);
308
309 /* Translate ellipse and point so that the ellipse is centered at */
310 /* the origin. Scale the point's coordinates to maintain the */
311 /* correct relative position to the scaled ellipse. */
312
313 vsub_(point, center, tmppnt);
314 vsclip_(&scale, tmppnt);
315
316 /* We want to reduce the problem to a two-dimensional one. We'll */
317 /* work in a coordinate system whose x- and y- axes are aligned with */
318 /* the semi-major and semi-minor axes of the input ellipse. The */
319 /* z-axis is picked to give us a right-handed system. We find the */
320 /* matrix that transforms coordinates to our new system using TWOVEC. */
321
322 twovec_(smajor, &c__1, sminor, &c__2, rotate);
323
324 /* Apply the coordinate transformation to our scaled input point. */
325
326 mxv_(rotate, tmppnt, tempv);
327 vequ_(tempv, tmppnt);
328
329 /* We must find the distance between the orthogonal projection of */
330 /* TMPPNT onto the x-y plane and the ellipse. The projection is */
331 /* just */
332
333 /* ( TMPPNT(1), TMPPNT(2), 0 ); */
334
335 /* we'll call this projection PRJPNT. */
336
337
338 vpack_(tmppnt, &tmppnt[1], &c_b10, prjpnt);
339
340 /* Now we're ready to find the distance between and a triaxial */
341 /* ellipsoid whose intersection with the x-y plane is the ellipse */
342 /* and whose third semi-axis lies on the z-axis. */
343
344 /* Because we've scaled the ellipse's axes so as to give the longer */
345 /* axis length 1, a length of 2.D0 suffices for the ellipsoid's */
346 /* z-semi-axis. */
347
348
349 /* Find the nearest point to PRJPNT on the ellipoid, PNEAR. */
350
351 d__1 = minlen / majlen;
352 nearpt_(prjpnt, &c_b11, &d__1, &c_b12, pnear, dist);
353
354 /* Scale the near point coordinates back to the original scale. */
355
356 vsclip_(&majlen, pnear);
357
358 /* Apply the required inverse rotation and translation to PNEAR. */
359 /* Compute the point-to-ellipse distance. */
360
361 mtxv_(rotate, pnear, tempv);
362 vadd_(tempv, center, pnear);
363 *dist = vdist_(pnear, point);
364 chkout_("NPELPT", (ftnlen)6);
365 return 0;
366 } /* npelpt_ */
367
368