1 /* zzgflong.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__15 = 15;
11 static integer c__7 = 7;
12 static integer c__0 = 0;
13 static integer c__1 = 1;
14 static doublereal c_b69 = 1.;
15 static doublereal c_b70 = 0.;
16
17 /* $Procedure ZZGFLONG ( GF, longitude solver ) */
zzgflong_(char * vecdef,char * method,char * target,char * ref,char * abcorr,char * obsrvr,char * dref,doublereal * dvec,char * crdsys,char * crdnam,char * relate,doublereal * refval,doublereal * tol,doublereal * adjust,U_fp udstep,U_fp udrefn,logical * rpt,U_fp udrepi,U_fp udrepu,U_fp udrepf,logical * bail,L_fp udbail,integer * mw,integer * nw,doublereal * work,doublereal * cnfine,doublereal * result,ftnlen vecdef_len,ftnlen method_len,ftnlen target_len,ftnlen ref_len,ftnlen abcorr_len,ftnlen obsrvr_len,ftnlen dref_len,ftnlen crdsys_len,ftnlen crdnam_len,ftnlen relate_len)18 /* Subroutine */ int zzgflong_(char *vecdef, char *method, char *target, char
19 *ref, char *abcorr, char *obsrvr, char *dref, doublereal *dvec, char *
20 crdsys, char *crdnam, char *relate, doublereal *refval, doublereal *
21 tol, doublereal *adjust, U_fp udstep, U_fp udrefn, logical *rpt, U_fp
22 udrepi, U_fp udrepu, U_fp udrepf, logical *bail, L_fp udbail, integer
23 *mw, integer *nw, doublereal *work, doublereal *cnfine, doublereal *
24 result, ftnlen vecdef_len, ftnlen method_len, ftnlen target_len,
25 ftnlen ref_len, ftnlen abcorr_len, ftnlen obsrvr_len, ftnlen dref_len,
26 ftnlen crdsys_len, ftnlen crdnam_len, ftnlen relate_len)
27 {
28 /* Initialized data */
29
30 static char ops[6*7] = "< " "= " "> " "LOCMIN" "ABSMIN" "LOC"
31 "MAX" "ABSMAX";
32 static doublereal y[3] = { 0.,1.,0. };
33
34 /* System generated locals */
35 integer work_dim1, work_dim2, work_offset, i__1, i__2, i__3, i__4;
36 doublereal d__1, d__2;
37
38 /* Builtin functions */
39 integer s_rnge(char *, integer, char *, integer), s_cmp(char *, char *,
40 ftnlen, ftnlen);
41 double cos(doublereal);
42 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
43 double sqrt(doublereal), sin(doublereal), atan2(doublereal, doublereal);
44
45 /* Local variables */
46 integer head, node, left, quad;
47 logical flip;
48 integer next;
49 extern /* Subroutine */ int zzgfcodc_(), zzgfcocd_();
50 extern /* Subroutine */ int zzgfcocg_(doublereal *, doublereal *);
51 extern /* Subroutine */ int zzgfcosd_();
52 extern /* Subroutine */ int zzgfcoin_(char *, char *, char *, char *,
53 char *, char *, char *, doublereal *, char *, char *, ftnlen,
54 ftnlen, ftnlen, ftnlen, ftnlen, ftnlen, ftnlen, ftnlen, ftnlen),
55 zzgfcosg_(doublereal *, doublereal *);
56 extern /* Subroutine */ int zzgfudlt_();
57 extern /* Subroutine */ int zzgfrelx_(U_fp, U_fp, U_fp, U_fp, S_fp, char *
58 , doublereal *, doublereal *, doublereal *, doublereal *, integer
59 *, integer *, doublereal *, logical *, U_fp, U_fp, U_fp, char *,
60 char *, logical *, L_fp, doublereal *, ftnlen, ftnlen, ftnlen);
61 integer i__;
62 extern integer cardd_(doublereal *);
63 integer n, s;
64 extern logical elemi_(integer *, integer *);
65 extern /* Subroutine */ int chkin_(char *, ftnlen), ucase_(char *, char *,
66 ftnlen, ftnlen), errch_(char *, char *, ftnlen, ftnlen), lnkan_(
67 integer *, integer *);
68 integer class__, compl;
69 logical found;
70 doublereal value;
71 integer right;
72 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen), copyd_(
73 doublereal *, doublereal *), repmi_(char *, char *, integer *,
74 char *, ftnlen, ftnlen, ftnlen);
75 integer total, f1, f2;
76 char rlist[32*7];
77 doublereal r2ovr2, start;
78 extern /* Subroutine */ int ljust_(char *, char *, ftnlen, ftnlen);
79 extern doublereal twopi_(void);
80 integer q1, q2, q3, q4;
81 extern /* Subroutine */ int bodc2s_(integer *, char *, ftnlen);
82 extern logical failed_(void);
83 extern doublereal pi_(void);
84 doublereal cv, et;
85 integer nl;
86 extern integer isrchc_(char *, integer *, char *, ftnlen, ftnlen),
87 lnknxt_(integer *, integer *), wncard_(doublereal *);
88 extern logical return_(void), smsgnd_(doublereal *, doublereal *);
89 char nrmcrd[32], nrmsys[32], prxcrd[32], prxfun[50], prxsys[32], rctrnm[
90 36], rptpre[80*2], rptsuf[80*2], tmplat[80], prxrel[6];
91 doublereal cmpval, extval, locref, loctol, prxval, sv, xrfval;
92 integer clssid, frcode, needwn[13], refctr;
93 doublereal alt, lat;
94 integer region[3], wh, bot, wwpool[26] /* was [2][13] */;
95 extern /* Subroutine */ int setmsg_(char *, ftnlen);
96 doublereal lon;
97 integer res;
98 extern /* Subroutine */ int errint_(char *, integer *, ftnlen), sigerr_(
99 char *, ftnlen), chkout_(char *, ftnlen), ssized_(integer *,
100 doublereal *), lnkini_(integer *, integer *);
101 integer top;
102 char uop[6];
103 extern /* Subroutine */ int cmprss_(char *, integer *, char *, char *,
104 ftnlen, ftnlen, ftnlen), scardd_(integer *, doublereal *);
105 integer wix[7];
106 extern /* Subroutine */ int namfrm_(char *, integer *, ftnlen), frinfo_(
107 integer *, integer *, integer *, integer *, logical *), recpgr_(
108 char *, doublereal *, doublereal *, doublereal *, doublereal *,
109 doublereal *, doublereal *, ftnlen), wninsd_(doublereal *,
110 doublereal *, doublereal *), wndifd_(doublereal *, doublereal *,
111 doublereal *), wnunid_(doublereal *, doublereal *, doublereal *),
112 lnkila_(integer *, integer *, integer *), wnintd_(doublereal *,
113 doublereal *, doublereal *), ssizei_(integer *, integer *),
114 insrti_(integer *, integer *), lnkfsl_(integer *, integer *,
115 integer *), zzgfcog_(doublereal *, doublereal *);
116 integer res1, res2;
117
118 /* $ Abstract */
119
120 /* SPICE Private routine intended solely for the support of SPICE */
121 /* routines. Users should not call this routine directly due */
122 /* to the volatile nature of this routine. */
123
124 /* This routine determines time windows when the longitude */
125 /* or right ascension of a specified vector satisfies a specified */
126 /* mathematical condition. */
127
128 /* $ Disclaimer */
129
130 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
131 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
132 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
133 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
134 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
135 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
136 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
137 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
138 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
139 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
140
141 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
142 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
143 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
144 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
145 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
146 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
147
148 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
149 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
150 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
151 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
152
153 /* $ Required_Reading */
154
155 /* SPK */
156 /* TIME */
157 /* NAIF_IDS */
158 /* FRAMES */
159
160 /* $ Keywords */
161
162 /* EPHEMERIS */
163 /* GEOMETRY */
164 /* PRIVATE */
165 /* SEARCH */
166
167 /* $ Declarations */
168 /* $ Abstract */
169
170 /* This file contains public, global parameter declarations */
171 /* for the SPICELIB Geometry Finder (GF) subsystem. */
172
173 /* $ Disclaimer */
174
175 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
176 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
177 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
178 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
179 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
180 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
181 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
182 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
183 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
184 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
185
186 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
187 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
188 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
189 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
190 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
191 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
192
193 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
194 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
195 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
196 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
197
198 /* $ Required_Reading */
199
200 /* GF */
201
202 /* $ Keywords */
203
204 /* GEOMETRY */
205 /* ROOT */
206
207 /* $ Restrictions */
208
209 /* None. */
210
211 /* $ Author_and_Institution */
212
213 /* N.J. Bachman (JPL) */
214 /* L.E. Elson (JPL) */
215 /* E.D. Wright (JPL) */
216
217 /* $ Literature_References */
218
219 /* None. */
220
221 /* $ Version */
222
223 /* - SPICELIB Version 2.0.0 29-NOV-2016 (NJB) */
224
225 /* Upgraded to support surfaces represented by DSKs. */
226
227 /* Bug fix: removed declaration of NVRMAX parameter. */
228
229 /* - SPICELIB Version 1.3.0, 01-OCT-2011 (NJB) */
230
231 /* Added NWILUM parameter. */
232
233 /* - SPICELIB Version 1.2.0, 14-SEP-2010 (EDW) */
234
235 /* Added NWPA parameter. */
236
237 /* - SPICELIB Version 1.1.0, 08-SEP-2009 (EDW) */
238
239 /* Added NWRR parameter. */
240 /* Added NWUDS parameter. */
241
242 /* - SPICELIB Version 1.0.0, 21-FEB-2009 (NJB) (LSE) (EDW) */
243
244 /* -& */
245
246 /* Root finding parameters: */
247
248 /* CNVTOL is the default convergence tolerance used by the */
249 /* high-level GF search API routines. This tolerance is */
250 /* used to terminate searches for binary state transitions: */
251 /* when the time at which a transition occurs is bracketed */
252 /* by two times that differ by no more than CNVTOL, the */
253 /* transition time is considered to have been found. */
254
255 /* Units are TDB seconds. */
256
257
258 /* NWMAX is the maximum number of windows allowed for user-defined */
259 /* workspace array. */
260
261 /* DOUBLE PRECISION WORK ( LBCELL : MW, NWMAX ) */
262
263 /* Currently no more than twelve windows are required; the three */
264 /* extra windows are spares. */
265
266 /* Callers of GFEVNT can include this file and use the parameter */
267 /* NWMAX to declare the second dimension of the workspace array */
268 /* if necessary. */
269
270
271 /* Callers of GFIDST should declare their workspace window */
272 /* count using NWDIST. */
273
274
275 /* Callers of GFSEP should declare their workspace window */
276 /* count using NWSEP. */
277
278
279 /* Callers of GFRR should declare their workspace window */
280 /* count using NWRR. */
281
282
283 /* Callers of GFUDS should declare their workspace window */
284 /* count using NWUDS. */
285
286
287 /* Callers of GFPA should declare their workspace window */
288 /* count using NWPA. */
289
290
291 /* Callers of GFILUM should declare their workspace window */
292 /* count using NWILUM. */
293
294
295 /* ADDWIN is a parameter used to expand each interval of the search */
296 /* (confinement) window by a small amount at both ends in order to */
297 /* accommodate searches using equality constraints. The loaded */
298 /* kernel files must accommodate these expanded time intervals. */
299
300
301 /* FRMNLN is a string length for frame names. */
302
303
304 /* FOVTLN -- maximum length for FOV string. */
305
306
307 /* Specify the character strings that are allowed in the */
308 /* specification of field of view shapes. */
309
310
311 /* Character strings that are allowed in the */
312 /* specification of occultation types: */
313
314
315 /* Occultation target shape specifications: */
316
317
318 /* Specify the number of supported occultation types and occultation */
319 /* type string length: */
320
321
322 /* Instrument field-of-view (FOV) parameters */
323
324 /* Maximum number of FOV boundary vectors: */
325
326
327 /* FOV shape parameters: */
328
329 /* circle */
330 /* ellipse */
331 /* polygon */
332 /* rectangle */
333
334
335 /* End of file gf.inc. */
336
337 /* $ Abstract */
338
339 /* SPICE private include file intended solely for the support of */
340 /* SPICE routines. Users should not include this routine in their */
341 /* source code due to the volatile nature of this file. */
342
343 /* This file contains private, global parameter declarations */
344 /* for the SPICELIB Geometry Finder (GF) subsystem. */
345
346 /* $ Disclaimer */
347
348 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
349 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
350 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
351 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
352 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
353 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
354 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
355 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
356 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
357 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
358
359 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
360 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
361 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
362 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
363 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
364 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
365
366 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
367 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
368 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
369 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
370
371 /* $ Required_Reading */
372
373 /* GF */
374
375 /* $ Keywords */
376
377 /* GEOMETRY */
378 /* ROOT */
379
380 /* $ Restrictions */
381
382 /* None. */
383
384 /* $ Author_and_Institution */
385
386 /* N.J. Bachman (JPL) */
387 /* E.D. Wright (JPL) */
388
389 /* $ Literature_References */
390
391 /* None. */
392
393 /* $ Version */
394
395 /* - SPICELIB Version 1.0.0, 17-FEB-2009 (NJB) (EDW) */
396
397 /* -& */
398
399 /* The set of supported coordinate systems */
400
401 /* System Coordinates */
402 /* ---------- ----------- */
403 /* Rectangular X, Y, Z */
404 /* Latitudinal Radius, Longitude, Latitude */
405 /* Spherical Radius, Colatitude, Longitude */
406 /* RA/Dec Range, Right Ascension, Declination */
407 /* Cylindrical Radius, Longitude, Z */
408 /* Geodetic Longitude, Latitude, Altitude */
409 /* Planetographic Longitude, Latitude, Altitude */
410
411 /* Below we declare parameters for naming coordinate systems. */
412 /* User inputs naming coordinate systems must match these */
413 /* when compared using EQSTR. That is, user inputs must */
414 /* match after being left justified, converted to upper case, */
415 /* and having all embedded blanks removed. */
416
417
418 /* Below we declare names for coordinates. Again, user */
419 /* inputs naming coordinates must match these when */
420 /* compared using EQSTR. */
421
422
423 /* Note that the RA parameter value below matches */
424
425 /* 'RIGHT ASCENSION' */
426
427 /* when extra blanks are compressed out of the above value. */
428
429
430 /* Parameters specifying types of vector definitions */
431 /* used for GF coordinate searches: */
432
433 /* All string parameter values are left justified, upper */
434 /* case, with extra blanks compressed out. */
435
436 /* POSDEF indicates the vector is defined by the */
437 /* position of a target relative to an observer. */
438
439
440 /* SOBDEF indicates the vector points from the center */
441 /* of a target body to the sub-observer point on */
442 /* that body, for a given observer and target. */
443
444
445 /* SOBDEF indicates the vector points from the center */
446 /* of a target body to the surface intercept point on */
447 /* that body, for a given observer, ray, and target. */
448
449
450 /* Number of workspace windows used by ZZGFREL: */
451
452
453 /* Number of additional workspace windows used by ZZGFLONG: */
454
455
456 /* Index of "existence window" used by ZZGFCSLV: */
457
458
459 /* Progress report parameters: */
460
461 /* MXBEGM, */
462 /* MXENDM are, respectively, the maximum lengths of the progress */
463 /* report message prefix and suffix. */
464
465 /* Note: the sum of these lengths, plus the length of the */
466 /* "percent complete" substring, should not be long enough */
467 /* to cause wrap-around on any platform's terminal window. */
468
469
470 /* Total progress report message length upper bound: */
471
472
473 /* End of file zzgf.inc. */
474
475 /* $ Brief_I/O */
476
477 /* VARIABLE I/O DESCRIPTION */
478 /* -------- --- -------------------------------------------------- */
479 /* LBCELL P Cell lower bound. */
480 /* VECDEF I Vector definition. */
481 /* METHOD I Computation method. */
482 /* TARGET I Target name. */
483 /* REF I Reference frame name. */
484 /* ABCORR I Aberration correction. */
485 /* OBSRVR I Observer name. */
486 /* DREF I Ray's direction vector frame. */
487 /* DVEC I Ray's direction vector. */
488 /* CRDSYS I Coordinate system name. */
489 /* CRDNAM I Coordinate name. */
490 /* RELATE I Relational operator. */
491 /* REFVAL I Reference value. */
492 /* TOL I Convergence tolerance. */
493 /* ADJUST I Absolute extremum adjustment value. */
494 /* UDSTEP I Step size routine. */
495 /* UDREFN I Search refinement routine. */
496 /* RPT I Progress report flag. */
497 /* UDREPI I Progress report initialization routine. */
498 /* UDREPU I Progress report update routine. */
499 /* UDREPF I Progress report termination routine. */
500 /* BAIL I Bail-out flag. */
501 /* UDBAIL I Bail-out status function. */
502 /* MW I Workspace window size. */
503 /* NW I Workspace window count. */
504 /* WORK I-O Workspace window array. */
505 /* CNFINE I Confinement window. */
506 /* RESULT O Result window. */
507
508 /* $ Detailed_Input */
509
510
511 /* VECDEF Every coordinate computed by this routine is a */
512 /* function of an underlying vector. VECDEF is a short */
513 /* string describing the means by which the vector of */
514 /* interest is defined. Only parameters from the Fortran */
515 /* INCLUDE file zzgf.inc should be used. Parameter names */
516 /* and meanings are: */
517
518 /* POSDEF Vector is position of */
519 /* target relative to observer. */
520
521 /* SOBDEF Vector is sub-observer */
522 /* point on target body. Vector */
523 /* points from target body */
524 /* center to sub-observer point. */
525 /* The target must be an extended */
526 /* body modeled as a triaxial */
527 /* ellipsoid. */
528
529 /* SINDEF Vector is ray-surface intercept */
530 /* point on target body. Vector */
531 /* points from target body */
532 /* center to sub-observer point. */
533 /* The target must be an extended */
534 /* body modeled as a triaxial */
535 /* ellipsoid. */
536
537 /* Case, leading and trailing blanks ARE significant */
538 /* in the string VECDEF. */
539
540
541 /* METHOD is a string specifying the computational method */
542 /* applicable to the vector of interest. When VECDEF */
543 /* is the parameter */
544
545 /* SOBDEF */
546
547 /* METHOD should be set to one of the values accepted */
548 /* by the SPICELIB routine SUBPNT. */
549
550 /* When VECDEF is the parameter */
551
552 /* SINDEF */
553
554 /* METHOD should be set to one of the values accepted */
555 /* by the SPICELIB routine SINCPT. */
556
557 /* METHOD is ignored if VECDEF is set to */
558
559 /* POSDEF */
560
561 /* Case, leading and trailing blanks are not significant */
562 /* in the string METHOD. */
563
564
565 /* TARGET is the name of the target object. */
566
567
568 /* REF is the name of the reference frame relative to which */
569 /* the vector of interest is specified. The specified */
570 /* condition applies to the specified coordinate of */
571 /* of this vector in frame REF. */
572
573 /* When geodetic coordinates are used, the reference */
574 /* ellipsoid is assumed to be that associated with */
575 /* the central body of the frame designated by REF. */
576 /* In this case, the central body of the frame must */
577 /* be an extended body. */
578
579 /* Case, leading and trailing blanks are not significant */
580 /* in the string REF. */
581
582
583 /* ABCORR indicates the aberration corrections to be applied to */
584 /* the state of the target body to account for one-way */
585 /* light time and stellar aberration. The orientation */
586 /* of the target body will also be corrected for one-way */
587 /* light time when light time corrections are requested. */
588
589 /* Supported aberration correction options for */
590 /* observation (case where radiation is received by */
591 /* observer at ET) are: */
592
593 /* 'NONE' No correction. */
594 /* 'LT' Light time only. */
595 /* 'LT+S' Light time and stellar aberration. */
596 /* 'CN' Converged Newtonian (CN) light time. */
597 /* 'CN+S' CN light time and stellar aberration. */
598
599 /* Supported aberration correction options for */
600 /* transmission (case where radiation is emitted from */
601 /* observer at ET) are: */
602
603 /* 'XLT' Light time only. */
604 /* 'XLT+S' Light time and stellar aberration. */
605 /* 'XCN' Converged Newtonian (CN) light time. */
606 /* 'XCN+S' CN light time and stellar aberration. */
607
608 /* For detailed information, see the geometry finder */
609 /* required reading, gf.req. Also see the header of */
610 /* SPKEZR, which contains a detailed discussion of */
611 /* aberration corrections. */
612
613 /* Case, leading and trailing blanks are not significant */
614 /* in the string ABCORR. */
615
616
617 /* OBSRVR is the name of the observer. */
618
619
620 /* DREF is the name of the reference frame relative to which a */
621 /* ray's direction vector is expressed. This may be any */
622 /* frame supported by the SPICE system, including */
623 /* built-in frames (documented in the Frames Required */
624 /* Reading) and frames defined by a loaded frame kernel */
625 /* (FK). The string DREF is case-insensitive, and leading */
626 /* and trailing blanks in DREF are not significant. */
627
628 /* When DREF designates a non-inertial frame, the */
629 /* orientation of the frame is evaluated at an epoch */
630 /* dependent on the frame's center and, if the center is */
631 /* not the observer, on the selected aberration */
632 /* correction. See the description of the direction */
633 /* vector DVEC for details. */
634
635
636 /* DVEC Ray direction vector emanating from the observer. The */
637 /* intercept with the target body's surface of the ray */
638 /* defined by the observer and DVEC is sought. */
639
640 /* DVEC is specified relative to the reference frame */
641 /* designated by DREF. */
642
643 /* Non-inertial reference frames are treated as follows: */
644 /* if the center of the frame is at the observer's */
645 /* location, the frame is evaluated at ET. If the frame's */
646 /* center is located elsewhere, then letting LTCENT be */
647 /* the one-way light time between the observer and the */
648 /* central body associated with the frame, the */
649 /* orientation of the frame is evaluated at ET-LTCENT, */
650 /* ET+LTCENT, or ET depending on whether the requested */
651 /* aberration correction is, respectively, for received */
652 /* radiation, transmitted radiation, or is omitted. */
653 /* LTCENT is computed using the method indicated by */
654 /* ABCORR. */
655
656
657 /* CRDSYS is the name of the coordinate system to which the */
658 /* coordinate of interest belongs. Allowed values are */
659 /* those defined in the GF Fortran INCLUDE file */
660
661 /* zzgf.inc. */
662
663 /* CRDSYS must refer to a system in which longitude */
664
665 /* or right ascension is a coordinate. Note that when */
666 /* geodetic coordinates are used, the reference ellipsoid */
667 /* is that associated with the central body of the */
668 /* reference frame designated by REF. The central body */
669 /* must be an extended body in this case. */
670
671 /* Case, leading and trailing blanks are not significant */
672 /* in the string CRDSYS. */
673
674
675 /* CRDNAM is the name of the coordinate of interest: this is */
676 /* the coordinate to which the specified condition */
677 /* applies. Supported coordinates are */
678
679 /* Planetocentric longitude */
680 /* Right ascension */
681
682 /* which are designated respectively by the parameters */
683
684 /* LONCRD */
685 /* RACRD */
686
687 /* See the INCLUDE file */
688
689 /* zzgf.inc */
690
691 /* for the declarations of these parameters. */
692
693 /* For the */
694
695 /* Latitudinal */
696 /* Geodetic */
697 /* Spherical */
698
699 /* coordinate systems, longitude lies in the range */
700
701 /* ( -pi, pi ] */
702
703 /* For the */
704
705 /* Cylindrical */
706 /* Planetographic */
707
708 /* coordinate systems, longitude lies in the range */
709
710 /* [ 0, 2*pi ) */
711
712 /* Right ascension lies in the range */
713
714 /* [ 0, 2*pi ) */
715
716 /* Case, leading and trailing blanks are not significant */
717 /* in the string CRDNAM. */
718
719
720 /* RELATE is a relational operator used to define a constraint */
721 /* on longitude or right ascension of the specified */
722 /* vector. The result window found by this routine */
723 /* indicates the time intervals where the constraint is */
724 /* satisfied. Supported values of RELATE and */
725 /* corresponding meanings are shown below: */
726
727 /* '>' Longitude or RA is greater than the */
728 /* reference value REFVAL. */
729
730 /* '=' Longitude or RA is equal to the reference */
731 /* value REFVAL. */
732
733 /* '<' Longitude or RA is less than the */
734 /* reference value REFVAL. */
735
736
737 /* 'ABSMAX' Longitude or RA is at an absolute maximum. */
738
739 /* 'ABSMIN' Longitude or RA is at an absolute */
740 /* minimum. */
741
742 /* 'LOCMAX' Longitude or RA is at a local maximum. */
743
744 /* 'LOCMIN' Longitude or RA is at a local minimum. */
745
746 /* The caller may indicate that the region of interest */
747 /* is the set of time intervals where the quantity is */
748 /* within a specified tolerance of an absolute extremum. */
749 /* The argument ADJUST (described below) is used to */
750 /* specify this tolerance. */
751
752 /* Local extrema are considered to exist only in the */
753 /* interiors of the intervals comprising the confinement */
754 /* window: a local extremum cannot exist at a boundary */
755 /* point of the confinement window. */
756
757 /* Case is not significant in the string RELATE. */
758
759
760 /* REFVAL is the reference value used to define equality or */
761 /* inequality conditions. */
762
763 /* REFVAL has units of radians. */
764
765 /* When the coordinate of interest is longitude, REFVAL */
766 /* is interpreted as though it were translated, if */
767 /* necessary, by an integer multiple of 2*pi to place it */
768 /* in the standard range for longitude: (-pi, pi]. */
769 /* Similarly, when the coordinate of interest is right */
770 /* ascension, REFVAL is interpreted as though it were */
771 /* translated, if necessary, by an integer multiple of */
772 /* 2*pi into the range [0, 2*pi). */
773
774 /* Example: suppose REFVAL is set to -4.5. Then the */
775 /* condition */
776
777 /* longitude equals REFVAL */
778
779 /* is interpreted as */
780
781 /* longitude equals -0.5 * pi */
782
783 /* so the solution window for this condition may well */
784 /* be non-empty. */
785
786 /* REFVAL is ignored if RELATE is not an equality or */
787 /* inequality operator. */
788
789
790 /* TOL is a tolerance value used to determine convergence of */
791 /* root-finding operations. TOL is measured in TDB */
792 /* seconds and is greater than zero. */
793
794
795 /* ADJUST The amount by which the coordinate is allowed to vary */
796 /* from an absolute extremum. ADJUST is not used for */
797 /* equality or inequality conditions. ADJUST must not be */
798 /* negative. */
799
800 /* If ADJUST is positive and a search for an absolute */
801 /* minimum is performed, the resulting schedule contains */
802 /* time intervals when the specified coordinate has */
803 /* values between */
804
805 /* ABSMIN */
806 /* and MIN ( ABSMIN + ADJUST, MX ) */
807
808 /* where MX is the maximum value of the coordinate's */
809 /* range. */
810
811 /* If the search is for an absolute maximum, the */
812 /* corresponding range is between */
813
814 /* MAX ( ABSMAX - ADJUST, MN ) */
815 /* and ABSMAX */
816
817 /* where MN is the minimum value of the coordinate's */
818 /* range. */
819
820
821 /* UDSTEP is a routine that computes a time step used to search */
822 /* for a transition of the state of the specified */
823 /* coordinate. In the context of this routine's */
824 /* algorithm, a "state transition" occurs where the */
825 /* coordinate's time derivative changes from negative to */
826 /* non-negative or vice versa. */
827
828 /* This routine relies on UDSTEP returning step sizes */
829 /* small enough so that state transitions within the */
830 /* confinement window are not overlooked. There must */
831 /* never be two roots A and B separated by less than */
832 /* STEP, where STEP is the minimum step size returned by */
833 /* UDSTEP for any value of ET in the interval [A, B]. */
834
835 /* The calling sequence for UDSTEP is: */
836
837 /* CALL UDSTEP ( ET, STEP ) */
838
839 /* where: */
840
841 /* ET is the input start time from which the */
842 /* algorithm is to search forward for a state */
843 /* transition. ET is expressed as seconds past */
844 /* J2000 TDB. ET is a DOUBLE PRECISION number. */
845
846 /* STEP is the output step size. STEP indicates */
847 /* how far to advance ET so that ET and */
848 /* ET+STEP may bracket a state transition and */
849 /* definitely do not bracket more than one */
850 /* state transition. STEP is a DOUBLE */
851 /* PRECISION number. Units are TDB seconds. */
852
853 /* If a constant step size is desired, the routine GFSTEP */
854 /* may be used. GFSTEP returns the step size that was set */
855 /* via the most recent call to GFSSTP. */
856
857
858 /* UDREFN is the name of the externally specified routine that */
859 /* computes a refinement in the times that bracket a */
860 /* transition point. In other words, once a pair of */
861 /* times have been detected such that the system is in */
862 /* different states at each of the two times, UDREFN */
863 /* selects an intermediate time which should be closer to */
864 /* the transition state than one of the two known times. */
865 /* The calling sequence for UDREFN is: */
866
867 /* CALL UDREFN ( T1, T2, S1, S2, T ) */
868
869 /* where the inputs are: */
870
871 /* T1 is a time when the system is in state S1. T1 */
872 /* is a DOUBLE PRECISION number. */
873
874 /* T2 is a time when the system is in state S2. T2 */
875 /* is a DOUBLE PRECISION number and is assumed */
876 /* to be larger than T1. */
877
878 /* S1 is the state of the system at time T1. */
879 /* S1 is a LOGICAL value. */
880
881 /* S2 is the state of the system at time T2. */
882 /* S2 is a LOGICAL value. */
883
884 /* UDREFN may use or ignore the S1 and S2 values. */
885
886 /* The output is: */
887
888 /* T is next time to check for a state transition. */
889 /* T is a DOUBLE PRECISION number between T1 and */
890 /* T2. */
891
892 /* If a simple bisection method is desired, the routine */
893 /* GFREFN may be used. This is the default option. */
894
895
896 /* RPT is a logical variable which controls whether the */
897 /* progress reporter is on or off; setting RPT */
898 /* to .TRUE. enables progress reporting. */
899
900
901 /* UDREPI is a user-defined subroutine that initializes a */
902 /* progress report. When progress reporting is */
903 /* enabled, UDREPI is called at the start of a search */
904 /* pass (see the implementation of ZZGFREL for details on */
905 /* search passes). The calling sequence of UDREPI is */
906
907 /* UDREPI ( CNFINE, RPTPRE, RPTSUF ) */
908
909 /* DOUBLE PRECISION CNFINE ( LBCELL : * ) */
910 /* CHARACTER*(*) RPTPRE */
911 /* CHARACTER*(*) RPTSUF */
912
913 /* where */
914
915 /* CNFINE */
916
917 /* is the confinement window passed into ZZGFREL, and */
918
919 /* RPTPRE */
920 /* RPTSUF */
921
922 /* are prefix and suffix strings used in the progress */
923 /* report: these strings are intended to bracket a */
924 /* representation of the fraction of work done. */
925
926 /* SPICELIB provides the default progress reporting */
927 /* initialization routine GFREPI. If GFREPI is used, then */
928 /* the progress reporting update and termination routines */
929 /* GFREPU and GFREPF must be used as well. */
930
931
932 /* UDREPU is a user-defined subroutine that updates the */
933 /* progress report for a search pass. The calling */
934 /* sequence of UDREPU is */
935
936 /* UDREPU (IVBEG, IVEND, ET ) */
937
938 /* DOUBLE PRECISION ET */
939 /* DOUBLE PRECISION IVBEG */
940 /* DOUBLE PRECISION IVEND */
941
942 /* where ET is an epoch belonging to the confinement */
943 /* window, IVBEG and IVEND are the start and stop times, */
944 /* respectively of the current confinement window */
945 /* interval. The ratio of the measure of the portion */
946 /* of CNFINE that precedes ET to the measure of CNFINE */
947 /* would be a logical candidate for the search's */
948 /* completion percentage; however the method of */
949 /* measurement is up to the user. */
950
951
952 /* UDREPF is a user-defined subroutine that finalizes a */
953 /* progress report. UDREPF has no arguments. */
954
955
956 /* BAIL is a logical flag indicating whether or not interrupt */
957 /* signal handling is enabled. Setting BAIL to .TRUE. */
958 /* enables interrupt signal handling: the GF system will */
959 /* then call UDBAIL to check for interrupt signals. */
960
961
962 /* UDBAIL is the name of a user defined logical function that */
963 /* checks to see whether an interrupt signal has been */
964 /* issued from, e.g. the keyboard. UDBAIL is used only */
965 /* when BAIL is set to .TRUE. If interrupt handling is */
966 /* not used, the SPICELIB function GFBAIL should be */
967 /* passed in as the actual bail-out function argument. */
968
969
970 /* MW is the cell size of the windows in the workspace array */
971 /* WORK. */
972
973
974 /* NW is the number of windows in the workspace array WORK. */
975 /* NW must be at least as large as the parameter NWMAX. */
976
977
978 /* WORK is an array used to store workspace windows. This */
979 /* array has dimensions ( LBCELL : MW, NW). */
980
981 /* The windows contained WORK that used by this routine */
982 /* are initialized here to have size MW. The other */
983 /* elements of WORK are not modified. */
984
985
986 /* CNFINE is a SPICE window that confines the bounds of the */
987 /* search. */
988
989 /* For coordinates defined by ray-target surface */
990 /* intercepts, the intercept and its time derivative are */
991 /* expected to be computable on the confinement window. */
992
993
994 /* RESULT is an initialized SPICE window. RESULT must be large */
995 /* enough to hold all of the intervals, within the */
996 /* confinement window, on which the specified condition */
997 /* is met. */
998
999 /* RESULT must be initialized by the caller via the */
1000 /* SPICELIB routine SSIZED. */
1001
1002 /* $ Detailed_Output */
1003
1004
1005 /* WORK has undefined contents on output, with the exception */
1006 /* of the windows occupying the range */
1007
1008 /* ( LBCELL : NW, EXWIDX : NWMAX ) */
1009
1010 /* which are not modified by this routine. */
1011
1012 /* RESULT is a SPICELIB window containing the intersection of */
1013 /* the confinement window and the set of time intervals */
1014 /* when the value of the specified coordinate satisfies */
1015 /* constraints specified by RELATE and ADJUST. */
1016
1017 /* For coordinates defined by ray-target surface */
1018 /* intercepts, RESULT is further restricted to the window */
1019 /* over which the intercept and its derivative with */
1020 /* respect to time are computable. See the description of */
1021 /* CNFINE above for details. */
1022
1023 /* Due to computational accuracy limitations, the */
1024 /* coordinate of interest *may not satisfy the */
1025 /* specified condition* at all points belonging to */
1026 /* RESULT. For example, if the caller specifies */
1027 /* a tolerance of 1.E-6 seconds and seeks the */
1028 /* solution set for the condition */
1029
1030 /* The planetocentric longitude of the geometric */
1031 /* earth-sun vector in the J2000 frame is greater */
1032 /* than or equal to zero */
1033
1034 /* the right endpoints of some intervals in RESULT may be */
1035 /* times that map to negative longitude values very close */
1036 /* to -pi radians. */
1037
1038 /* The user (of SPICE API routines dependent on this */
1039 /* routine) may wish to contract RESULT using WNCOND in */
1040 /* order to guarantee that the specified condition */
1041 /* is satisfied on RESULT. Selection of a suitable */
1042 /* contraction value is dependent on the user's */
1043 /* requirements and the specific problem to be solved. */
1044
1045 /* $ Parameters */
1046
1047 /* LBCELL is the SPICELIB cell lower bound. */
1048
1049 /* $ Exceptions */
1050
1051 /* 1) In order for this routine to produce correct results, */
1052 /* the external step size routine UDGSTP must return step sizes */
1053 /* appropriate for the problem at hand. Step sizes that */
1054 /* are too large may cause this routine to miss roots; */
1055 /* step sizes that are too small may cause this routine to */
1056 /* run unacceptably slowly and in some cases, find spurious */
1057 /* roots. */
1058
1059 /* This routine does not diagnose invalid step sizes, */
1060 /* except that if the step size is non-positive, the error */
1061 /* SPICE(VALUEOUTOFRANGE) is signaled. */
1062
1063 /* 2) In order for this routine to produce correct results, */
1064 /* the convergence tolerance TOL must be appropriate for the */
1065 /* problem at hand. The error in any interval endpoint */
1066 /* contained in RESULT should be expected to be no smaller */
1067 /* than TOL; depending on the behavior of the coordinate */
1068 /* and the condition, the error could be much larger. For */
1069 /* example, for some functions, finding correct, unique */
1070 /* extrema is notoriously difficult. */
1071
1072 /* The user should keep in mind that the minimum separation */
1073 /* between successive values of ET is about 1.E-7 seconds */
1074 /* for SPICE platforms and values of ET not extremely close to */
1075 /* J2000. */
1076
1077 /* This routine does not diagnose invalid tolerance values, */
1078 /* except that if the tolerance is non-positive, the error */
1079 /* SPICE(VALUEOUTOFRANGE) is signaled. */
1080
1081 /* 3) A negative value for ADJUST causes the routine to signal */
1082 /* the error SPICE(VALUEOUTOFRANGE). A non-zero value for ADJUST */
1083 /* when RELATE has any value other than "ABSMIN" or "ABSMAX", */
1084 /* causes the routine to signal the error SPICE(INVALIDVALUE). */
1085
1086 /* 4) If the operator string RELATE doesn't contain a recognized */
1087 /* value, the error SPICE(NOTRECOGNIZED) is signaled. */
1088
1089 /* 5) If any error occurs while initializing the coordinate */
1090 /* utility package, the error will be diagnosed by routines */
1091 /* in the call tree of ZZGFCOIN. */
1092
1093 /* 6) If any error occurs while performing computations */
1094 /* to determine if a quantity of interest is decreasing */
1095 /* at a specified time, the error will be diagnosed by */
1096 /* routines in the call tree of this routine. */
1097
1098 /* 7) If any error occurs while performing computations */
1099 /* to determine if a quantity of interest is less than a */
1100 /* specified reference value at a specified time, the error will */
1101 /* be diagnosed by routines in the call tree of this routine. */
1102
1103 /* 8) If an error (typically cell overflow) occurs while performing */
1104 /* window arithmetic, the error will be diagnosed by */
1105 /* routines in the call trees of window routines called by */
1106 /* this routine. */
1107
1108 /* 9) Due to numerical errors, in particular, */
1109
1110 /* - Truncation error in time values */
1111 /* - Finite tolerance value */
1112 /* - Errors in computed geometric quantities */
1113
1114 /* it is *normal* that the condition of interest is not */
1115 /* satisfied on the entire result window. */
1116
1117 /* The result window may need to be contracted slightly by the */
1118 /* caller to achieve desired results, in particular to remove */
1119 /* times where discontinuities of longitude or right ascension */
1120 /* are crossed. */
1121
1122 /* 10) Most relational conditions involving longitude or */
1123 /* right ascension make sense only when latitude or declination */
1124 /* is bounded away from +/- pi/2 radians. Users should */
1125 /* select the confinement window accordingly. */
1126
1127 /* 11) The user must take care when searching for an extremum */
1128 /* (ABSMAX, ABSMIN, LOCMAX, LOCMIN) of LONGITUDE or */
1129 /* RIGHT ASCENSION values. Since these quantities are cyclical, */
1130 /* rather than monotonically increasing or decreasing, an */
1131 /* extremum may be hard to interpret. In particular, if an */
1132 /* extremum is found near the cycle boundary (- PI for */
1133 /* longitude, 2 PI for RIGHT ASCENSION) it may not be */
1134 /* numerically reasonable. For example, the search for times */
1135 /* when a longitude coordinate is at its absolute maximum may */
1136 /* result in a time when the longitude value is - PI, due to */
1137 /* roundoff error. */
1138
1139 /* $ Files */
1140
1141 /* This routine doesn't directly participate in SPICE kernel loading */
1142 /* or unloading. However, a variety of SPICE kernels must be loaded */
1143 /* in order for this routine to work: */
1144
1145 /* - Since all coordinate computations supported by this routine */
1146 /* depend on observer-target vectors, at a minimum, SPK files */
1147 /* providing ephemeris data enabling computation of these */
1148 /* vectors are required. */
1149
1150 /* - If non-inertial reference frames are used, then PCK */
1151 /* files, frame kernels, C-kernels, and SCLK kernels may be */
1152 /* needed. */
1153
1154 /* - If the coordinate of interest is defined in terms of a target */
1155 /* surface point, then (currently) a PCK providing radii for a */
1156 /* triaxial shape model must be loaded. */
1157
1158 /* - If geodetic coordinates are used, then a PCK providing radii */
1159 /* for a triaxial shape model must be loaded. */
1160
1161 /* See the Files section of GFEVNT's header for further information. */
1162
1163
1164 /* $ Particulars */
1165
1166 /* Since this is a private SPICELIB routine, the header comments */
1167 /* make many references to the routine's implementation. This */
1168 /* is done to help the maintenance programmer understand the */
1169 /* routine; however, these comments may themselves need to be */
1170 /* updated if the GF subsystem implementation changes. */
1171
1172 /* This routine determines time windows when the longitude or right */
1173 /* ascension of a specified vector satisfies a specified */
1174 /* mathematical condition. This routine can (in some cases, by */
1175 /* means of multiple calls) answer questions such as */
1176
1177 /* When does the moon pass over the earth's prime meridian? */
1178
1179 /* Given a time window when the geodetic latitude of the MGS */
1180 /* spacecraft relative to the IAU_MARS frame is between -30 : +30 */
1181 /* degrees, when within this window is the planetographic */
1182 /* longitude of the spacecraft between 15 and 16 degrees? */
1183
1184 /* For brevity, throughout this routine, we'll refer to the vector */
1185 /* whose longitude or right ascension is of interest as "the vector" */
1186 /* or "the vector of interest." We'll also call the longitude or */
1187 /* right ascension "the coordinate" or "the coordinate of interest." */
1188
1189 /* A note concerning processing speed: the algorithm used by this */
1190 /* routine takes a "divide and conquer" approach that involves, in */
1191 /* many cases, multiple calls to the low-level GF root finding */
1192 /* routines. So the user can expect most longitude or right */
1193 /* ascension computations to be relatively slow. Using a */
1194 /* confinement window that is more easily computed, say one */
1195 /* involving latitude constraints, can be very helpful. */
1196
1197 /* $ Examples */
1198
1199 /* See usage in GFEVNT. */
1200
1201 /* $ Restrictions */
1202
1203 /* 1) The interface and functionality of this routine may change */
1204 /* without notice. This routine should be called only by */
1205 /* SPICELIB routines. */
1206
1207 /* 2) Root-finding problems of the sort solved by this routine are, */
1208 /* when a computer is involved, replete with mathematical */
1209 /* complications. We've tried to cover all the angles in the */
1210 /* Detailed_Input, Detailed_Output, and Exceptions header */
1211 /* sections. No doubt some issues remain unaddressed. Correct */
1212 /* usage of this routine depends in good measure on the user */
1213 /* posing "reasonable" problems to solve. */
1214
1215 /* 3) The kernel files to be used by ZZGFLONG must be loaded */
1216 /* (normally via the SPICELIB routine FURNSH) before ZZGFLONG is */
1217 /* called. */
1218
1219 /* 4) This routine has the side effect of re-initializing the */
1220 /* coordinate quantity utility package. Callers may themselves */
1221 /* need to re-initialize the coordinate quantity utility package */
1222 /* after calling this routine. */
1223
1224 /* $ Literature_References */
1225
1226 /* None. */
1227
1228 /* $ Author_and_Institution */
1229
1230 /* E.D. Wright (JPL) */
1231 /* N.J. Bachman (JPL) */
1232
1233 /* $ Version */
1234
1235 /* - SPICELIB Version 2.1.0 04-APR-2011 (EDW) */
1236
1237 /* Replaced use of rooutines ZZGFREL with ZZGFRELX. ZZGFCOIN */
1238 /* argument list edited to remove the unneeded argument REFVAL. */
1239
1240 /* The 2.1.0 changes should not affect the numerical results */
1241 /* of GF coordinate computations. */
1242
1243 /* - SPICELIB Version 2.0.0 12-MAY-2009 (NJB) */
1244
1245 /* Upgraded to support targets and observers having */
1246 /* no names associated with their ID codes. */
1247
1248 /* - SPICELIB Version 1.0.0 23-FEB-2009 (NJB) (EDW) */
1249
1250 /* -& */
1251
1252 /* SPICELIB functions */
1253
1254
1255 /* Entry points in the coordinate utility package. */
1256 /* We have the usual GF entry points for the coordinate, plus */
1257 /* utilities for the cosine and sine of the coordinate. */
1258
1259 /* Names and meanings: */
1260
1261 /* ZZGFCODC Is coordinate decreasing? */
1262 /* ZZGFCOG Get coordinate value. */
1263 /* ZZGFCOCD Is cosine of the coordinate decreasing? */
1264 /* ZZGFCOCG Get cosine of the coordinate value. */
1265 /* ZZGFCOSD Is sine of the coordinate decreasing? */
1266 /* ZZGFCOSG Get sine of the coordinate value. */
1267
1268
1269 /* The routine to test UDFUNC < REFVAL. */
1270
1271
1272 /* Local parameters */
1273
1274
1275
1276 /* Margin for branch cut avoidance. Units are radians: */
1277
1278
1279 /* Margin for local extrema search. Units are radians: */
1280
1281
1282 /* Short alias for LBCELL: */
1283
1284
1285 /* Number of supported comparison operators: */
1286
1287
1288 /* Assorted string lengths: */
1289
1290 /* Maximum body name length: */
1291
1292
1293 /* NAMLEN is the maximum length of both a frame name and of */
1294 /* any kernel pool variable name. */
1295
1296
1297 /* OPLEN is the maximum string length for comparison operators. */
1298 /* OPLEN may grow if new comparisons are added. */
1299
1300
1301 /* FUNLEN is the length of the function name string. */
1302
1303
1304 /* CRDLEN is the maximum length of a coordinate name. */
1305
1306
1307 /* SYSLEN is the maximum length of a coordinate system name. */
1308
1309
1310 /* RPTLEN is the maximum length of a progress reporter message. */
1311
1312
1313 /* Local variables */
1314
1315
1316 /* Saved variables */
1317
1318
1319 /* Initial values */
1320
1321 /* Parameter adjustments */
1322 work_dim1 = *mw + 6;
1323 work_dim2 = *nw;
1324 work_offset = work_dim1 - 5;
1325
1326 /* Function Body */
1327
1328 /* Standard SPICE error handling. */
1329
1330 if (return_()) {
1331 return 0;
1332 } else {
1333 chkin_("ZZGFLONG", (ftnlen)8);
1334 }
1335
1336 /* Overview */
1337 /* ======== */
1338
1339
1340 /* Terminology */
1341 /* ----------- */
1342
1343 /* - Proxy function */
1344
1345 /* In many cases, instead of finding a time window */
1346 /* where the coordinate of interest satisfies a specified */
1347 /* condition, we'll find a time window where a second, related */
1348 /* function satisfies a related condition. We'll call this */
1349 /* second function the "proxy function." */
1350
1351 /* The proxy function will be one that is "better behaved" */
1352 /* than the original in the domain of interest. For */
1353 /* example, when searching for times when longitude is */
1354 /* equal to pi radians, we may instead intersect the */
1355 /* confinement window with a window on which cosine of */
1356 /* longitude is negative, and then within that more */
1357 /* restricted intersection, find the times when the sine */
1358 /* of longitude is zero. In this example sine(longitude) */
1359 /* is a proxy function for longitude. */
1360
1361 /* - Resolution of a function */
1362
1363 /* Below we'll refer to the "resolution" of a proxy function. */
1364 /* In order to find roots accurately, it's necessary for */
1365 /* a proxy function to change a by a reasonable amount */
1366 /* when the function it represents changes. Mathematically, */
1367 /* the magnitude of the derivative of the proxy function */
1368 /* with respect to the function it represents should not */
1369 /* be too much less than 1. An example of a *bad* choice */
1370 /* of a proxy function would be to use cosine of longitude */
1371 /* as a proxy function for longitude in a confinement */
1372 /* window in which longitude is close to zero. This */
1373 /* choice would lead to considerable loss of accuracy. On */
1374 /* the other hand, sine of longitude would be a reasonable */
1375 /* proxy function for this case. */
1376
1377 /* - The unit circle */
1378
1379 /* In the discussion below, we'll freely associate angular */
1380 /* coordinates with locations on the unit circle. For example, */
1381 /* we might say "longitude is in the upper half of the unit */
1382 /* circle." */
1383
1384 /* - Window aliases */
1385
1386 /* We're going to make extensive use workspace windows. */
1387 /* In many cases, we'll need to reuse various windows for */
1388 /* different purposes at different times. So instead */
1389 /* of using mnemonic parameter names for window indices, */
1390 /* we'll use variables we call window aliases. For example, */
1391 /* when we want to use the 8th workspace window to hold */
1392 /* the window of times when longitude is in the upper half */
1393 /* of the unit circle, we'll set the alias UPPER equal to */
1394 /* 8, so we can refer to the window by */
1395
1396 /* WORK( LB, UPPER ) */
1397
1398 /* and keep track of what we're using the window for. */
1399
1400 /* Some of the aliases aren't wonderful names: we use */
1401 /* F1, F2, etc. to represent "free" window 1, 2, etc. */
1402
1403
1404 /* Algorithm */
1405 /* --------- */
1406
1407 /* - Equality */
1408
1409 /* We use sine or cosine of the coordinate as proxy functions. */
1410 /* The proxy function having the better resolution is */
1411 /* selected. For example, to find times when right ascension */
1412 /* is 2*pi/3, we search for the times when cosine of right */
1413 /* ascension is equal to -1/2. Since these searches can produce */
1414 /* spurious roots, we cast out any such roots after completing */
1415 /* the search. */
1416
1417
1418 /* - Local extrema */
1419
1420 /* We first find local extrema in the right and left half */
1421 /* circles, using longitude as a proxy function on the right */
1422 /* half and right ascension on the left. */
1423
1424
1425 /* - Absolute extrema */
1426
1427 /* We deal with absolute extrema before inequalities because */
1428 /* this allows us to use the code (later in this routine) for */
1429 /* inequality relations when the user specifies a non-zero */
1430 /* ADJUST value. When ADJUST is non-zero, having the actual */
1431 /* extreme value in hand, we can easily solve for the window */
1432 /* in which the coordinate is greater than */
1433
1434 /* <absolute maximum> - ADJUST */
1435
1436 /* or less than */
1437
1438 /* <absolute minimum> + ADJUST */
1439
1440 /* Below, "Searching in a region" means that we find the */
1441 /* window when the coordinate is in the region (and of course */
1442 /* in the confinement window), then use this window as the */
1443 /* confinement window. */
1444
1445 /* Finding absolute extrema is a matter of successively */
1446 /* searching for extrema in different parts of the unit */
1447 /* circle. For example, when we search for an absolute */
1448 /* maximum of longitude, we first search in the second */
1449 /* quadrant, then if we find nothing, the right half circle, */
1450 /* then if we find nothing, the fourth quadrant. */
1451
1452 /* We always use longitude as a proxy function on the right */
1453 /* half circle and right ascension as a proxy function on */
1454 /* the left half circle. */
1455
1456
1457 /* - Inequality */
1458
1459 /* In general, we use proxy functions and break up the unit */
1460 /* circle into regions where the proxy functions are single */
1461 /* valued. The exact solution approach depends on where the */
1462 /* reference value is. For example, to find the window on */
1463 /* which longitude is less than 3*pi/4, we first search */
1464 /* for the solution in the second quadrant. We then */
1465 /* combine this result window with the window of times */
1466 /* when longitude is in the right half circle, and with */
1467 /* the window of times when longitude is in the third */
1468 /* quadrant. */
1469
1470
1471 /* Code layout */
1472 /* ----------- */
1473
1474 /* We've tried to arrange the code to minimize calls to */
1475 /* ZZGFREL, primarily because these expensive in terms of */
1476 /* run time. They also take up a lot of space. */
1477
1478 /* The code starts out by re-formulating the constraint, */
1479 /* if necessary, as one applying to planetocentric longitude */
1480 /* or right ascension. This simplifies the subsequent logic. */
1481
1482 /* Equality searches are handled before the rest. The routine */
1483 /* exits after processing a search having an equality constraint. */
1484
1485 /* Searches for local extrema are handled next. Again, the */
1486 /* routine exits after processing these types of searches. */
1487
1488 /* The next portion of the code is devoted to dealing with */
1489 /* absolute extrema. If the search is for absolute extrema and */
1490 /* AJDUST is non-zero, we use the results from this portion of */
1491 /* the code to set up an inequality search, which is done below. */
1492
1493 /* After the portion of the code dealing with absolute extrema */
1494 /* with ADJUST equal to zero, we perform setup functions to */
1495 /* prepare to call ZZGFREL. In general, what's happening here is */
1496 /* that we're deciding what regions of the unit circle we're */
1497 /* going to use in our solution, and we prepare to find windows */
1498 /* when the coordinate is in the various regions of interest. */
1499 /* This setup code includes assignment of window aliases, */
1500 /* selection of proxy functions, and setting flags indicating */
1501 /* which windows corresponding to search regions must be */
1502 /* computed. */
1503
1504 /* Next, the windows corresponding to times when the coordinate */
1505 /* is in the selected regions are found using ZZGFREL. */
1506
1507
1508 /* Check the workspace window count. */
1509
1510 if (*nw < 15) {
1511 setmsg_("Workspace window count was # but must be at least #.", (
1512 ftnlen)52);
1513 errint_("#", nw, (ftnlen)1);
1514 errint_("#", &c__15, (ftnlen)1);
1515 sigerr_("SPICE(TOOFEWWINDOWS)", (ftnlen)20);
1516 chkout_("ZZGFLONG", (ftnlen)8);
1517 return 0;
1518 }
1519
1520 /* We can't initialize the whole workspace, but we can initialize */
1521 /* the windows we actually own. Do so. */
1522
1523 for (i__ = 1; i__ <= 7; ++i__) {
1524 ssized_(mw, &work[(i__1 = (i__ + 5) * work_dim1 - 5 - work_offset) <
1525 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
1526 i__1, "zzgflong_", (ftnlen)1287)]);
1527 }
1528
1529 /* Initialize the workspace window pool. Set up the parallel */
1530 /* array of window indices. */
1531
1532 lnkini_(&c__7, wwpool);
1533 for (i__ = 1; i__ <= 7; ++i__) {
1534 wix[(i__1 = i__ - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix", i__1,
1535 "zzgflong_", (ftnlen)1297)] = i__ + 5;
1536 }
1537
1538 /* Get an upper case, left-justified version of the */
1539 /* requested comparison operation. */
1540
1541 ljust_(relate, uop, relate_len, (ftnlen)6);
1542 ucase_(uop, uop, (ftnlen)6, (ftnlen)6);
1543
1544 /* Reject bad operators. */
1545
1546 /* Use the original operator string in the error message. */
1547
1548 i__ = isrchc_(uop, &c__7, ops, (ftnlen)6, (ftnlen)6);
1549 if (i__ == 0) {
1550 setmsg_("The comparison operator, # is not recognized. Supported qu"
1551 "antities are: <, =, >, LOCMIN, ABSMIN, LOCMAX, ABSMAX.", (
1552 ftnlen)113);
1553 errch_("#", relate, (ftnlen)1, relate_len);
1554 sigerr_("SPICE(NOTRECOGNIZED)", (ftnlen)20);
1555 chkout_("ZZGFLONG", (ftnlen)8);
1556 return 0;
1557 }
1558
1559 /* Make sure TOL is positive. */
1560
1561 if (*tol <= 0.) {
1562 setmsg_("TOL was #; must be positive.", (ftnlen)28);
1563 errdp_("#", tol, (ftnlen)1);
1564 sigerr_("SPICE(VALUEOUTOFRANGE)", (ftnlen)22);
1565 chkout_("ZZGFLONG", (ftnlen)8);
1566 return 0;
1567 }
1568
1569 /* We'll use a local tolerance equal to 1/5 of the input value. */
1570 /* This will allow us to keep the total round-off error within */
1571 /* the desired tolerance. */
1572
1573 /* Computing MAX */
1574 d__1 = 1e-7, d__2 = *tol / 10.;
1575 loctol = max(d__1,d__2);
1576
1577 /* Make sure ADJUST is non-negative. */
1578
1579 if (*adjust < 0.) {
1580 setmsg_("ADJUST was #; must be non-negative.", (ftnlen)35);
1581 errdp_("#", adjust, (ftnlen)1);
1582 sigerr_("SPICE(VALUEOUTOFRANGE)", (ftnlen)22);
1583 chkout_("ZZGFLONG", (ftnlen)8);
1584 return 0;
1585 }
1586
1587 /* Confirm ADJUST equals zero unless UOP (RELATE) has value */
1588 /* "ABSMAX" or "ABSMIN." */
1589
1590 if (s_cmp(uop, "ABSMIN", (ftnlen)6, (ftnlen)6) != 0 && s_cmp(uop, "ABSMAX"
1591 , (ftnlen)6, (ftnlen)6) != 0) {
1592 if (*adjust != 0.) {
1593 setmsg_("ADJUST should have value zero for all comparison operat"
1594 "ors except ABSMAX and ABSMIN", (ftnlen)83);
1595 sigerr_("SPICE(INVALIDVALUE)", (ftnlen)19);
1596 chkout_("ZZGFLONG", (ftnlen)8);
1597 return 0;
1598 }
1599 }
1600
1601 /* Get an upper case, left-justified, compressed versions of the */
1602 /* coordinate system and coordinate names. */
1603
1604 ljust_(crdsys, nrmsys, crdsys_len, (ftnlen)32);
1605 cmprss_(" ", &c__0, nrmsys, nrmsys, (ftnlen)1, (ftnlen)32, (ftnlen)32);
1606 ucase_(nrmsys, nrmsys, (ftnlen)32, (ftnlen)32);
1607 ljust_(crdnam, nrmcrd, crdnam_len, (ftnlen)32);
1608 cmprss_(" ", &c__1, nrmcrd, nrmcrd, (ftnlen)1, (ftnlen)32, (ftnlen)32);
1609 ucase_(nrmcrd, nrmcrd, (ftnlen)32, (ftnlen)32);
1610
1611 /* Make an initial call to the coordinate utility initialization */
1612 /* routine to invoke error checking. We don't want to have */
1613 /* to duplicate the checking here. Later, when necessary, we'll */
1614 /* re-initialize the utilities. */
1615
1616 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec, nrmsys,
1617 nrmcrd, vecdef_len, method_len, target_len, ref_len, abcorr_len,
1618 obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
1619 if (failed_()) {
1620 chkout_("ZZGFLONG", (ftnlen)8);
1621 return 0;
1622 }
1623
1624 /* We've done the basic error checking. Empty the result window and */
1625 /* return now if the confinement window is empty. */
1626
1627 if (wncard_(cnfine) == 0) {
1628 scardd_(&c__0, result);
1629 chkout_("ZZGFLONG", (ftnlen)8);
1630 return 0;
1631 }
1632
1633 /* Initialize the total number of search passes performed. */
1634
1635 total = 0;
1636
1637 /* To eliminate special cases, we'll check for inequality */
1638 /* constraints that are always met or can't be met. */
1639
1640 if (s_cmp(nrmsys, "CYLINDRICAL", (ftnlen)32, (ftnlen)11) == 0 || s_cmp(
1641 nrmsys, "PLANETOGRAPHIC", (ftnlen)32, (ftnlen)14) == 0 || s_cmp(
1642 nrmsys, "RA/DEC", (ftnlen)32, (ftnlen)6) == 0) {
1643 if (cos(*refval) == 1.) {
1644
1645 /* The reference value lies on the branch cut at 0. */
1646
1647 if (s_cmp(uop, "<", (ftnlen)6, (ftnlen)1) == 0) {
1648
1649 /* These coordinates can never be less than zero. */
1650
1651 scardd_(&c__0, result);
1652 chkout_("ZZGFLONG", (ftnlen)8);
1653 return 0;
1654 } else if (s_cmp(uop, ">", (ftnlen)6, (ftnlen)1) == 0) {
1655
1656 /* The solution is the whole confinement window. This */
1657 /* is because the inequality operators really act like */
1658 /* '>=' and '<=' operators, and because we assume the */
1659 /* quantity is increasing or decreasing except on a */
1660 /* set of measure zero. */
1661
1662 copyd_(cnfine, result);
1663 chkout_("ZZGFLONG", (ftnlen)8);
1664 return 0;
1665 }
1666 }
1667 } else if (s_cmp(nrmsys, "GEODETIC", (ftnlen)32, (ftnlen)8) == 0 || s_cmp(
1668 nrmsys, "LATITUDINAL", (ftnlen)32, (ftnlen)11) == 0 || s_cmp(
1669 nrmsys, "SPHERICAL", (ftnlen)32, (ftnlen)9) == 0) {
1670 if (cos(*refval) == -1.) {
1671
1672 /* The reference value lies on the branch cut at pi. */
1673
1674 if (s_cmp(uop, "<", (ftnlen)6, (ftnlen)1) == 0) {
1675
1676 /* The solution is the whole confinement window. */
1677
1678 copyd_(cnfine, result);
1679 chkout_("ZZGFLONG", (ftnlen)8);
1680 return 0;
1681 } else if (s_cmp(uop, ">", (ftnlen)6, (ftnlen)1) == 0) {
1682
1683 /* These coordinates can never be greater */
1684 /* than pi. */
1685
1686 scardd_(&c__0, result);
1687 chkout_("ZZGFLONG", (ftnlen)8);
1688 return 0;
1689 }
1690 }
1691 }
1692
1693 /* At this point, we make some adjustments to simplify the */
1694 /* remaining code. We map the input coordinate system to */
1695 /* either "latitudinal" or "RA/DEC" and modify the */
1696 /* constraint if the original system is "planetographic." */
1697 /* The longitude coordinate is renamed accordingly, if necessary. */
1698 /* The mapping is as follows: */
1699
1700 /* Spherical ( longitude range is (-pi, pi] ) -> Latitudinal */
1701
1702 /* Cylindrical ( longitude range is [0, 2*pi] ) -> RA/Dec */
1703 /* Longitude -> RA */
1704
1705 /* Planetographic ( longitude range is [0, 2*pi] ) -> RA/Dec */
1706 /* Longitude -> RA */
1707
1708
1709 /* For planetographic coordinates, if the longitude is positive */
1710 /* west, and since REFVAL does not lie on the branch cut, we can */
1711 /* make the following additional adjustments: */
1712
1713 /* Input relational operator Transformed operator */
1714 /* ------------------------- -------------------- */
1715 /* ABSMAX ABSMIN */
1716 /* ABSMAX - ADJUST ABSMIN + ADJUST */
1717 /* ABSMIN ABSMAX */
1718 /* ABSMIN + AJDUST ABSMAX - ADJUST */
1719 /* LOCMAX LOCMIN */
1720 /* LOCMIN LOCMAX */
1721 /* < REFVAL > 2*pi - REFVAL */
1722 /* > REFVAL < 2*pi - REFVAL */
1723 /* = REFVAL = 2*pi - REFVAL */
1724
1725
1726 xrfval = *refval;
1727 if (s_cmp(nrmsys, "SPHERICAL", (ftnlen)32, (ftnlen)9) == 0) {
1728 s_copy(nrmsys, "LATITUDINAL", (ftnlen)32, (ftnlen)11);
1729 xrfval = *refval;
1730 } else if (s_cmp(nrmsys, "CYLINDRICAL", (ftnlen)32, (ftnlen)11) == 0) {
1731 s_copy(nrmsys, "RA/DEC", (ftnlen)32, (ftnlen)6);
1732 s_copy(nrmcrd, "RIGHT ASCENSION", (ftnlen)32, (ftnlen)15);
1733 xrfval = *refval;
1734 } else if (s_cmp(nrmsys, "PLANETOGRAPHIC", (ftnlen)32, (ftnlen)14) == 0) {
1735 s_copy(nrmsys, "RA/DEC", (ftnlen)32, (ftnlen)6);
1736 s_copy(nrmcrd, "RIGHT ASCENSION", (ftnlen)32, (ftnlen)15);
1737
1738 /* If the planetographic coordinates are positive West, we'll */
1739 /* need to transform the constraint and reference value. */
1740
1741 /* Get the name of the central body of frame REF. */
1742
1743 /* NOTE: We omit error checking here because ZZGFCOIN has done */
1744 /* it already. */
1745
1746 namfrm_(ref, &frcode, ref_len);
1747 frinfo_(&frcode, &refctr, &class__, &clssid, &found);
1748 if (failed_()) {
1749 chkout_("ZZGFLONG", (ftnlen)8);
1750 return 0;
1751 }
1752 if (! found) {
1753 setmsg_("FRINFO didn't find data for frame # which has frame ID "
1754 "code #. This frame should have been validated by ZZGFCOI"
1755 "N.", (ftnlen)113);
1756 errch_("#", ref, (ftnlen)1, ref_len);
1757 errint_("#", &frcode, (ftnlen)1);
1758 sigerr_("SPICE(BUG)", (ftnlen)10);
1759 chkout_("ZZGFLONG", (ftnlen)8);
1760 return 0;
1761 }
1762 bodc2s_(&refctr, rctrnm, (ftnlen)36);
1763
1764 /* Find the longitude of the +Y axis. If this longitude */
1765 /* is greater than pi, the sense is positive West. Note */
1766 /* that we don't need to use realistic values of the */
1767 /* equatorial radius and flattening factor: 1 and 0, */
1768 /* respectively, are just fine. */
1769
1770 recpgr_(rctrnm, y, &c_b69, &c_b70, &lon, &lat, &alt, (ftnlen)36);
1771
1772 /* Planetographic longitude ranges from 0 to 2*pi, so */
1773 /* longitudes corresponding to positive Y values are */
1774 /* in the range pi to 2*pi. */
1775
1776 if (lon > pi_()) {
1777
1778 /* Planetographic longitude for the frame center is positive */
1779 /* West. */
1780
1781 /* Note that no action is required to modify non-zero */
1782 /* extremum adjustment values. */
1783
1784 if (s_cmp(uop, "ABSMAX", (ftnlen)6, (ftnlen)6) == 0) {
1785 s_copy(uop, "ABSMIN", (ftnlen)6, (ftnlen)6);
1786 } else if (s_cmp(uop, "ABSMIN", (ftnlen)6, (ftnlen)6) == 0) {
1787 s_copy(uop, "ABSMAX", (ftnlen)6, (ftnlen)6);
1788 } else if (s_cmp(uop, "LOCMAX", (ftnlen)6, (ftnlen)6) == 0) {
1789 s_copy(uop, "LOCMIN", (ftnlen)6, (ftnlen)6);
1790 } else if (s_cmp(uop, "LOCMIN", (ftnlen)6, (ftnlen)6) == 0) {
1791 s_copy(uop, "LOCMAX", (ftnlen)6, (ftnlen)6);
1792 } else if (s_cmp(uop, "=", (ftnlen)6, (ftnlen)1) == 0) {
1793 xrfval = twopi_() - *refval;
1794 } else if (s_cmp(uop, "<", (ftnlen)6, (ftnlen)1) == 0) {
1795 s_copy(uop, ">", (ftnlen)6, (ftnlen)1);
1796 xrfval = twopi_() - *refval;
1797 } else if (s_cmp(uop, ">", (ftnlen)6, (ftnlen)1) == 0) {
1798 s_copy(uop, "<", (ftnlen)6, (ftnlen)1);
1799 xrfval = twopi_() - *refval;
1800 } else {
1801
1802 /* We shouldn't get here. */
1803
1804 setmsg_("Unexpected UOP value: #", (ftnlen)23);
1805 errch_("#", uop, (ftnlen)1, (ftnlen)6);
1806 sigerr_("SPICE(BUG)", (ftnlen)10);
1807 chkout_("ZZGFLONG", (ftnlen)8);
1808 return 0;
1809 }
1810 } else {
1811
1812 /* Longitude is positive East, so we treat */
1813 /* the constraint as though the coordinate were RA. */
1814
1815 xrfval = *refval;
1816 }
1817 }
1818
1819 /* From this point on, we use: */
1820
1821 /* Coordinate system: NRMSYS */
1822 /* Coordinate: NRMCRD */
1823 /* Operator: UOP */
1824 /* Reference value: XRFVAL */
1825
1826
1827 /* The result window must be initialized by the caller of the GF */
1828 /* system (usually a user application). We simply empty the result */
1829 /* window here. */
1830
1831 scardd_(&c__0, result);
1832
1833 /* We use the constant 0.5 * 2**0.5 quite a bit. Create a */
1834 /* "macro" variable for it. */
1835
1836 r2ovr2 = sqrt(2.) / 2.;
1837
1838 /* Set the progress report suffix strings. */
1839
1840 s_copy(rptsuf, "done.", (ftnlen)80, (ftnlen)5);
1841 s_copy(rptsuf + 80, "done.", (ftnlen)80, (ftnlen)5);
1842
1843 /* Case: '=' */
1844
1845 if (s_cmp(uop, "=", (ftnlen)6, (ftnlen)1) == 0) {
1846
1847 /* Equality constraints are the simplest to handle, so we'll get */
1848 /* them out of the way now. Our approach is to use sine or cosine */
1849 /* as proxy functions; we'll select the proxy function with the */
1850 /* highest resolution at the reference value. For the proxy */
1851 /* function f, our proxy constraint is */
1852
1853 /* f(x) = f(XRFVAL) */
1854
1855 /* This may yield spurious roots; we'll delete these after we've */
1856 /* done our search. */
1857
1858 /* Find the sine and cosine of the reference value. We'll use */
1859 /* these both to locate the quadrant of the reference value and */
1860 /* to have continuously differentiable functions to work with. */
1861 /* Note that if the original reference value is not in the */
1862 /* standard range, this presents no problem. */
1863
1864 cv = cos(xrfval);
1865 sv = sin(xrfval);
1866
1867 /* Decide which proxy function to use. */
1868
1869 if (abs(sv) >= r2ovr2) {
1870
1871 /* The reference value lies in the top or bottom quarter of */
1872 /* the unit circle. The "comparison value" CMPVAL will be */
1873 /* used later to delete solutions with matching sines but */
1874 /* non-matching cosines. */
1875
1876 s_copy(prxfun, "COS", (ftnlen)50, (ftnlen)3);
1877 prxval = cv;
1878 cmpval = sv;
1879 } else {
1880 s_copy(prxfun, "SIN", (ftnlen)50, (ftnlen)3);
1881 prxval = sv;
1882 cmpval = cv;
1883 }
1884
1885 /* Set up the progress reporting prefix strings. We have one */
1886 /* ZZGFREL call which performs two passes. */
1887
1888 s_copy(rptpre, "Coordinate pass 1 of 2", (ftnlen)80, (ftnlen)22);
1889 s_copy(rptpre + 80, "Coordinate pass 2 of 2", (ftnlen)80, (ftnlen)22);
1890
1891 /* Allocate a workspace window. */
1892
1893 lnkan_(wwpool, &node);
1894 f1 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
1895 i__1, "zzgflong_", (ftnlen)1749)];
1896
1897 /* Make sure the coordinate utilities have been initialized */
1898 /* with the actual values we'll use for our search. */
1899
1900 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
1901 nrmsys, nrmcrd, vecdef_len, method_len, target_len, ref_len,
1902 abcorr_len, obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
1903
1904 /* Now we're ready to compute the window in which our proxy */
1905 /* function satisfies the proxy constraint. */
1906
1907 if (s_cmp(prxfun, "SIN", (ftnlen)50, (ftnlen)3) == 0) {
1908
1909 /* Find the window where the sine of the coordinate satisfies */
1910 /* the proxy constraint. */
1911
1912 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcosd_, (U_fp)
1913 zzgfudlt_, (S_fp)zzgfcosg_, "=", &prxval, &loctol, &c_b70,
1914 cnfine, mw, nw, work, rpt, (U_fp)udrepi, (U_fp)udrepu, (
1915 U_fp)udrepf, rptpre, rptsuf, bail, (L_fp)udbail, &work[(
1916 i__1 = f1 * work_dim1 - 5 - work_offset) < work_dim1 *
1917 work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work", i__1,
1918 "zzgflong_", (ftnlen)1769)], (ftnlen)1, (ftnlen)80, (
1919 ftnlen)80);
1920 } else {
1921
1922 /* Find the window where the cosine of the coordinate */
1923 /* satisfies the proxy constraint. */
1924
1925 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcocd_, (U_fp)
1926 zzgfudlt_, (S_fp)zzgfcocg_, "=", &prxval, &loctol, &c_b70,
1927 cnfine, mw, nw, work, rpt, (U_fp)udrepi, (U_fp)udrepu, (
1928 U_fp)udrepf, rptpre, rptsuf, bail, (L_fp)udbail, &work[(
1929 i__1 = f1 * work_dim1 - 5 - work_offset) < work_dim1 *
1930 work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work", i__1,
1931 "zzgflong_", (ftnlen)1784)], (ftnlen)1, (ftnlen)80, (
1932 ftnlen)80);
1933 }
1934 if (failed_()) {
1935 chkout_("ZZGFLONG", (ftnlen)8);
1936 return 0;
1937 }
1938
1939 /* Handle interrupts if necessary. */
1940
1941 if (*bail) {
1942 if ((*udbail)()) {
1943 chkout_("ZZGFLONG", (ftnlen)8);
1944 return 0;
1945 }
1946 }
1947
1948 /* Remove any spurious results. */
1949
1950 n = cardd_(&work[(i__1 = f1 * work_dim1 - 5 - work_offset) <
1951 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
1952 i__1, "zzgflong_", (ftnlen)1813)]);
1953 i__1 = n;
1954 for (i__ = 1; i__ <= i__1; i__ += 2) {
1955 start = work[(i__2 = i__ + f1 * work_dim1 - work_offset) <
1956 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge("work",
1957 i__2, "zzgflong_", (ftnlen)1817)];
1958 if (s_cmp(prxfun, "SIN", (ftnlen)50, (ftnlen)3) == 0) {
1959
1960 /* Get the cosine of the coordinate at the interval start */
1961 /* time. If this cosine has the same sign as the cosine of */
1962 /* the reference value, we have a winner. Note that the */
1963 /* cosines of spurious values won't ever be close to the */
1964 /* correct values, so round-off isn't an issue. */
1965
1966 zzgfcocg_(&start, &value);
1967 } else {
1968
1969 /* Same deal, but here we're using sines. */
1970
1971 zzgfcosg_(&start, &value);
1972 }
1973 if (smsgnd_(&cmpval, &value)) {
1974
1975 /* This is a winner. */
1976
1977 wninsd_(&start, &start, result);
1978 }
1979 }
1980
1981 /* All done. */
1982
1983 chkout_("ZZGFLONG", (ftnlen)8);
1984 return 0;
1985 }
1986
1987 /* Case: local minimum or local maximum */
1988
1989 if (s_cmp(uop, "LOCMAX", (ftnlen)6, (ftnlen)6) == 0 || s_cmp(uop, "LOCMIN"
1990 , (ftnlen)6, (ftnlen)6) == 0) {
1991
1992 /* This algorithm uses 4 ZZGFREL calls, 2 of which perform */
1993 /* 2 passes and 2 of which perform 1 pass. */
1994
1995 s_copy(rptsuf, "done.", (ftnlen)80, (ftnlen)5);
1996 s_copy(rptsuf + 80, "done.", (ftnlen)80, (ftnlen)5);
1997
1998 /* Empty the result window. */
1999
2000 scardd_(&c__0, result);
2001
2002 /* We'll first find two windows covering the left and right */
2003 /* halves of the unit circle, with both halves extended */
2004 /* slightly to ensure no roots are missed. We start by */
2005 /* finding the window on which the cosine of the coordinate */
2006 /* is less than cos(LCXMRG) (which is a small, positive number). */
2007
2008 lnkan_(wwpool, &node);
2009 left = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2010 i__1, "zzgflong_", (ftnlen)1880)];
2011 s_copy(rptpre, "Coordinate pass 1 of 6", (ftnlen)80, (ftnlen)22);
2012 s_copy(rptpre + 80, "Coordinate pass 2 of 6", (ftnlen)80, (ftnlen)22);
2013 s_copy(prxrel, "<", (ftnlen)6, (ftnlen)1);
2014 prxval = cos(1e-12);
2015 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
2016 nrmsys, nrmcrd, vecdef_len, method_len, target_len, ref_len,
2017 abcorr_len, obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
2018 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcocd_, (U_fp)
2019 zzgfudlt_, (S_fp)zzgfcocg_, prxrel, &prxval, &loctol, &c_b70,
2020 cnfine, mw, nw, work, rpt, (U_fp)udrepi, (U_fp)udrepu, (U_fp)
2021 udrepf, rptpre, rptsuf, bail, (L_fp)udbail, &work[(i__1 =
2022 left * work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 &&
2023 0 <= i__1 ? i__1 : s_rnge("work", i__1, "zzgflong_", (ftnlen)
2024 1892)], (ftnlen)6, (ftnlen)80, (ftnlen)80);
2025
2026 /* Handle interrupts if necessary. */
2027
2028 if (*bail) {
2029 if ((*udbail)()) {
2030 chkout_("ZZGFLONG", (ftnlen)8);
2031 return 0;
2032 }
2033 }
2034
2035 /* Now search for the time period when the cosine of the */
2036 /* coordinate is greater than -cos(LCXMRG). We can save some time */
2037 /* by searching within the window designated by LEFT for the */
2038 /* complement of this window and then complementing the result of */
2039 /* that search. */
2040
2041 lnkan_(wwpool, &node);
2042 compl = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2043 i__1, "zzgflong_", (ftnlen)1919)];
2044 lnkan_(wwpool, &node);
2045 right = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2046 i__1, "zzgflong_", (ftnlen)1922)];
2047 s_copy(rptpre, "Coordinate pass 3 of 6", (ftnlen)80, (ftnlen)22);
2048 s_copy(rptpre + 80, "Coordinate pass 4 of 6", (ftnlen)80, (ftnlen)22);
2049 s_copy(prxrel, "<", (ftnlen)6, (ftnlen)1);
2050 prxval = -cos(1e-12);
2051 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
2052 nrmsys, nrmcrd, vecdef_len, method_len, target_len, ref_len,
2053 abcorr_len, obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
2054 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcocd_, (U_fp)
2055 zzgfudlt_, (S_fp)zzgfcocg_, prxrel, &prxval, &loctol, &c_b70,
2056 &work[(i__1 = left * work_dim1 - 5 - work_offset) < work_dim1
2057 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work", i__1, "zzgf"
2058 "long_", (ftnlen)1935)], mw, nw, work, rpt, (U_fp)udrepi, (
2059 U_fp)udrepu, (U_fp)udrepf, rptpre, rptsuf, bail, (L_fp)udbail,
2060 &work[(i__2 = compl * work_dim1 - 5 - work_offset) <
2061 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge("work",
2062 i__2, "zzgflong_", (ftnlen)1935)], (ftnlen)6, (ftnlen)80, (
2063 ftnlen)80);
2064
2065 /* Handle interrupts if necessary. */
2066
2067 if (*bail) {
2068 if ((*udbail)()) {
2069 chkout_("ZZGFLONG", (ftnlen)8);
2070 return 0;
2071 }
2072 }
2073
2074 /* WORK(LB,COMPL) contains the complement of the window */
2075 /* we want. */
2076
2077 wndifd_(cnfine, &work[(i__1 = compl * work_dim1 - 5 - work_offset) <
2078 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2079 i__1, "zzgflong_", (ftnlen)1958)], &work[(i__2 = right *
2080 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 && 0 <=
2081 i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (ftnlen)1958)]
2082 );
2083
2084 /* We're now going to find local extrema of the coordinate in the */
2085 /* windows indexed by LEFT and RIGHT. */
2086
2087 for (i__ = 1; i__ <= 2; ++i__) {
2088 if (i__ == 1) {
2089
2090 /* The sector we're searching is indexed by LEFT. */
2091 /* We'll use RA as a proxy function, since RA has no */
2092 /* singularity on the left half circle. */
2093
2094 s = left;
2095 s_copy(prxsys, "RA/DEC", (ftnlen)32, (ftnlen)6);
2096 s_copy(prxcrd, "RIGHT ASCENSION", (ftnlen)32, (ftnlen)15);
2097 lnkan_(wwpool, &node);
2098 res1 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge(
2099 "wix", i__1, "zzgflong_", (ftnlen)1977)];
2100 res = res1;
2101 s_copy(rptpre, "Coordinate pass 5 of 6", (ftnlen)80, (ftnlen)
2102 22);
2103 s_copy(rptpre + 80, " ", (ftnlen)80, (ftnlen)1);
2104 } else {
2105 s = right;
2106 s_copy(prxsys, "LATITUDINAL", (ftnlen)32, (ftnlen)11);
2107 s_copy(prxcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9);
2108 lnkan_(wwpool, &node);
2109 res2 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge(
2110 "wix", i__1, "zzgflong_", (ftnlen)1990)];
2111 res = res2;
2112 s_copy(rptpre, "Coordinate pass 6 of 6", (ftnlen)80, (ftnlen)
2113 22);
2114 s_copy(rptpre + 80, " ", (ftnlen)80, (ftnlen)1);
2115 }
2116 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
2117 prxsys, prxcrd, vecdef_len, method_len, target_len,
2118 ref_len, abcorr_len, obsrvr_len, dref_len, (ftnlen)32, (
2119 ftnlen)32);
2120 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcodc_, (U_fp)
2121 zzgfudlt_, (S_fp)zzgfcog_, uop, &c_b70, &loctol, &c_b70, &
2122 work[(i__1 = s * work_dim1 - 5 - work_offset) < work_dim1
2123 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work", i__1,
2124 "zzgflong_", (ftnlen)2003)], mw, nw, work, rpt, (U_fp)
2125 udrepi, (U_fp)udrepu, (U_fp)udrepf, rptpre, rptsuf, bail,
2126 (L_fp)udbail, &work[(i__2 = res * work_dim1 - 5 -
2127 work_offset) < work_dim1 * work_dim2 && 0 <= i__2 ? i__2 :
2128 s_rnge("work", i__2, "zzgflong_", (ftnlen)2003)], (
2129 ftnlen)6, (ftnlen)80, (ftnlen)80);
2130
2131 /* Handle interrupts if necessary. */
2132
2133 if (*bail) {
2134 if ((*udbail)()) {
2135 chkout_("ZZGFLONG", (ftnlen)8);
2136 return 0;
2137 }
2138 }
2139 }
2140
2141 /* Combine the contributions of both searches in RESULT. */
2142
2143 wnunid_(&work[(i__1 = res1 * work_dim1 - 5 - work_offset) < work_dim1
2144 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work", i__1, "zzgf"
2145 "long_", (ftnlen)2027)], &work[(i__2 = res2 * work_dim1 - 5 -
2146 work_offset) < work_dim1 * work_dim2 && 0 <= i__2 ? i__2 :
2147 s_rnge("work", i__2, "zzgflong_", (ftnlen)2027)], result);
2148
2149 /* End of the LOCMIN and LOCMAX cases. RESULT is set. */
2150
2151 chkout_("ZZGFLONG", (ftnlen)8);
2152 return 0;
2153 }
2154
2155 /* The remaining operators are: ABSMAX, ABSMIN, '<', '>'. */
2156
2157 /* Initialize the window aliases. A value of zero indicates the */
2158 /* corresponding region hasn't been computed. */
2159
2160 top = 0;
2161 bot = 0;
2162 right = 0;
2163 left = 0;
2164 q1 = 0;
2165 q2 = 0;
2166 q3 = 0;
2167 q4 = 0;
2168 s = 0;
2169 wh = 0;
2170 f1 = 0;
2171 f2 = 0;
2172
2173 /* If we have an absolute extremum or inequality relation, */
2174 /* we'll need to find times when the coordinate is in the */
2175 /* various quadrants. We'll start out by setting up windows */
2176 /* for the times when the coordinate is in the top and right */
2177 /* halves of the unit circle. */
2178
2179 /* The ZZGFREL call below involves two passes. */
2180
2181 if (s_cmp(uop, "ABSMAX", (ftnlen)6, (ftnlen)6) == 0 || s_cmp(uop, "ABSMIN"
2182 , (ftnlen)6, (ftnlen)6) == 0) {
2183 if (*adjust == 0.) {
2184 s_copy(tmplat, "Coordinate pass # of 7", (ftnlen)80, (ftnlen)22);
2185 } else {
2186 s_copy(tmplat, "Coordinate pass # of 7-9", (ftnlen)80, (ftnlen)24)
2187 ;
2188 }
2189 } else {
2190
2191 /* Ordinary inequality searches use 8 passes. */
2192
2193 s_copy(tmplat, "Coordinate pass # of 8", (ftnlen)80, (ftnlen)22);
2194 }
2195 for (i__ = 1; i__ <= 2; ++i__) {
2196 repmi_(tmplat, "#", &i__, rptpre + ((i__1 = i__ - 1) < 2 && 0 <= i__1
2197 ? i__1 : s_rnge("rptpre", i__1, "zzgflong_", (ftnlen)2083)) *
2198 80, (ftnlen)80, (ftnlen)1, (ftnlen)80);
2199 }
2200
2201 /* Find the window where the sine of the coordinate is greater than */
2202 /* the sine of the branch cut avoidance tolerance. */
2203
2204 /* Make sure the coordinate utilities have been initialized */
2205 /* with the actual values we'll use for our search. */
2206
2207 lnkan_(wwpool, &node);
2208 head = node;
2209 top = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix", i__1,
2210 "zzgflong_", (ftnlen)2095)];
2211 prxval = sin(1e-11);
2212 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec, nrmsys,
2213 nrmcrd, vecdef_len, method_len, target_len, ref_len, abcorr_len,
2214 obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
2215 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcosd_, (U_fp)zzgfudlt_, (
2216 S_fp)zzgfcosg_, ">", &prxval, &loctol, &c_b70, cnfine, mw, nw,
2217 work, rpt, (U_fp)udrepi, (U_fp)udrepu, (U_fp)udrepf, rptpre,
2218 rptsuf, bail, (L_fp)udbail, &work[(i__1 = top * work_dim1 - 5 -
2219 work_offset) < work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2220 "work", i__1, "zzgflong_", (ftnlen)2103)], (ftnlen)1, (ftnlen)80,
2221 (ftnlen)80);
2222
2223 /* 2 passes done. */
2224
2225 total = 2;
2226 if (*bail) {
2227 if ((*udbail)()) {
2228 chkout_("ZZGFLONG", (ftnlen)8);
2229 return 0;
2230 }
2231 }
2232
2233 /* Find the window where the sine of the coordinate is less than */
2234 /* the negative of the sine of the branch cut avoidance tolerance. */
2235
2236 /* Make sure the coordinate utilities have been initialized */
2237 /* with the actual values we'll use for our search. */
2238
2239 /* The ZZGFREL call below involves two passes. */
2240
2241 for (i__ = 1; i__ <= 2; ++i__) {
2242 i__2 = total + i__;
2243 repmi_(tmplat, "#", &i__2, rptpre + ((i__1 = i__ - 1) < 2 && 0 <=
2244 i__1 ? i__1 : s_rnge("rptpre", i__1, "zzgflong_", (ftnlen)
2245 2134)) * 80, (ftnlen)80, (ftnlen)1, (ftnlen)80);
2246 }
2247 lnkan_(wwpool, &node);
2248 lnkila_(&head, &node, wwpool);
2249 bot = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix", i__1,
2250 "zzgflong_", (ftnlen)2140)];
2251 prxval = -sin(1e-11);
2252 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec, nrmsys,
2253 nrmcrd, vecdef_len, method_len, target_len, ref_len, abcorr_len,
2254 obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
2255 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcosd_, (U_fp)zzgfudlt_, (
2256 S_fp)zzgfcosg_, "<", &prxval, &loctol, &c_b70, cnfine, mw, nw,
2257 work, rpt, (U_fp)udrepi, (U_fp)udrepu, (U_fp)udrepf, rptpre,
2258 rptsuf, bail, (L_fp)udbail, &work[(i__1 = bot * work_dim1 - 5 -
2259 work_offset) < work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2260 "work", i__1, "zzgflong_", (ftnlen)2149)], (ftnlen)1, (ftnlen)80,
2261 (ftnlen)80);
2262
2263 /* 4 passes done. */
2264
2265 total += 2;
2266 if (*bail) {
2267 if ((*udbail)()) {
2268 chkout_("ZZGFLONG", (ftnlen)8);
2269 return 0;
2270 }
2271 }
2272
2273 /* Find the window where the cosine of the coordinate is */
2274 /* greater than zero. */
2275
2276
2277 /* The ZZGFREL call below involves two passes. */
2278
2279 for (i__ = 1; i__ <= 2; ++i__) {
2280 i__2 = total + i__;
2281 repmi_(tmplat, "#", &i__2, rptpre + ((i__1 = i__ - 1) < 2 && 0 <=
2282 i__1 ? i__1 : s_rnge("rptpre", i__1, "zzgflong_", (ftnlen)
2283 2178)) * 80, (ftnlen)80, (ftnlen)1, (ftnlen)80);
2284 }
2285
2286 /* We'll keep all of the allocated nodes linked together. */
2287 /* Since the order of the nodes is unimportant, we insert */
2288 /* each new node following the head node; this is non-standard */
2289 /* but ensures the list head doesn't change until we delete */
2290 /* nodes from the list. */
2291
2292 lnkan_(wwpool, &node);
2293 lnkila_(&head, &node, wwpool);
2294 right = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2295 i__1, "zzgflong_", (ftnlen)2190)];
2296 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec, nrmsys,
2297 nrmcrd, vecdef_len, method_len, target_len, ref_len, abcorr_len,
2298 obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
2299 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcocd_, (U_fp)zzgfudlt_, (
2300 S_fp)zzgfcocg_, ">", &c_b70, &loctol, &c_b70, cnfine, mw, nw,
2301 work, rpt, (U_fp)udrepi, (U_fp)udrepu, (U_fp)udrepf, rptpre,
2302 rptsuf, bail, (L_fp)udbail, &work[(i__1 = right * work_dim1 - 5 -
2303 work_offset) < work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2304 "work", i__1, "zzgflong_", (ftnlen)2196)], (ftnlen)1, (ftnlen)80,
2305 (ftnlen)80);
2306
2307 /* 6 passes done. */
2308
2309 total += 2;
2310 if (*bail) {
2311 if ((*udbail)()) {
2312 chkout_("ZZGFLONG", (ftnlen)8);
2313 return 0;
2314 }
2315 }
2316 if (failed_()) {
2317 chkout_("ZZGFLONG", (ftnlen)8);
2318 return 0;
2319 }
2320
2321 /* Now find the absolute extremum, if this was requested. */
2322
2323 if (s_cmp(uop, "ABSMAX", (ftnlen)6, (ftnlen)6) == 0 || s_cmp(uop, "ABSMIN"
2324 , (ftnlen)6, (ftnlen)6) == 0) {
2325
2326 /* If we're looking for an absolute extremum and the */
2327 /* adjustment value is 0, each ZZGFREL call executes */
2328 /* on search pass; otherwise these calls execute two */
2329 /* search passes. */
2330
2331 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2332
2333 /* We need windows when the coordinate is in quadrants 2 and */
2334 /* 3. We can derive these from the windows TOP and RIGHT */
2335 /* without additional searches. */
2336
2337 lnkan_(wwpool, &node);
2338 lnkila_(&head, &node, wwpool);
2339 q2 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2340 i__1, "zzgflong_", (ftnlen)2242)];
2341 lnkan_(wwpool, &node);
2342 lnkila_(&head, &node, wwpool);
2343 q3 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2344 i__1, "zzgflong_", (ftnlen)2246)];
2345
2346 /* Compute windows for the second and third quadrants. Note */
2347 /* that these windows are bounded away from the branch cut */
2348 /* at pi radians, since windows TOP and BOT have been */
2349 /* trimmed. */
2350
2351 wndifd_(&work[(i__1 = top * work_dim1 - 5 - work_offset) <
2352 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2353 i__1, "zzgflong_", (ftnlen)2254)], &work[(i__2 = right *
2354 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 && 0
2355 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
2356 ftnlen)2254)], &work[(i__3 = q2 * work_dim1 - 5 -
2357 work_offset) < work_dim1 * work_dim2 && 0 <= i__3 ? i__3 :
2358 s_rnge("work", i__3, "zzgflong_", (ftnlen)2254)]);
2359 wndifd_(&work[(i__1 = bot * work_dim1 - 5 - work_offset) <
2360 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2361 i__1, "zzgflong_", (ftnlen)2255)], &work[(i__2 = right *
2362 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 && 0
2363 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
2364 ftnlen)2255)], &work[(i__3 = q3 * work_dim1 - 5 -
2365 work_offset) < work_dim1 * work_dim2 && 0 <= i__3 ? i__3 :
2366 s_rnge("work", i__3, "zzgflong_", (ftnlen)2255)]);
2367 if (s_cmp(uop, "ABSMAX", (ftnlen)6, (ftnlen)6) == 0) {
2368 region[0] = q2;
2369 region[1] = right;
2370 region[2] = q3;
2371 } else {
2372 region[0] = q3;
2373 region[1] = right;
2374 region[2] = q2;
2375 }
2376 } else if (s_cmp(nrmcrd, "RIGHT ASCENSION", (ftnlen)32, (ftnlen)15) ==
2377 0) {
2378
2379 /* We need windows when the coordinate is in quadrants 1 and */
2380 /* 4, and the window when the coordinate is in the left half */
2381 /* of the unit circle. We can derive these from the windows */
2382 /* TOP and RIGHT without additional searches. */
2383
2384 lnkan_(wwpool, &node);
2385 lnkila_(&head, &node, wwpool);
2386 q1 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2387 i__1, "zzgflong_", (ftnlen)2278)];
2388 lnkan_(wwpool, &node);
2389 lnkila_(&head, &node, wwpool);
2390 left = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge(
2391 "wix", i__1, "zzgflong_", (ftnlen)2282)];
2392 lnkan_(wwpool, &node);
2393 lnkila_(&head, &node, wwpool);
2394 q4 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2395 i__1, "zzgflong_", (ftnlen)2286)];
2396
2397 /* Compute windows for the first and fourth quadrants. Note */
2398 /* that these windows are bounded away from the branch cut */
2399 /* at pi radians, since windows TOP and BOT have been */
2400 /* trimmed. Also compute the window LEFT, which is the */
2401 /* complement of window RIGHT. */
2402
2403 wnintd_(&work[(i__1 = right * work_dim1 - 5 - work_offset) <
2404 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2405 i__1, "zzgflong_", (ftnlen)2295)], &work[(i__2 = top *
2406 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 && 0
2407 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
2408 ftnlen)2295)], &work[(i__3 = q1 * work_dim1 - 5 -
2409 work_offset) < work_dim1 * work_dim2 && 0 <= i__3 ? i__3 :
2410 s_rnge("work", i__3, "zzgflong_", (ftnlen)2295)]);
2411 wnintd_(&work[(i__1 = right * work_dim1 - 5 - work_offset) <
2412 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2413 i__1, "zzgflong_", (ftnlen)2296)], &work[(i__2 = bot *
2414 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 && 0
2415 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
2416 ftnlen)2296)], &work[(i__3 = q4 * work_dim1 - 5 -
2417 work_offset) < work_dim1 * work_dim2 && 0 <= i__3 ? i__3 :
2418 s_rnge("work", i__3, "zzgflong_", (ftnlen)2296)]);
2419 wndifd_(cnfine, &work[(i__1 = right * work_dim1 - 5 - work_offset)
2420 < work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2421 "work", i__1, "zzgflong_", (ftnlen)2297)], &work[(i__2 =
2422 left * work_dim1 - 5 - work_offset) < work_dim1 *
2423 work_dim2 && 0 <= i__2 ? i__2 : s_rnge("work", i__2,
2424 "zzgflong_", (ftnlen)2297)]);
2425 if (s_cmp(uop, "ABSMAX", (ftnlen)6, (ftnlen)6) == 0) {
2426 region[0] = q4;
2427 region[1] = left;
2428 region[2] = q1;
2429 } else {
2430 region[0] = q1;
2431 region[1] = left;
2432 region[2] = q4;
2433 }
2434 } else {
2435
2436 /* We're not expecting to see a coordinate other than */
2437 /* longitude or RA here. */
2438
2439 setmsg_("Unexpected coordinate # (0)", (ftnlen)27);
2440 errch_("#", nrmcrd, (ftnlen)1, (ftnlen)32);
2441 sigerr_("SPICE(BUG)", (ftnlen)10);
2442 chkout_("ZZGFLONG", (ftnlen)8);
2443 return 0;
2444 }
2445
2446 /* Now search the list of regions for the specified */
2447 /* extremum. */
2448
2449 found = FALSE_;
2450 i__ = 1;
2451 while(i__ <= 3 && ! found) {
2452
2453 /* Search region I. Set the reference and adjustment */
2454 /* values to 0 for this search. */
2455
2456 /* The ZZGFREL call below executes 1 pass, since it's */
2457 /* doing an absolute extremum search with 0 adjustment */
2458 /* value (even if ADJUST is non-zero). */
2459
2460 i__1 = total + 1;
2461 repmi_(tmplat, "#", &i__1, rptpre, (ftnlen)80, (ftnlen)1, (ftnlen)
2462 80);
2463 s_copy(rptpre + 80, " ", (ftnlen)80, (ftnlen)1);
2464 scardd_(&c__0, result);
2465
2466 /* Perform our searches with functions that have no branch */
2467 /* cuts near the region boundaries. */
2468
2469 if (region[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
2470 "region", i__1, "zzgflong_", (ftnlen)2347)] == q1 ||
2471 region[(i__2 = i__ - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge(
2472 "region", i__2, "zzgflong_", (ftnlen)2347)] == q4 ||
2473 region[(i__3 = i__ - 1) < 3 && 0 <= i__3 ? i__3 : s_rnge(
2474 "region", i__3, "zzgflong_", (ftnlen)2347)] == right) {
2475 s_copy(prxsys, "LATITUDINAL", (ftnlen)32, (ftnlen)11);
2476 s_copy(prxcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9);
2477 } else {
2478 s_copy(prxsys, "RA/DEC", (ftnlen)32, (ftnlen)6);
2479 s_copy(prxcrd, "RIGHT ASCENSION", (ftnlen)32, (ftnlen)15);
2480 }
2481 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
2482 prxsys, prxcrd, vecdef_len, method_len, target_len,
2483 ref_len, abcorr_len, obsrvr_len, dref_len, (ftnlen)32, (
2484 ftnlen)32);
2485 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcodc_, (U_fp)
2486 zzgfudlt_, (S_fp)zzgfcocg_, uop, &c_b70, &loctol, &c_b70,
2487 &work[(i__2 = region[(i__1 = i__ - 1) < 3 && 0 <= i__1 ?
2488 i__1 : s_rnge("region", i__1, "zzgflong_", (ftnlen)2362)]
2489 * work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 &&
2490 0 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
2491 ftnlen)2362)], mw, nw, work, rpt, (U_fp)udrepi, (U_fp)
2492 udrepu, (U_fp)udrepf, rptpre, rptsuf, bail, (L_fp)udbail,
2493 result, (ftnlen)6, (ftnlen)80, (ftnlen)80);
2494
2495 /* ZZGFREL will have performed a pass only if the confinement */
2496 /* window was non-empty. */
2497
2498 if (cardd_(&work[(i__2 = region[(i__1 = i__ - 1) < 3 && 0 <= i__1
2499 ? i__1 : s_rnge("region", i__1, "zzgflong_", (ftnlen)2375)
2500 ] * work_dim1 - 5 - work_offset) < work_dim1 * work_dim2
2501 && 0 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
2502 ftnlen)2375)]) > 0) {
2503
2504 /* Another pass has been completed. */
2505
2506 ++total;
2507 }
2508 if (*bail) {
2509 if ((*udbail)()) {
2510 chkout_("ZZGFLONG", (ftnlen)8);
2511 return 0;
2512 }
2513 }
2514 if (wncard_(result) > 0) {
2515
2516 /* We found an extremum. We don't have to search further. */
2517
2518 found = TRUE_;
2519 } else {
2520 ++i__;
2521 }
2522 }
2523 if (*adjust == 0.) {
2524
2525 /* The result we have is the final result. */
2526
2527 chkout_("ZZGFLONG", (ftnlen)8);
2528 return 0;
2529 }
2530
2531 /* This is the case of an absolute extremum search with */
2532 /* non-zero adjustment value. */
2533
2534 /* We'll need to obtain the extreme value. */
2535
2536 et = result[6];
2537 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
2538 nrmsys, nrmcrd, vecdef_len, method_len, target_len, ref_len,
2539 abcorr_len, obsrvr_len, dref_len, (ftnlen)32, (ftnlen)32);
2540 zzgfcog_(&et, &extval);
2541
2542 /* Re-set the operator and reference value to enable */
2543 /* us to conduct an inequality search. */
2544
2545 if (s_cmp(uop, "ABSMAX", (ftnlen)6, (ftnlen)6) == 0) {
2546 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2547 /* Computing MAX */
2548 d__1 = extval - *adjust, d__2 = -pi_();
2549 xrfval = max(d__1,d__2);
2550 } else {
2551 /* Computing MAX */
2552 d__1 = extval - *adjust;
2553 xrfval = max(d__1,0.);
2554 }
2555 s_copy(uop, ">", (ftnlen)6, (ftnlen)1);
2556 } else {
2557 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2558 /* Computing MIN */
2559 d__1 = extval + *adjust, d__2 = pi_();
2560 xrfval = min(d__1,d__2);
2561 } else {
2562 /* Computing MIN */
2563 d__1 = extval + *adjust, d__2 = twopi_();
2564 xrfval = min(d__1,d__2);
2565 }
2566 s_copy(uop, "<", (ftnlen)6, (ftnlen)1);
2567 }
2568 }
2569
2570 /* Case: inequality */
2571
2572 /* Searches for absolute extrema with non-zero adjustment values */
2573 /* also use this code block. */
2574
2575 if (s_cmp(uop, "<", (ftnlen)6, (ftnlen)1) == 0 || s_cmp(uop, ">", (ftnlen)
2576 6, (ftnlen)1) == 0) {
2577
2578 /* We'll find the window when the coordinate is less than */
2579 /* the reference value. If the relation is '>', we'll */
2580 /* complement the result. Let FLIP indicate whether */
2581 /* we need to take the complement of our result at the */
2582 /* end of the search. */
2583
2584 if (s_cmp(uop, ">", (ftnlen)6, (ftnlen)1) == 0) {
2585 s_copy(uop, "<", (ftnlen)6, (ftnlen)1);
2586 flip = TRUE_;
2587 } else {
2588 flip = FALSE_;
2589 }
2590
2591 /* We'll need the sine and cosine of the reference value. */
2592
2593 cv = cos(xrfval);
2594 sv = sin(xrfval);
2595
2596 /* Determine the quadrant QUAD of the reference value. */
2597
2598 locref = atan2(sv, cv);
2599 if (locref < -pi_() / 2) {
2600 quad = 3;
2601 } else if (locref < 0.) {
2602 quad = 4;
2603 } else if (locref < pi_() / 2) {
2604 quad = 1;
2605 } else {
2606 quad = 2;
2607 }
2608
2609 /* Create a list of region windows to compute. The order */
2610 /* of list items is significant: the regions will */
2611 /* be computed in the order in which they're listed. */
2612
2613 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2614 nl = 2;
2615 s_copy(rlist, "Q2", (ftnlen)32, (ftnlen)2);
2616 s_copy(rlist + 32, "Q3", (ftnlen)32, (ftnlen)2);
2617 } else {
2618 nl = 3;
2619 s_copy(rlist, "LEFT", (ftnlen)32, (ftnlen)4);
2620 s_copy(rlist + 32, "Q1", (ftnlen)32, (ftnlen)2);
2621 s_copy(rlist + 64, "Q4", (ftnlen)32, (ftnlen)2);
2622 }
2623
2624 /* Compute all of the region windows. */
2625
2626 /* We make use of the fact that windows TOP and RIGHT */
2627 /* have already been computed. */
2628
2629 i__1 = nl;
2630 for (i__ = 1; i__ <= i__1; ++i__) {
2631 if (s_cmp(rlist + (((i__2 = i__ - 1) < 7 && 0 <= i__2 ? i__2 :
2632 s_rnge("rlist", i__2, "zzgflong_", (ftnlen)2528)) << 5),
2633 "LEFT", (ftnlen)32, (ftnlen)4) == 0 && left == 0) {
2634 lnkan_(wwpool, &node);
2635 lnkila_(&head, &node, wwpool);
2636 left = wix[(i__2 = node - 1) < 7 && 0 <= i__2 ? i__2 : s_rnge(
2637 "wix", i__2, "zzgflong_", (ftnlen)2532)];
2638 wndifd_(cnfine, &work[(i__2 = right * work_dim1 - 5 -
2639 work_offset) < work_dim1 * work_dim2 && 0 <= i__2 ?
2640 i__2 : s_rnge("work", i__2, "zzgflong_", (ftnlen)2534)
2641 ], &work[(i__3 = left * work_dim1 - 5 - work_offset) <
2642 work_dim1 * work_dim2 && 0 <= i__3 ? i__3 : s_rnge(
2643 "work", i__3, "zzgflong_", (ftnlen)2534)]);
2644 } else if (s_cmp(rlist + (((i__2 = i__ - 1) < 7 && 0 <= i__2 ?
2645 i__2 : s_rnge("rlist", i__2, "zzgflong_", (ftnlen)2536))
2646 << 5), "Q1", (ftnlen)32, (ftnlen)2) == 0 && q1 == 0) {
2647 if (q1 == 0) {
2648 lnkan_(wwpool, &node);
2649 lnkila_(&head, &node, wwpool);
2650 q1 = wix[(i__2 = node - 1) < 7 && 0 <= i__2 ? i__2 :
2651 s_rnge("wix", i__2, "zzgflong_", (ftnlen)2542)];
2652 }
2653 wnintd_(&work[(i__2 = right * work_dim1 - 5 - work_offset) <
2654 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge(
2655 "work", i__2, "zzgflong_", (ftnlen)2546)], &work[(
2656 i__3 = top * work_dim1 - 5 - work_offset) < work_dim1
2657 * work_dim2 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
2658 "zzgflong_", (ftnlen)2546)], &work[(i__4 = q1 *
2659 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2
2660 && 0 <= i__4 ? i__4 : s_rnge("work", i__4, "zzgflong_"
2661 , (ftnlen)2546)]);
2662 } else if (s_cmp(rlist + (((i__2 = i__ - 1) < 7 && 0 <= i__2 ?
2663 i__2 : s_rnge("rlist", i__2, "zzgflong_", (ftnlen)2549))
2664 << 5), "Q2", (ftnlen)32, (ftnlen)2) == 0 && q2 == 0) {
2665 lnkan_(wwpool, &node);
2666 lnkila_(&head, &node, wwpool);
2667 q2 = wix[(i__2 = node - 1) < 7 && 0 <= i__2 ? i__2 : s_rnge(
2668 "wix", i__2, "zzgflong_", (ftnlen)2553)];
2669 wndifd_(&work[(i__2 = top * work_dim1 - 5 - work_offset) <
2670 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge(
2671 "work", i__2, "zzgflong_", (ftnlen)2555)], &work[(
2672 i__3 = right * work_dim1 - 5 - work_offset) <
2673 work_dim1 * work_dim2 && 0 <= i__3 ? i__3 : s_rnge(
2674 "work", i__3, "zzgflong_", (ftnlen)2555)], &work[(
2675 i__4 = q2 * work_dim1 - 5 - work_offset) < work_dim1 *
2676 work_dim2 && 0 <= i__4 ? i__4 : s_rnge("work", i__4,
2677 "zzgflong_", (ftnlen)2555)]);
2678 } else if (s_cmp(rlist + (((i__2 = i__ - 1) < 7 && 0 <= i__2 ?
2679 i__2 : s_rnge("rlist", i__2, "zzgflong_", (ftnlen)2558))
2680 << 5), "Q3", (ftnlen)32, (ftnlen)2) == 0 && q3 == 0) {
2681
2682 /* Note: we need the bottom window in order to compute Q3! */
2683
2684 lnkan_(wwpool, &node);
2685 lnkila_(&head, &node, wwpool);
2686 q3 = wix[(i__2 = node - 1) < 7 && 0 <= i__2 ? i__2 : s_rnge(
2687 "wix", i__2, "zzgflong_", (ftnlen)2564)];
2688 wndifd_(&work[(i__2 = bot * work_dim1 - 5 - work_offset) <
2689 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge(
2690 "work", i__2, "zzgflong_", (ftnlen)2566)], &work[(
2691 i__3 = right * work_dim1 - 5 - work_offset) <
2692 work_dim1 * work_dim2 && 0 <= i__3 ? i__3 : s_rnge(
2693 "work", i__3, "zzgflong_", (ftnlen)2566)], &work[(
2694 i__4 = q3 * work_dim1 - 5 - work_offset) < work_dim1 *
2695 work_dim2 && 0 <= i__4 ? i__4 : s_rnge("work", i__4,
2696 "zzgflong_", (ftnlen)2566)]);
2697 } else if (s_cmp(rlist + (((i__2 = i__ - 1) < 7 && 0 <= i__2 ?
2698 i__2 : s_rnge("rlist", i__2, "zzgflong_", (ftnlen)2569))
2699 << 5), "Q4", (ftnlen)32, (ftnlen)2) == 0 && q4 == 0) {
2700
2701 /* NOTE: We need the bottom window in order to compute Q4! */
2702
2703 lnkan_(wwpool, &node);
2704 lnkila_(&head, &node, wwpool);
2705 q4 = wix[(i__2 = node - 1) < 7 && 0 <= i__2 ? i__2 : s_rnge(
2706 "wix", i__2, "zzgflong_", (ftnlen)2575)];
2707 wnintd_(&work[(i__2 = right * work_dim1 - 5 - work_offset) <
2708 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge(
2709 "work", i__2, "zzgflong_", (ftnlen)2577)], &work[(
2710 i__3 = bot * work_dim1 - 5 - work_offset) < work_dim1
2711 * work_dim2 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
2712 "zzgflong_", (ftnlen)2577)], &work[(i__4 = q4 *
2713 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2
2714 && 0 <= i__4 ? i__4 : s_rnge("work", i__4, "zzgflong_"
2715 , (ftnlen)2577)]);
2716 }
2717 }
2718 if (failed_()) {
2719 chkout_("ZZGFLONG", (ftnlen)8);
2720 return 0;
2721 }
2722
2723 /* Now decide the sector and proxy function we'll use to */
2724 /* search for the time when the reference value is hit. */
2725
2726 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2727 if (quad == 1) {
2728 s = right;
2729 s_copy(prxfun, "LONGITUDE", (ftnlen)50, (ftnlen)9);
2730 } else if (quad == 2) {
2731 s = q2;
2732 s_copy(prxfun, "RIGHT ASCENSION", (ftnlen)50, (ftnlen)15);
2733 } else if (quad == 3) {
2734 s = q3;
2735 s_copy(prxfun, "RIGHT ASCENSION", (ftnlen)50, (ftnlen)15);
2736 } else {
2737 s = right;
2738 s_copy(prxfun, "LONGITUDE", (ftnlen)50, (ftnlen)9);
2739 }
2740 } else {
2741 if (quad == 1) {
2742 s = q1;
2743 s_copy(prxfun, "LONGITUDE", (ftnlen)50, (ftnlen)9);
2744 } else if (quad == 2) {
2745 s = left;
2746 s_copy(prxfun, "RIGHT ASCENSION", (ftnlen)50, (ftnlen)15);
2747 } else if (quad == 3) {
2748 s = left;
2749 s_copy(prxfun, "RIGHT ASCENSION", (ftnlen)50, (ftnlen)15);
2750 } else {
2751 s = q4;
2752 s_copy(prxfun, "LONGITUDE", (ftnlen)50, (ftnlen)9);
2753 }
2754 }
2755
2756 /* Set the proxy reference value based on the input */
2757 /* reference value and the choice of proxy function. */
2758
2759 if (s_cmp(prxfun, "LONGITUDE", (ftnlen)50, (ftnlen)9) == 0) {
2760 prxval = atan2(sv, cv);
2761 } else {
2762 prxval = atan2(sv, cv);
2763 if (prxval < 0.) {
2764 prxval += twopi_();
2765 }
2766 }
2767
2768 /* We're going to need additional windows in order to search */
2769 /* quadrant Q. At this point, we're going to de-allocate all */
2770 /* windows except those needed for the upcoming searches. */
2771
2772 /* Create the set NEEDWN of the windows we need to retain. */
2773
2774 ssizei_(&c__7, needwn);
2775 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2776 insrti_(&q2, needwn);
2777 insrti_(&q3, needwn);
2778 insrti_(&right, needwn);
2779 } else {
2780 insrti_(&q1, needwn);
2781 insrti_(&q4, needwn);
2782 insrti_(&left, needwn);
2783 }
2784
2785 /* Now delete all windows not referenced by NEEDWN. */
2786
2787 node = head;
2788 while(node > 0) {
2789
2790 /* Find the next node in the list. */
2791
2792 next = lnknxt_(&node, wwpool);
2793 if (! elemi_(&wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 :
2794 s_rnge("wix", i__1, "zzgflong_", (ftnlen)2694)], needwn))
2795 {
2796
2797 /* Delete NODE; update HEAD if we deleted the head node. */
2798
2799 lnkfsl_(&node, &node, wwpool);
2800 if (head == node) {
2801 head = next;
2802 }
2803 }
2804
2805 /* Prepare to look at the next node. */
2806
2807 node = next;
2808 }
2809 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2810
2811 /* This is a longitude search. */
2812
2813 /* For each quadrant, identify or compute the window on which */
2814 /* the constraint is automatically satisfied. Store the result */
2815 /* in workspace window F1. If this window is empty, set F1 to */
2816 /* 0. */
2817
2818 if (quad == 1) {
2819 f1 = q3;
2820 } else if (quad == 2) {
2821 lnkan_(wwpool, &node);
2822 lnkila_(&head, &node, wwpool);
2823 f1 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge(
2824 "wix", i__1, "zzgflong_", (ftnlen)2731)];
2825 wnunid_(&work[(i__1 = q3 * work_dim1 - 5 - work_offset) <
2826 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2827 "work", i__1, "zzgflong_", (ftnlen)2733)], &work[(
2828 i__2 = right * work_dim1 - 5 - work_offset) <
2829 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge(
2830 "work", i__2, "zzgflong_", (ftnlen)2733)], &work[(
2831 i__3 = f1 * work_dim1 - 5 - work_offset) < work_dim1 *
2832 work_dim2 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
2833 "zzgflong_", (ftnlen)2733)]);
2834 } else if (quad == 3) {
2835 f1 = 0;
2836 } else {
2837
2838 /* QUAD is 4. */
2839
2840 f1 = q3;
2841 }
2842 } else {
2843
2844 /* We're working with RA. */
2845
2846 if (quad == 1) {
2847 f1 = 0;
2848 } else if (quad == 2) {
2849 f1 = q1;
2850 } else if (quad == 3) {
2851 f1 = q1;
2852 } else {
2853
2854 /* QUAD is 4. */
2855
2856 lnkan_(wwpool, &node);
2857 lnkila_(&head, &node, wwpool);
2858 f1 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge(
2859 "wix", i__1, "zzgflong_", (ftnlen)2770)];
2860 wnunid_(&work[(i__1 = left * work_dim1 - 5 - work_offset) <
2861 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2862 "work", i__1, "zzgflong_", (ftnlen)2772)], &work[(
2863 i__2 = q1 * work_dim1 - 5 - work_offset) < work_dim1 *
2864 work_dim2 && 0 <= i__2 ? i__2 : s_rnge("work", i__2,
2865 "zzgflong_", (ftnlen)2772)], &work[(i__3 = f1 *
2866 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2
2867 && 0 <= i__3 ? i__3 : s_rnge("work", i__3, "zzgflong_"
2868 , (ftnlen)2772)]);
2869 }
2870 }
2871 if (failed_()) {
2872 chkout_("ZZGFLONG", (ftnlen)8);
2873 return 0;
2874 }
2875
2876 /* Search sector S to find times when the relation */
2877
2878 /* PRXFUN PRXREL PRXVAL */
2879
2880 /* holds. */
2881
2882 /* Allocate window F2 to hold the result of the search. */
2883
2884
2885 for (i__ = 1; i__ <= 2; ++i__) {
2886 i__2 = total + i__;
2887 repmi_(tmplat, "#", &i__2, rptpre + ((i__1 = i__ - 1) < 2 && 0 <=
2888 i__1 ? i__1 : s_rnge("rptpre", i__1, "zzgflong_", (ftnlen)
2889 2794)) * 80, (ftnlen)80, (ftnlen)1, (ftnlen)80);
2890 }
2891 lnkan_(wwpool, &node);
2892 lnkila_(&head, &node, wwpool);
2893 f2 = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2894 i__1, "zzgflong_", (ftnlen)2800)];
2895 scardd_(&c__0, &work[(i__1 = f2 * work_dim1 - 5 - work_offset) <
2896 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2897 i__1, "zzgflong_", (ftnlen)2802)]);
2898 if (s_cmp(prxfun, "LONGITUDE", (ftnlen)50, (ftnlen)9) == 0) {
2899
2900 /* Initialize the proxy search in sector S, then perform the */
2901 /* search. */
2902
2903 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
2904 "LATITUDINAL", "LONGITUDE", vecdef_len, method_len,
2905 target_len, ref_len, abcorr_len, obsrvr_len, dref_len, (
2906 ftnlen)11, (ftnlen)9);
2907 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcodc_, (U_fp)
2908 zzgfudlt_, (S_fp)zzgfcog_, "<", &prxval, &loctol, &c_b70,
2909 &work[(i__1 = s * work_dim1 - 5 - work_offset) <
2910 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2911 i__1, "zzgflong_", (ftnlen)2814)], mw, nw, work, rpt, (
2912 U_fp)udrepi, (U_fp)udrepu, (U_fp)udrepf, rptpre, rptsuf,
2913 bail, (L_fp)udbail, &work[(i__2 = f2 * work_dim1 - 5 -
2914 work_offset) < work_dim1 * work_dim2 && 0 <= i__2 ? i__2 :
2915 s_rnge("work", i__2, "zzgflong_", (ftnlen)2814)], (
2916 ftnlen)1, (ftnlen)80, (ftnlen)80);
2917 } else {
2918
2919 /* Initialize the proxy search in sector S, then perform the */
2920 /* search. */
2921
2922 zzgfcoin_(vecdef, method, target, ref, abcorr, obsrvr, dref, dvec,
2923 "RA/DEC", "RIGHT ASCENSION", vecdef_len, method_len,
2924 target_len, ref_len, abcorr_len, obsrvr_len, dref_len, (
2925 ftnlen)6, (ftnlen)15);
2926 zzgfrelx_((U_fp)udstep, (U_fp)udrefn, (U_fp)zzgfcodc_, (U_fp)
2927 zzgfudlt_, (S_fp)zzgfcog_, "<", &prxval, &loctol, &c_b70,
2928 &work[(i__1 = s * work_dim1 - 5 - work_offset) <
2929 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2930 i__1, "zzgflong_", (ftnlen)2833)], mw, nw, work, rpt, (
2931 U_fp)udrepi, (U_fp)udrepu, (U_fp)udrepf, rptpre, rptsuf,
2932 bail, (L_fp)udbail, &work[(i__2 = f2 * work_dim1 - 5 -
2933 work_offset) < work_dim1 * work_dim2 && 0 <= i__2 ? i__2 :
2934 s_rnge("work", i__2, "zzgflong_", (ftnlen)2833)], (
2935 ftnlen)1, (ftnlen)80, (ftnlen)80);
2936 }
2937
2938 /* 7 + 0:2 passes done for adjusted extrema. */
2939
2940 if (*bail) {
2941 if ((*udbail)()) {
2942 chkout_("ZZGFLONG", (ftnlen)8);
2943 return 0;
2944 }
2945 }
2946
2947 /* Combine the contents of windows F1 and F2 to obtain */
2948 /* the result. */
2949
2950 if (f1 != 0) {
2951 wnunid_(&work[(i__1 = f1 * work_dim1 - 5 - work_offset) <
2952 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2953 i__1, "zzgflong_", (ftnlen)2860)], &work[(i__2 = f2 *
2954 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2 && 0
2955 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
2956 ftnlen)2860)], result);
2957 } else {
2958 copyd_(&work[(i__1 = f2 * work_dim1 - 5 - work_offset) <
2959 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
2960 i__1, "zzgflong_", (ftnlen)2864)], result);
2961 }
2962
2963 /* Last step: complement the result if necessary. */
2964
2965 if (flip) {
2966
2967 /* Create the window relative to which we'll find */
2968 /* the complement of RESULT. The window we seek */
2969 /* is not CNFINE, but rather a union of windows */
2970 /* that avoids the branch cut. */
2971
2972 lnkan_(wwpool, &node);
2973 wh = wix[(i__1 = node - 1) < 7 && 0 <= i__1 ? i__1 : s_rnge("wix",
2974 i__1, "zzgflong_", (ftnlen)2880)];
2975 if (s_cmp(nrmcrd, "LONGITUDE", (ftnlen)32, (ftnlen)9) == 0) {
2976 wnunid_(&work[(i__1 = q2 * work_dim1 - 5 - work_offset) <
2977 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2978 "work", i__1, "zzgflong_", (ftnlen)2884)], &work[(
2979 i__2 = right * work_dim1 - 5 - work_offset) <
2980 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge(
2981 "work", i__2, "zzgflong_", (ftnlen)2884)], &work[(
2982 i__3 = f2 * work_dim1 - 5 - work_offset) < work_dim1 *
2983 work_dim2 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
2984 "zzgflong_", (ftnlen)2884)]);
2985 wnunid_(&work[(i__1 = q3 * work_dim1 - 5 - work_offset) <
2986 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2987 "work", i__1, "zzgflong_", (ftnlen)2885)], &work[(
2988 i__2 = f2 * work_dim1 - 5 - work_offset) < work_dim1 *
2989 work_dim2 && 0 <= i__2 ? i__2 : s_rnge("work", i__2,
2990 "zzgflong_", (ftnlen)2885)], &work[(i__3 = wh *
2991 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2
2992 && 0 <= i__3 ? i__3 : s_rnge("work", i__3, "zzgflong_"
2993 , (ftnlen)2885)]);
2994 } else {
2995 wnunid_(&work[(i__1 = q1 * work_dim1 - 5 - work_offset) <
2996 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
2997 "work", i__1, "zzgflong_", (ftnlen)2887)], &work[(
2998 i__2 = left * work_dim1 - 5 - work_offset) <
2999 work_dim1 * work_dim2 && 0 <= i__2 ? i__2 : s_rnge(
3000 "work", i__2, "zzgflong_", (ftnlen)2887)], &work[(
3001 i__3 = f2 * work_dim1 - 5 - work_offset) < work_dim1 *
3002 work_dim2 && 0 <= i__3 ? i__3 : s_rnge("work", i__3,
3003 "zzgflong_", (ftnlen)2887)]);
3004 wnunid_(&work[(i__1 = q4 * work_dim1 - 5 - work_offset) <
3005 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge(
3006 "work", i__1, "zzgflong_", (ftnlen)2888)], &work[(
3007 i__2 = f2 * work_dim1 - 5 - work_offset) < work_dim1 *
3008 work_dim2 && 0 <= i__2 ? i__2 : s_rnge("work", i__2,
3009 "zzgflong_", (ftnlen)2888)], &work[(i__3 = wh *
3010 work_dim1 - 5 - work_offset) < work_dim1 * work_dim2
3011 && 0 <= i__3 ? i__3 : s_rnge("work", i__3, "zzgflong_"
3012 , (ftnlen)2888)]);
3013 }
3014
3015 /* We use F2 as a temporary window index, since F2 is */
3016 /* guaranteed to exist at this point and is distinct from WH. */
3017
3018 wndifd_(&work[(i__1 = wh * work_dim1 - 5 - work_offset) <
3019 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
3020 i__1, "zzgflong_", (ftnlen)2895)], result, &work[(i__2 =
3021 f2 * work_dim1 - 5 - work_offset) < work_dim1 * work_dim2
3022 && 0 <= i__2 ? i__2 : s_rnge("work", i__2, "zzgflong_", (
3023 ftnlen)2895)]);
3024 copyd_(&work[(i__1 = f2 * work_dim1 - 5 - work_offset) <
3025 work_dim1 * work_dim2 && 0 <= i__1 ? i__1 : s_rnge("work",
3026 i__1, "zzgflong_", (ftnlen)2896)], result);
3027 }
3028 }
3029 chkout_("ZZGFLONG", (ftnlen)8);
3030 return 0;
3031 } /* zzgflong_ */
3032
3033