1 /* zzsglatx.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_b5 = 0.;
11
12 /* $Procedure ZZSGLATX ( Line segment latitude extent ) */
zzsglatx_(doublereal * p1,doublereal * p2,doublereal * minlat,doublereal * minp,doublereal * maxlat,doublereal * maxp)13 /* Subroutine */ int zzsglatx_(doublereal *p1, doublereal *p2, doublereal *
14 minlat, doublereal *minp, doublereal *maxlat, doublereal *maxp)
15 {
16 /* Initialized data */
17
18 static doublereal z__[3] = { 0.,0.,1. };
19
20 extern doublereal vdot_(doublereal *, doublereal *);
21 extern /* Subroutine */ int vsub_(doublereal *, doublereal *, doublereal *
22 ), vequ_(doublereal *, doublereal *);
23 doublereal r__, t[3];
24 extern /* Subroutine */ int chkin_(char *, ftnlen), vcrss_(doublereal *,
25 doublereal *, doublereal *);
26 extern logical vzero_(doublereal *);
27 integer nxpts;
28 doublereal plane2[4];
29 extern /* Subroutine */ int nvc2pl_(doublereal *, doublereal *,
30 doublereal *);
31 extern logical failed_(void);
32 doublereal crease[3];
33 extern /* Subroutine */ int reclat_(doublereal *, doublereal *,
34 doublereal *, doublereal *);
35 doublereal normal[3];
36 extern logical opsgnd_(doublereal *, doublereal *);
37 extern /* Subroutine */ int vhatip_(doublereal *), chkout_(char *, ftnlen)
38 ;
39 doublereal dp1, dp2;
40 extern /* Subroutine */ int inrypl_(doublereal *, doublereal *,
41 doublereal *, integer *, doublereal *);
42 extern logical return_(void);
43 doublereal dir[3], lat, lon, lat1, lat2;
44
45 /* $ Abstract */
46
47 /* SPICE Private routine intended solely for the support of SPICE */
48 /* routines. Users should not call this routine directly due */
49 /* to the volatile nature of this routine. */
50
51 /* Find the latitude extent (extrema of latitude) of a line segment. */
52 /* Latitude is defined in the planetocentric sense. */
53
54 /* $ Disclaimer */
55
56 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
57 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
58 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
59 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
60 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
61 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
62 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
63 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
64 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
65 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
66
67 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
68 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
69 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
70 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
71 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
72 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
73
74 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
75 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
76 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
77 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
78
79 /* $ Required_Reading */
80
81 /* None. */
82
83 /* $ Keywords */
84
85 /* EXTREMA */
86 /* GEOMETRY */
87 /* LATITUDE */
88 /* LINE */
89
90 /* $ Declarations */
91 /* $ Brief_I/O */
92
93 /* Variable I/O Description */
94 /* -------- --- -------------------------------------------------- */
95 /* P1 I First line segment endpoint. */
96 /* P2 I Second line segment endpoint. */
97 /* MINLAT O Minimum latitude on segment. */
98 /* MINP O Point where minimum latitude is attained. */
99 /* MAXLAT O Maximum latitude on segment. */
100 /* MAXP O Point where maximum latitude is attained. */
101 /* UBPL P SPICE plane upper bound. */
102
103 /* $ Detailed_Input */
104
105 /* P1, */
106 /* P2 are endpoints of a given line segment. */
107
108 /* $ Detailed_Output */
109
110
111 /* MINLAT is the minimum planetocentric (latitudinal) latitude */
112 /* that occurs on the line segment defined by the input */
113 /* endpoints. Units are radians. */
114
115 /* MINP is a point on the line segment at which the latitude */
116 /* MINLAT is attained. Note that in some cases, the */
117 /* minimum latitude can occur at more than one point. */
118
119 /* MAXLAT is the maximum planetocentric (latitudinal) latitude */
120 /* that occurs on the line segment defined by the input */
121 /* endpoints. Units are radians. */
122
123 /* MAXP is a point on the line segment at which the latitude */
124 /* MAXLAT is attained. Note that in some cases, the */
125 /* maximum latitude can occur at more than one point. */
126
127 /* $ Parameters */
128
129 /* None. */
130
131 /* $ Exceptions */
132
133 /* 1) If an error occurs when this routine attempts to construct a */
134 /* SPICE plane, the error will be diagnosed by routines in the */
135 /* call tree of this routine. */
136
137 /* 2) If an error occurs when this routine attempts to compute a */
138 /* ray-plane intersection, the error will be diagnosed by */
139 /* routines in the call tree of this routine. */
140
141 /* $ Files */
142
143 /* None. */
144
145 /* $ Particulars */
146
147
148 /* In this routine, "latitude" of a point refers to planetocentric */
149 /* latitude: the angle between the X-Y plane and the ray emanating */
150 /* from the origin and passing through the point. */
151
152 /* The term "latitude extents" on a line segment refers to the */
153 /* minimum and maximum values of latitude attained on that segment. */
154 /* The subset of the line segment at which a latitude extremum is */
155 /* attained can be an endpoint, a single interior point of the */
156 /* segment, or the whole segment. */
157
158 /* There are several geometric cases that can occur: */
159
160 /* Endpoints of segment are linearly dependent */
161 /* =========================================== */
162
163 /* The latitude extrema occur at the segment endpoints. Note that */
164 /* in this case the entire segment is tangent to one or both */
165 /* nappes of a vertical cone. The latitude values attained on the */
166 /* segment can be any non-empty subset of */
167
168 /* { minimum latitude, 0, maximum latitude } */
169
170
171 /* Segment lies in X-Y plane */
172 /* ========================= */
173
174 /* The latitude extrema are both zero. */
175
176
177 /* Segment intersects Z-axis in a single non-origin point */
178 /* ====================================================== */
179
180 /* One latitude extremum will occur on the Z-axis; the other */
181 /* extremum will occur at a segment endpoint. This case may be */
182 /* viewed as a degenerate version of the "local extremum exists */
183 /* in segment" case. */
184
185
186 /* Local extremum exists in segment */
187 /* ================================ */
188
189 /* If an extremum occurs at a single point T in the interior of */
190 /* the line segment, we call this a local extremum. Presuming the */
191 /* segment does not intersect the Z-axis, then T is a tangent */
192 /* point on a vertical cone C0 having its apex at the origin. All */
193 /* points, excluding the origin, on the nappe of the cone */
194 /* containing T have latitude equal to that of T. */
195
196 /* Let P0 be the plane that contains the line segment and the */
197 /* origin. If the endpoints of the line segment are linearly */
198 /* independent, P0 is uniquely defined, and P0 is tangent to the */
199 /* cone C0 at T. */
200
201 /* Let N0 be a normal vector of P0. Let P1 be a plane containing */
202 /* N0 and the Z-axis. If N0 is not parallel to the Z-axis, P1 is */
203 /* uniquely defined, and the point T lies in the plane P1 */
204 /* containing the Z-axis and N0: the intersection of the segment */
205 /* with P1 is T. */
206
207 /* Three of the cases excluded here */
208
209 /* - Segment endpoints are linearly dependent */
210
211 /* - Segment intersects the Z-axis, not at the origin */
212
213 /* - Normal to P0 is parallel to Z-axis (segment lies in */
214 /* X-Y plane) */
215
216 /* are all discussed above. The remaining case is discussed */
217 /* below. */
218
219
220 /* Local extremum occurs in extension of segment */
221 /* ============================================= */
222
223 /* A local extremum may exist on the line containing the segment, */
224 /* outside of the segment. In this case the segment cannot */
225 /* contain a point where a local extremum occurs. The extrema */
226 /* of latitude on the segment occur at the endpoints. */
227
228
229 /* $ Examples */
230
231 /* 1) One application of this routine is to assist in subsetting */
232 /* a type 2 DSK segment; this task requires determination of */
233 /* a set of triangular plates that intersect a given */
234 /* longitude/latitude rectangle. Note that this set cannot be */
235 /* found by examination of vertex latitudes alone; latitude */
236 /* extents of edges are needed. */
237
238 /* 2) This routine is also used to find the maximum latitude */
239 /* change that can occur between two points on a given ray. */
240
241 /* $ Restrictions */
242
243 /* None. */
244
245 /* $ Literature_References */
246
247 /* None. */
248
249 /* $ Author_and_Institution */
250
251 /* N.J. Bachman (JPL) */
252
253 /* $ Version */
254
255 /* - SPICELIB Version 1.0.0, 29-SEP-2016 (NJB) */
256
257 /* Original version 06-DEC-2012 (NJB) */
258
259 /* -& */
260 /* $ Index_Entries */
261
262 /* find planetocentric latitude extent of a line segment */
263 /* find planetocentric latitude range of a line segment */
264 /* find planetocentric latitude extrema of a line segment */
265
266 /* -& */
267
268 /* SPICELIB functions */
269
270
271 /* Local variables */
272
273
274 /* Saved variables */
275
276
277 /* Initial values */
278
279
280 /* Standard SPICE error handling. */
281
282 if (return_()) {
283 return 0;
284 }
285 chkin_("ZZSGLATX", (ftnlen)8);
286
287 /* Start by computing latitude at the segment's endpoints. */
288
289 reclat_(p1, &r__, &lon, &lat1);
290 reclat_(p2, &r__, &lon, &lat2);
291
292 /* Initialize the outputs using latitudes of the endpoints. */
293 /* If there are interior extrema, we'll update these outputs */
294 /* as needed. */
295
296 if (lat1 <= lat2) {
297 *minlat = lat1;
298 *maxlat = lat2;
299 vequ_(p1, minp);
300 vequ_(p2, maxp);
301 } else {
302 *minlat = lat2;
303 *maxlat = lat1;
304 vequ_(p2, minp);
305 vequ_(p1, maxp);
306 }
307
308 /* We want to work with the plane containing the origin, P1, and P2. */
309 /* We'll call this plane PLANE1. First see whether P1 and P2 are */
310 /* linearly independent. */
311
312 vcrss_(p1, p2, normal);
313 if (vzero_(normal)) {
314
315 /* We have a special case: P1 and P2 define a line passing */
316 /* through the origin. The latitude extrema lie on the */
317 /* segment endpoints, and possibly at every point on the */
318 /* segment. We've already computed the outputs. */
319
320 chkout_("ZZSGLATX", (ftnlen)8);
321 return 0;
322 }
323
324 /* At this point we know that NORMAL is non-zero. Convert it */
325 /* to a unit vector. */
326
327 vhatip_(normal);
328
329 /* Let ALPHA be the non-negative angle between PLANE1 and the X-Y */
330 /* plane. Then ALPHA and -ALPHA are, respectively, the maximum and */
331 /* minimum possible latitudes attained on the input segment. */
332 /* However, these values are not necessarily attained on the */
333 /* segment; we'll need to perform further analysis to find out. We */
334 /* don't need to compute ALPHA, but we'll refer to it in the */
335 /* discussion below. */
336
337 /* The next step is to find the normal vector to the plane defined */
338 /* by Z and NORMAL. We'll call this plane PLANE2. This plane might */
339 /* not exist if NORMAL and Z are linearly dependent. If PLANE2 */
340 /* does exist, the X-Y plane and PLANE1 intersect in a "crease" */
341 /* that is normal to PLANE2. */
342
343 vcrss_(z__, normal, crease);
344 if (vzero_(crease)) {
345
346 /* Z and NORMAL are linearly dependent; PLANE1 coincides (up to */
347 /* round-off error) with the X-Y plane. We've already computed */
348 /* the outputs. */
349
350 chkout_("ZZSGLATX", (ftnlen)8);
351 return 0;
352 }
353
354 /* At this point we know CREASE is non-zero. Convert */
355 /* it to a unit vector. */
356
357 vhatip_(crease);
358
359 /* By construction, CREASE is orthogonal to NORMAL. PLANE2 */
360 /* cuts PLANE1 in a line L passing through the origin. If */
361 /* the line segment has an interior latitude extremum, */
362 /* the point T where that extremum is attained lies on L. */
363 /* The segment is tangent at T to a nappe of a cone, centered on */
364 /* the Z-axis and having its apex at the origin, for which */
365 /* the half-angle is (pi/2)-ALPHA. The point T lies in PLANE2 */
366 /* since L is contained in PLANE2. */
367
368 /* If a single tangent point T exists in the interior of the */
369 /* segment, then the endpoints must be on opposite sides of PLANE2. */
370 /* See whether this is the case. */
371
372 dp1 = vdot_(p1, crease);
373 dp2 = vdot_(p2, crease);
374 if (opsgnd_(&dp1, &dp2)) {
375
376 /* The segment crosses PLANE2 at an interior point; this */
377 /* point is where the extremum occurs. Solve for the */
378 /* intersection. */
379
380 /* CREASE is guaranteed to be a unit vector. A zero input */
381 /* vector is the only cause for which NVC2PL will signal */
382 /* an error. Therefore we don't check FAILED after the */
383 /* following call. */
384
385 nvc2pl_(crease, &c_b5, plane2);
386 vsub_(p2, p1, dir);
387 inrypl_(p1, dir, plane2, &nxpts, t);
388 if (failed_()) {
389 chkout_("ZZSGLATX", (ftnlen)8);
390 return 0;
391 }
392 if (nxpts == 1) {
393
394 /* This is the normal case: we have one intersection of the */
395 /* segment with PLANE2. Update whichever of the extrema is */
396 /* superseded by the interior value. */
397
398 /* Note that this case can occur when NORMAL is orthogonal to */
399 /* Z, making ALPHA equal to pi/2. The nappes are degenerate in */
400 /* this case, consisting of the positive and negative Z-axes. */
401 /* This degenerate case occurs when the segment intersects the */
402 /* Z-axis in a point other than the origin, and the endpoints */
403 /* are linearly independent. */
404
405 /* This is not a special case computationally. */
406
407 reclat_(t, &r__, &lon, &lat);
408 if (lat > *maxlat) {
409 *maxlat = lat;
410 vequ_(t, maxp);
411 } else if (lat < *minlat) {
412 *minlat = lat;
413 vequ_(t, minp);
414 }
415
416 /* There can be only one local extremum, so we're done. */
417
418 }
419
420 /* If NXPTS is not 1, then even though the endpoints are on */
421 /* opposite sides of PLANE2, either the segment was found to lie */
422 /* in PLANE2 or no intersection was found. This situation must be */
423 /* due to finite precision arithmetic error. We'll make do with */
424 /* the extrema already found. */
425 }
426
427 /* We reach this point if we found a local extremum or if any of the */
428 /* following are true: */
429
430 /* 1) The segment misses PLANE2 altogether, in which case */
431 /* there's no tangency point. */
432
433 /* 2) One endpoint lies on PLANE2 and one endpoint does not. */
434
435 /* 3) Both endpoints lie in PLANE2. Then both endpoints lie */
436 /* in L, so we should have found them to be linearly */
437 /* dependent. This situation must be due to finite precision */
438 /* arithmetic error. */
439
440 /* In all of the numbered cases the extrema occur at the endpoints. */
441 /* and have been found already. In all cases, the outputs are set. */
442
443 chkout_("ZZSGLATX", (ftnlen)8);
444 return 0;
445 } /* zzsglatx_ */
446
447