1 /* zztanslv.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__1000 = 1000;
11
12 /* $Procedure ZZTANSLV ( Private --- tangent point solver ) */
zztanslv_(S_fp udcond,S_fp udstep,S_fp udrefn,logical * cstep,doublereal * step,doublereal * start,doublereal * finish,doublereal * tol,doublereal * result,doublereal * points,logical * endflg)13 /* Subroutine */ int zztanslv_(S_fp udcond, S_fp udstep, S_fp udrefn, logical
14 *cstep, doublereal *step, doublereal *start, doublereal *finish,
15 doublereal *tol, doublereal *result, doublereal *points, logical *
16 endflg)
17 {
18 /* System generated locals */
19 integer i__1;
20 doublereal d__1, d__2;
21
22 /* Builtin functions */
23 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
24
25 /* Local variables */
26 integer room;
27 extern /* Subroutine */ int vequ_(doublereal *, doublereal *);
28 doublereal curx, svdx;
29 extern /* Subroutine */ int zzwninsd_(doublereal *, doublereal *, char *,
30 doublereal *, ftnlen);
31 logical s;
32 doublereal begin, t;
33 extern /* Subroutine */ int chkin_(char *, ftnlen), errdp_(char *,
34 doublereal *, ftnlen);
35 extern integer sized_(doublereal *);
36 integer nloop;
37 doublereal xstep, x1, x2;
38 logical state1, state2;
39 extern logical failed_(void);
40 integer to;
41 extern doublereal brcktd_(doublereal *, doublereal *, doublereal *);
42 doublereal maxmag;
43 extern doublereal touchd_(doublereal *);
44 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
45 ftnlen);
46 logical cursta, instat, savsta;
47 extern /* Subroutine */ int setmsg_(char *, ftnlen), errint_(char *,
48 integer *, ftnlen);
49 extern logical return_(void);
50 char contxt[256];
51 doublereal xpoint[3];
52 logical prvset;
53 doublereal trnstn, prvpnt[3];
54
55 /* $ Abstract */
56
57 /* SPICE Private routine intended solely for the support of SPICE */
58 /* routines. Users should not call this routine directly due */
59 /* to the volatile nature of this routine. */
60
61 /* This routine finds tangent points of rays on a target surface, */
62 /* where the rays are confined to a specified half-plane. It may */
63 /* be used for limb and terminator computations. */
64
65 /* $ Disclaimer */
66
67 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
68 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
69 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
70 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
71 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
72 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
73 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
74 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
75 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
76 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
77
78 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
79 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
80 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
81 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
82 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
83 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
84
85 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
86 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
87 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
88 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
89
90 /* $ Required_Reading */
91
92 /* None. */
93
94 /* $ Keywords */
95
96 /* ROOT */
97 /* SEARCH */
98 /* WINDOWS */
99
100 /* $ Declarations */
101 /* $ Brief_I/O */
102
103 /* VARIABLE I/O DESCRIPTION */
104 /* -------- --- -------------------------------------------------- */
105 /* UDCOND I Name of the routine that compares the current */
106 /* state and ray intercept. */
107 /* UDSTEP I Name of the routine that computes a step */
108 /* UDREFN I Name of the routine that computes a refined input. */
109 /* CSTEP I Logical indicating constant step size. */
110 /* STEP I Constant step size for finding geometric events. */
111 /* START I Beginning of the search interval. */
112 /* FINISH I End of the search interval. */
113 /* TOL I Maximum error in detection of state transitions. */
114 /* RESULT I-O SPICE window containing results. */
115 /* POINTS O Array of points associated with transitions. */
116 /* ENDFLG O Endpoint transition flags. */
117
118 /* $ Detailed_Input */
119
120 /* For the purpose of solving for ray tangency points on a target */
121 /* body, the independent variable can be considered to be angular */
122 /* separation of a ray from an axis vector. For a limb computation, */
123 /* the axis vector points from the target body center toward the */
124 /* observer. For a terminator computation, the axis vector points */
125 /* from the target body center towards the center of the */
126 /* illumination source. The "system state" for these computations is */
127 /* the condition of the ray intersecting the target body: if an */
128 /* intersection exists, the state is "true." */
129
130 /* The discussion below is more general; we'll use the terms */
131 /* "abscissa" or "x-value" rather than "angle" to describe the */
132 /* independent variable. The "system state" is simply a boolean */
133 /* function of the independent variable. */
134
135 /* The first three inputs to this routine are names of subroutines */
136 /* that this routine will call. */
137
138 /* These routines should meet the following specifications. */
139
140 /* UDCOND the routine that determines if the system state */
141 /* satisfies some constraint at a given independent */
142 /* variable value X. */
143
144 /* The calling sequence: */
145
146 /* CALL UDCOND ( X, IN_CON, POINT ) */
147
148 /* where: */
149
150 /* X a double precision value at which to */
151 /* evaluate the state. */
152
153 /* IN_CON a logical value indicating whether */
154 /* or not the quantity satisfies the */
155 /* constraint at X (TRUE) or not (FALSE). */
156
157 /* POINT is a 3-vector associated with X. POINT */
158 /* is defined if and only if IN_CON is .TRUE. */
159
160
161 /* UDSTEP the routine that computes a step in an attempt to */
162 /* find a transition of the state of the specified */
163 /* coordinate. In the context of this routine's algorithm, */
164 /* a "state transition" occurs where the geometric state */
165 /* changes from being in the desired geometric condition */
166 /* event to not, or vice versa. */
167
168 /* This routine relies on UDSTEP returning step sizes */
169 /* small enough so that state transitions within the */
170 /* interval [START, FINISH] are not overlooked. There */
171 /* must never be two roots A and B separated by less than */
172 /* STEP, where STEP is the minimum step size returned by */
173 /* UDSTEP for any value of X in the interval [A, B]. */
174
175 /* The calling sequence for UDSTEP is: */
176
177 /* CALL UDSTEP ( X, STEP ) */
178
179 /* where: */
180
181 /* X a double precision value from which the */
182 /* algorithm is to search forward for a state */
183 /* transition. */
184
185 /* STEP is the output step size. STEP indicates how */
186 /* far to advance X so that X and X+STEP may */
187 /* bracket a state transition and definitely */
188 /* do not bracket more than one state */
189 /* transition. */
190
191 /* If a constant step size is desired, the routine */
192
193 /* GFSTEP */
194
195 /* may be used. This is the default option. If using */
196 /* GFSTEP, the step size must be set by calling */
197
198 /* GFSSTP(STEP) */
199
200 /* prior to calling this routine. */
201
202
203 /* UDREFN the routine that computes a refinement in the abscissa */
204 /* values that bracket a transition point. In other */
205 /* words, once a pair of abscissa values have been */
206 /* detected such that the system is in different states */
207 /* at each of the two values, UDREFN selects an */
208 /* intermediate abscissa value which should be closer to */
209 /* the transition state than one of the two known X */
210 /* values. The calling sequence for UDREFN is: */
211
212 /* CALL UDREFN ( X1, X2, S1, S2, T ) */
213
214 /* where the inputs are: */
215
216 /* X1 an X (abscissa) value at which the system is */
217 /* in state S1. */
218
219 /* X2 an X value at which the system is in state */
220 /* S2. X2 is assumed to be larger than X1. */
221
222 /* S1 a logical indicating the state of the system */
223 /* at X1. */
224
225 /* S2 a logical indicating the state of the system */
226 /* at X2. */
227
228 /* UDREFN may use or ignore the S1 and S2 values. */
229
230 /* The output is: */
231
232 /* T an X value to check for a state transition */
233 /* between X1 and X2. */
234
235 /* If a simple bisection method is desired, the routine */
236 /* GFREFN may be used. This is the default option. */
237
238 /* CSTEP is a logical indicating whether or not the step size */
239 /* used in searching is constant. If it is, the value */
240 /* STEP is used. Note that even if UDSTEP has the value */
241 /* GFSTEP, i.e. the public, constant step routine, CSTEP */
242 /* should still be .FALSE., in which case STEP is ignored. */
243
244 /* STEP is the step size to be used in the search. STEP must */
245 /* be short enough for a search using this step size to */
246 /* locate the intervals where the geometric event */
247 /* function is monotone increasing or decreasing. */
248 /* However, STEP must not be *too* short, or the search */
249 /* will take an unreasonable amount of time. */
250
251 /* The choice of STEP affects the completeness but not */
252 /* the precision of solutions found by this routine; */
253 /* precision is controlled by the convergence the */
254 /* tolerance, TOL. */
255
256 /* START is the beginning of the interval over which the state */
257 /* is to be detected. */
258
259 /* FINISH is the end of the interval over which the state is */
260 /* to be detected. */
261
262 /* TOL is a tolerance value used to determine convergence of */
263 /* root-finding operations. TOL is measured in the units */
264 /* of the independent variable and is greater than zero. */
265
266 /* RESULT is an initialized SPICE window. RESULT must be large */
267 /* enough to hold all of the intervals found by the */
268 /* search. */
269
270 /* $ Detailed_Output */
271
272 /* RESULT is a SPICE window containing the results of the */
273 /* search. With the exception of the first and last */
274 /* endpoints of the window, the endpoints of the */
275 /* intervals in RESULT always are abscissa values of */
276 /* state transitions. The first and last endpoints may or */
277 /* may not correspond to state transitions. */
278
279 /* The first and last endpoints are abscissa values of */
280 /* state transitions if and only if the state function is */
281 /* .FALSE. at those points. */
282
283 /* Note that, in the special case where the state function */
284 /* is .TRUE. at the first endpoint and .FALSE. in a */
285 /* half-neighborhood to the right of that endpoint, it is */
286 /* possible for this function to find a transition at that */
287 /* endpoint and assign it to the right endpoint of the */
288 /* degenerate interval */
289
290 /* [ left endpoint, left endpoint ] */
291
292 /* Analogously, it is possible to find a state transition */
293 /* at the last endpoint if the state function is .TRUE. */
294 /* at that endpoint and .FALSE. in a half-neighborhood */
295 /* to the left of that point. */
296
297 /* The output ENDFLG indicates whether the first and */
298 /* last endpoints of RESULT are transitions. */
299
300
301 /* POINTS is an array of 3-vectors associated with the endpoints */
302 /* of the intervals of RESULT. Elements */
303
304 /* POINTS(J,I), J = 1 .. 3 */
305
306 /* constitute a vector associated with */
307
308 /* RESULT(I) */
309
310 /* The first and last vectors of POINTS are valid if */
311 /* and only if the corresponding elements of ENDFLG */
312 /* are .TRUE. */
313
314 /* Presuming this routine is used to solve for tangent */
315 /* points on a target body, the vectors contained in */
316 /* POINTS, when valid, are such tangent points. */
317
318 /* POINTS must be declared with size at least 3 times */
319 /* that of RESULT. */
320
321
322 /* ENDFLG is an array of two logical flags that indicate */
323 /* whether state transitions occur at the initial */
324 /* and final endpoints of the result window. Element */
325 /* 1 of this array is .TRUE. if and only if */
326 /* there is a state transition at the first endpoint of */
327 /* RESULT (in element RESULT(1)); element 2 is .TRUE. */
328 /* if and only if there is a state transition at the */
329 /* last element of RESULT. */
330
331
332 /* $ Parameters */
333
334 /* LBCELL is the SPICELIB cell lower bound. */
335
336 /* $ Exceptions */
337
338 /* 1) If TOL is negative, the error SPICE(INVALIDTOLERANCE) */
339 /* will be signaled. */
340
341 /* 2) If START +/- TOL is indistinguishable from START or */
342 /* FINISH +/- TOL is indistinguishable from FINISH, the */
343 /* error SPICE(INVALIDTOLERANCE) will be signaled. */
344
345 /* 3) If START is greater than FINISH, the error */
346 /* SPICE(BOUNDSOUTOFORDER) will be signaled. */
347
348 /* 4) If a constant step is used, the step must be positive and */
349 /* have magnitude large enough so that, when added to a value in */
350 /* the range [START, FINISH], it will yield a distinct value. */
351 /* Otherwise, the error SPICE(INVALIDCONSTSTEP) will be */
352 /* signaled. */
353
354 /* 5) If the step function is used and it returns a value less than */
355 /* the preceding value, the error SPICE(INVALIDSTEP) will */
356 /* be signaled. */
357
358 /* 6) If the inner convergence loop fails to converge to TOL within */
359 /* MXLOOP iterations, the error SPICE(NOCONVERGENCE) will be */
360 /* signaled. */
361
362 /* 7) If the POINTS array doesn't have enough room to store */
363 /* the points associated with state transitions, the error */
364 /* SPICE(ARRAYTOOSMALL) will be signaled. */
365
366 /* 8) If the result window doesn't have enough room to store */
367 /* the abscissas associated with state transitions, the error */
368 /* will be signaled by a routine in the call tree of this */
369 /* routine. */
370
371 /* $ Files */
372
373 /* Kernels used by this routine are those needed by the input */
374 /* routines */
375
376 /* UDCOND */
377 /* UDGETP */
378 /* UDSTEP */
379 /* UDREFN */
380
381 /* $ Particulars */
382
383 /* This routine supports limb and terminator point detection. */
384
385 /* $ Examples */
386
387 /* None. */
388
389 /* $ Restrictions */
390
391 /* It is important that the user understand how the routines UDCOND, */
392 /* UDSTEP and UDREFN are to be used and that the calling sequences */
393 /* match precisely with the descriptions given here. */
394
395 /* $ Literature_References */
396
397 /* None. */
398
399 /* $ Author_and_Institution */
400
401 /* N.J. Bachman (JPL) */
402 /* L.S. Elson (JPL) */
403 /* W.L. Taber (JPL) */
404 /* I.M. Underwood (JPL) */
405 /* E.D. Wright (JPL) */
406
407 /* $ Version */
408
409 /* - SPICELIB Version 1.0.0, 30-JUN-2016 (NJB) (EDW) */
410
411 /* Updated logic for placement of points in output POINTS */
412 /* array. If the first element of RESULT is equal to BEGIN, */
413 /* then space will be reserved in the first element of POINTS, */
414 /* so as to keep the output points synced with the elements */
415 /* of RESULT. */
416
417 /* Updated short error messages. */
418
419 /* Updated header I/O sections. */
420
421 /* 12-FEB-2016 (NJB) (EDW) */
422
423 /* Derived from ZZGFSOLV Version 1.1.0, 24-OCT-2010 (EDW) */
424
425 /* -& */
426 /* $ Index_Entries */
427
428 /* find tangent points on target */
429
430 /* -& */
431
432 /* SPICELIB functions */
433
434
435 /* Local variables */
436
437
438 /* The maximum number of search loop iterations to execute. */
439 /* The default refinement method is bisection, a very slow */
440 /* method to convergence. Since 2**1000 ~ 10**301, */
441 /* 1000 loop iterations represents enough effort to assume */
442 /* either the search will not converge or that the refinement */
443 /* function operates slower than would bisection, in which */
444 /* case the user should use the default GFREFN function. */
445
446
447 /* Standard SPICE error handling. */
448
449 if (return_()) {
450 return 0;
451 }
452 chkin_("ZZTANSLV", (ftnlen)8);
453
454 /* Check the convergence tolerance. */
455
456 if (*tol <= 0.) {
457 setmsg_("Tolerance must be positive but was #.", (ftnlen)37);
458 errdp_("#", tol, (ftnlen)1);
459 sigerr_("SPICE(INVALIDTOLERANCE)", (ftnlen)23);
460 chkout_("ZZTANSLV", (ftnlen)8);
461 return 0;
462 }
463
464 /* Make sure that START is not greater than FINISH. Signal an */
465 /* error for START > FINISH. */
466
467 if (*start > *finish) {
468 setmsg_("Bad input interval: START = # > FINISH = #.", (ftnlen)43);
469 errdp_("#", start, (ftnlen)1);
470 errdp_("#", finish, (ftnlen)1);
471 sigerr_("SPICE(BOUNDSOUTOFORDER)", (ftnlen)23);
472 chkout_("ZZTANSLV", (ftnlen)8);
473 return 0;
474 }
475
476 /* Make sure that TOL is not too small, i.e. that neither */
477 /* START + TOL nor START - TOL equals START. */
478
479 d__1 = *start - *tol;
480 d__2 = *start + *tol;
481 if (touchd_(&d__1) == *start || touchd_(&d__2) == *start) {
482 setmsg_("TOL has value #1. This value is too small to distinguish ST"
483 "ART - TOL or START + TOL from START, #2.", (ftnlen)99);
484 errdp_("#1", tol, (ftnlen)2);
485 errdp_("#2", start, (ftnlen)2);
486 sigerr_("SPICE(INVALIDTOLERANCE)", (ftnlen)23);
487 chkout_("ZZTANSLV", (ftnlen)8);
488 return 0;
489 }
490
491 /* Make sure that TOL is not too small, i.e. that neither */
492 /* FINISH + TOL nor FINISH - TOL equals FINISH. */
493
494 d__1 = *finish - *tol;
495 d__2 = *finish + *tol;
496 if (touchd_(&d__1) == *finish || touchd_(&d__2) == *finish) {
497 setmsg_("TOL has value #1. This value is too small to distinguish FI"
498 "NISH - TOL or FINISH + TOL from FINISH, #2.", (ftnlen)102);
499 errdp_("#1", tol, (ftnlen)2);
500 errdp_("#2", finish, (ftnlen)2);
501 sigerr_("SPICE(INVALIDTOLERANCE)", (ftnlen)23);
502 chkout_("ZZTANSLV", (ftnlen)8);
503 return 0;
504 }
505
506 /* Make sure that STEP is not too small: it must be greater */
507 /* than TOL. */
508
509 if (*cstep) {
510 if (*step <= 0.) {
511 setmsg_("STEP has value #1. The search step must be positive.", (
512 ftnlen)52);
513 errdp_("#1", step, (ftnlen)2);
514 sigerr_("SPICE(INVALIDCONSTSTEP)", (ftnlen)23);
515 chkout_("ZZTANSLV", (ftnlen)8);
516 return 0;
517 }
518 /* Computing MAX */
519 d__1 = abs(*start), d__2 = abs(*finish);
520 maxmag = max(d__1,d__2);
521 d__1 = maxmag + *step;
522 if (touchd_(&d__1) == maxmag) {
523 setmsg_("STEP has value #1. This value is too small to guarantee"
524 " that the search will advance.", (ftnlen)85);
525 errdp_("#1", step, (ftnlen)2);
526 sigerr_("SPICE(INVALIDCONSTSTEP)", (ftnlen)23);
527 chkout_("ZZTANSLV", (ftnlen)8);
528 return 0;
529 }
530 }
531
532 /* This algorithm determines those intervals when a given state is */
533 /* observed to occur within a specified search interval. */
534
535 /* Pairs of X values are recorded. The first member of each pair */
536 /* denotes the X value at which the system changes to the state of */
537 /* interest. The second denotes a transition out of that state. */
538
539 /* If the state is .TRUE. at the beginning of the interval, the */
540 /* beginning of the X interval will be recorded. This may or may not */
541 /* be a transition point. */
542
543 /* Similarly if the state is .TRUE. at the end of the interval, the */
544 /* end of the interval will be recorded. Again, this may or may not */
545 /* be a transition point. */
546
547 /* Initially the current X value is the beginning of the search */
548 /* interval. */
549
550 curx = *start;
551 to = 1;
552 room = sized_(result);
553 prvset = FALSE_;
554
555 /* Determine if the state at the current X value satisfies the */
556 /* constraint. */
557
558 (*udcond)(&curx, &cursta, xpoint);
559 if (failed_()) {
560 chkout_("ZZTANSLV", (ftnlen)8);
561 return 0;
562 }
563 if (cursta) {
564 vequ_(xpoint, prvpnt);
565 prvset = TRUE_;
566 }
567
568 /* If the system is in the state of interest, record the initial */
569 /* X value of the search interval. The variable BEGIN will be */
570 /* used to store the starting point of an interval over which */
571 /* the state is .TRUE. */
572
573 if (cursta) {
574 instat = TRUE_;
575 begin = curx;
576 endflg[0] = FALSE_;
577
578 /* BEGIN will be the first element of RESULT, presuming */
579 /* a state transition is found later. We'll shift the */
580 /* pointer for the output point so the Ith point will */
581 /* correspond to the Ith element of RESULT. */
582
583 /* We don't have to check ROOM yet because we're not */
584 /* inserting anything into POINTS. */
585
586 ++to;
587 --room;
588 } else {
589 instat = FALSE_;
590 endflg[0] = TRUE_;
591 }
592
593 /* If the step size is constant, use the value supplied. */
594
595 if (*cstep) {
596 xstep = *step;
597 }
598
599 /* Save the current X value and state. */
600
601 svdx = curx;
602 savsta = cursta;
603
604 /* Once initializations have been performed keep working */
605 /* until the search interval has been exhausted. */
606
607 /* While the last X value precedes the end of the interval: */
608
609 while(svdx < *finish) {
610
611 /* Attempt to bracket a state change. */
612
613 /* Using the current window and internally stored information */
614 /* about the current state, select a new current X. */
615
616 if (! (*cstep)) {
617 (*udstep)(&curx, &xstep);
618 if (failed_()) {
619 chkout_("ZZTANSLV", (ftnlen)8);
620 return 0;
621 }
622 }
623
624 /* Add the X step to the current X. Make sure that the */
625 /* X does not move beyond the end of the search interval. */
626
627 /* Computing MIN */
628 d__2 = curx + xstep;
629 d__1 = touchd_(&d__2);
630 curx = min(d__1,*finish);
631
632 /* Compute the state at CURX. */
633
634 (*udcond)(&curx, &cursta, xpoint);
635 if (failed_()) {
636 chkout_("ZZTANSLV", (ftnlen)8);
637 return 0;
638 }
639 if (cursta) {
640 vequ_(xpoint, prvpnt);
641 prvset = TRUE_;
642 }
643
644 /* While the state remains unchanged and the interval has not */
645 /* been completely searched ... */
646
647 while(savsta == cursta && svdx < *finish) {
648
649 /* Save the current X and state. */
650
651 svdx = curx;
652 savsta = cursta;
653
654 /* Compute a new current X so that we will not step */
655 /* past the end of the interval. */
656
657 if (! (*cstep)) {
658 (*udstep)(&curx, &xstep);
659 if (failed_()) {
660 chkout_("ZZTANSLV", (ftnlen)8);
661 return 0;
662 }
663 }
664 /* Computing MIN */
665 d__2 = curx + xstep;
666 d__1 = touchd_(&d__2);
667 curx = min(d__1,*finish);
668
669 /* Compute the current state. */
670
671 (*udcond)(&curx, &cursta, xpoint);
672 if (failed_()) {
673 chkout_("ZZTANSLV", (ftnlen)8);
674 return 0;
675 }
676 if (cursta) {
677
678 /* Save the associated vector for the X value CURX. In */
679 /* normal usage, XPOINT is a surface intercept point. */
680
681 vequ_(xpoint, prvpnt);
682 prvset = TRUE_;
683 }
684
685 /* Loop back to see if the state has changed. */
686
687 }
688
689 /* At this point, SVDX and CURX are the X-values at the previous */
690 /* and latest steps, respectively. SAVSTA and CURSTA are the */
691 /* states at these X-values, respectively. */
692
693 /* If we have detected a state change and not merely run out */
694 /* of the search interval... */
695
696 if (savsta != cursta) {
697
698 /* Call the previous state STATE1. */
699 /* Call the current state STATE2. */
700
701 /* Let X1 be the X value at state STATE1. */
702 /* Let X2 be the X value at state STATE2. */
703
704 /* Save the current X. */
705
706 state1 = savsta;
707 state2 = cursta;
708 x1 = svdx;
709 x2 = curx;
710
711 /* Make sure that X1 is not greater than X2. Signal an */
712 /* error for X1 > X2. */
713
714 if (x1 > x2) {
715 setmsg_("Bad x interval result: X1 = # > X2 = #.", (ftnlen)39)
716 ;
717 errdp_("#", &x1, (ftnlen)1);
718 errdp_("#", &x2, (ftnlen)1);
719 sigerr_("SPICE(INVALIDSTEP)", (ftnlen)18);
720 chkout_("ZZTANSLV", (ftnlen)8);
721 return 0;
722 }
723
724 /* Update the saved X and state values to those on the */
725 /* right side of the bracketing interval. We'll use these */
726 /* values for the next bracketing step after a root is */
727 /* found. */
728
729 svdx = curx;
730 savsta = cursta;
731
732 /* X1 and X2 bracket the X value of transition. Squeeze this */
733 /* interval down until it is less than some tolerance in */
734 /* length. Do it as described below... */
735
736 /* Loop while the difference between the X values X1 and X2 */
737 /* exceeds a specified tolerance. */
738
739 nloop = 0;
740 for(;;) { /* while(complicated condition) */
741 d__1 = x2 - x1;
742 if (!(touchd_(&d__1) > *tol))
743 break;
744 ++nloop;
745
746 /* This loop count error exists to catch pathologies */
747 /* in the refinement function. The default bisection */
748 /* refinement will converge before 1000 iterations if */
749 /* a convergence is numerically possible. Any other */
750 /* refinement function should require fewer iterations */
751 /* compared to bisection. If not, the user should */
752 /* probably use bisection. */
753
754 if (nloop >= 1000) {
755 setmsg_("Loop run exceeds maximum loop count. Unable to "
756 "converge to TOL value #1 within MXLOOP value #2 "
757 "iterations.", (ftnlen)106);
758 errdp_("#1", tol, (ftnlen)2);
759 errint_("#2", &c__1000, (ftnlen)2);
760 sigerr_("SPICE(NOCONVERGENCE)", (ftnlen)20);
761 chkout_("ZZTANSLV", (ftnlen)8);
762 return 0;
763 }
764
765 /* Select an X value T, between X1 and X2 (possibly based */
766 /* on the state values). */
767
768 (*udrefn)(&x1, &x2, &state1, &state2, &t);
769
770 /* Check for an error signal. The default refinement */
771 /* routine, GFREFN, does not include error checks. */
772
773 if (failed_()) {
774 chkout_("ZZTANSLV", (ftnlen)8);
775 return 0;
776 }
777
778 /* Check whether T is between X1 and X2. If */
779 /* not then assume that we have gone as far as */
780 /* we can in refining our estimate of the transition */
781 /* point. Set X1 and X2 equal to T. */
782
783 t = brcktd_(&t, &x1, &x2);
784 if (t == x1) {
785
786 /* This assignment may break the invariant that */
787 /* the state at X2 is STATE2. This is allowed */
788 /* because we'll exit the loop immediately. */
789
790 x2 = t;
791 } else if (t == x2) {
792
793 /* This assignment may break the invariant that */
794 /* the state at X1 is STATE1. This is allowed */
795 /* because we'll exit the loop immediately. */
796 x1 = t;
797 } else {
798
799 /* Compute the state at X value T. If this state, S, */
800 /* equals STATE1, set X1 to T, otherwise set X2 to T. */
801
802 (*udcond)(&t, &s, xpoint);
803 if (s) {
804
805 /* Save the latest point associated with a */
806 /* .TRUE. state. */
807
808 vequ_(xpoint, prvpnt);
809 prvset = TRUE_;
810 }
811
812 /* Narrow the interval. Either increase X1 or decrease */
813 /* X2 by setting one of these endpoints to T. Maintain */
814 /* the invariant that the state is STATE1 at X1 and */
815 /* STATE2 at X2. */
816
817 if (s == state1) {
818 x1 = t;
819 } else {
820 x2 = t;
821 }
822 }
823 }
824
825 /* Let TRNSTN be the midpoint of [X1, X2]. Record this */
826 /* abscissa value as marking the transition from STATE1 to */
827 /* STATE2. */
828
829 d__1 = (x1 + x2) * .5;
830 trnstn = brcktd_(&d__1, &x1, &x2);
831
832 /* In state-of-interest or not? INSTAT indicates that STATE1 */
833 /* was .TRUE. We record intervals where the state is .TRUE. */
834 /* when we detect the right hand endpoints of these intervals. */
835
836 if (instat) {
837
838 /* We were in the state of interest. TRNSTN marks the point */
839 /* on the X-axis when the state changed to .FALSE. We need */
840 /* to record the interval from BEGIN to FINISH and note */
841 /* that the state has become .FALSE. */
842
843 /* Add an interval starting at BEGIN and ending at TRNSTN */
844 /* to the result window. */
845
846 s_copy(contxt, "Adding interval [BEGIN,TRNSTN] to RESULT. TR"
847 "NSTN represents time of passage out of the state-of-"
848 "interest.", (ftnlen)256, (ftnlen)105);
849 zzwninsd_(&begin, &trnstn, contxt, result, (ftnlen)256);
850 if (failed_()) {
851 chkout_("ZZTANSLV", (ftnlen)8);
852 return 0;
853 }
854 } else {
855
856 /* The previous state was .FALSE. As a result TRNSTN marks */
857 /* the point where the state becomes .TRUE. Note that we */
858 /* have transitioned to the state of interest and record */
859 /* the X-value at which the transition occurred. */
860
861 begin = trnstn;
862 }
863
864 /* A transition occurred either from from in-state to */
865 /* out-of-state or the inverse. Reverse the value of the */
866 /* INSTAT flag to signify the transition event. */
867
868 instat = ! instat;
869
870 /* For all state transitions, record the last point found */
871 /* by the state function. */
872
873 if (room > 0) {
874
875 /* Add the last point found during the transition search to */
876 /* the POINTS array. */
877
878 if (prvset) {
879 vequ_(prvpnt, &points[to * 3 - 3]);
880 ++to;
881 --room;
882 prvset = FALSE_;
883 } else {
884 setmsg_("PRVPNT should always be set when a transition i"
885 "s detected. We found a transition at #, but PRVS"
886 "ET indicates we don't have a previous point save"
887 "d.", (ftnlen)145);
888 errdp_("#", &trnstn, (ftnlen)1);
889 sigerr_("SPICE(BUG)", (ftnlen)10);
890 chkout_("ZZTANSLV", (ftnlen)8);
891 return 0;
892 }
893 } else {
894
895 /* We ran out of room in the output point array. Note that */
896 /* this error can occur before the result window insertion */
897 /* fails, since that insertion takes place when the state */
898 /* becomes .FALSE. */
899
900 setmsg_("Out of room in the POINTS array. Room is assumed to"
901 " be adequate for SIZED(RESULT) 3-vectors; this size "
902 "is #.", (ftnlen)108);
903 i__1 = sized_(result);
904 errint_("#", &i__1, (ftnlen)1);
905 sigerr_("SPICE(ARRAYTOOSMALL)", (ftnlen)20);
906 chkout_("ZZTANSLV", (ftnlen)8);
907 return 0;
908 }
909
910 /* That's it for this detection of state change. */
911
912 }
913
914 /* Continue if the search interval extends to the right */
915 /* of the latest step. */
916
917 /* SVDX and SAVSTA are already set to the values at the */
918 /* right side of the bracketing interval. */
919
920 }
921
922 /* Check if in-state at this abscissa value (FINISH). INSTAT is the */
923 /* latest state value. If so record the interval. */
924
925 if (instat) {
926
927 /* The state is .TRUE. at FINISH. */
928
929 /* Add an interval starting at BEGIN and ending at FINISH to the */
930 /* window. */
931
932 s_copy(contxt, "Adding interval [BEGIN,FINISH] to RESULT. FINISH rep"
933 "resents end of the search interval.", (ftnlen)256, (ftnlen)87)
934 ;
935 zzwninsd_(&begin, finish, contxt, result, (ftnlen)256);
936 endflg[1] = FALSE_;
937 } else {
938 endflg[1] = TRUE_;
939 }
940 chkout_("ZZTANSLV", (ftnlen)8);
941 return 0;
942 } /* zztanslv_ */
943
944