/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:48 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "ssfind.h" #include void /*FUNCTION*/ ssfind( float xt[], long ix1, long ix2, float x, long *lefti, long *mode) { long int ihi, ilo, middle, step; /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-03-30 SSFIND Krogh Removed Fortran 90 comment. *>> 1994-11-11 CLL Replaced Fortran 90 DO WHILE's with GOTO's for F77. *>> 1994-10-20 SSFIND Krogh Changes to use M77CON *>> 1992-11-12 SSFIND C. L. Lawson, JPL *>> 1992-10-23 C. L. Lawson, JPL *>> 1988-03-16 C. L. Lawson, JPL * * Require IX1 < IX2. * Require the values in XT() indexed from IX1 to IX2 to be * nondecreasing. * Require XT(IX1) < XT(IX1+1). * Require XT(IX2-1) < XT(IX2). * Let A = XT(IX1) and B = XT(IX2). * The closed interval, [A, B], is regarded as being partitioned * partitioned into IX2-IX1 disjoint subintervals, with all but the * last being half-open, and the last one being closed: * [XT(I), XT(I+1)), I = IX1, ..., IX2-2 * [XT(IX2-1), XT(IX2)] * Some of these intervals may have zero length, but not the first * and last ones. * This subroutine identifies the location of X with respect to these * subintervals by setting LEFTI and MODE as follows: * If X is contained in one of these segments of nonzero length, * set MODE = 0 and set LEFTI to be the index of the left end of the * segment containing X. Thus LEFTI will satisfy * IX1 .le. LEFTI .le. IX2-1. * If X < A, set MODE = -1 and LEFTI = IX1. * If X > B, set MODE = +1 and LEFTI = IX2-1. * Issue an error message and stop if * X < XT(IX1+1) and XT(IX1) .ge. XT(IX1+1) or if * X .ge. XT(IX2-1) and XT(IX2-1) .ge. XT(IX2) * ------------------------------------------------------------------ * Method * * Saves the value of LEFTI returned on previous call in ILO. * Starts by checking to see if X is in this same segment. * If so, finishes quickly. We assume this will be a frequently * occurring case. * If not, searches either to the left or right from this previous * segment, as appropriate. During this search the stride is doubled * until a bracketing value is found. Then we finish with a binary * search between the bracketing points. * ------------------------------------------------------------------ * Based on subroutine INTERV on pp. 92-93 of A PRACTICAL GUIDE TO * SPLINES by Carl De Boor, Springer-Verlag, 1978. * Current version by C. L. Lawson, JPL, March 1988. * ------------------------------------------------------------------ * XT() [in] This subr will access XT(i) for i = IX1, ..., IX2. These * values must be nondecreasing. Repeated values are permitted * except for the first two and last two. * IX1, IX2 [in] Indices specifying the portion of the array XT() to * be considered by this subr. Require IX1 < IX2. * X [in] Value to be looked up. * LEFTI [inout] On entry LEFTI must contain an integer value. If it * in the range [IX1, IX2-1] the search will start with this index. * Otherwise the search will start with IX1 or IX2-1. * On return LEFTI will be in [IX1, IX2-1], and identifies * the interval (of nonzero length) from X(LEFTI) to * X(LEFTI+1) as the reference subinterval for X. * LEFTI = IX1 means X < XT(IX1+1) * IX1 < LEFTI < IX2-1 means XT(LEFTI) .le. X .lt. XT(LEFTI+1) * LEFTI = IX2-1 means XT(IX2-1) .le. X * MODE [out] Set to * -1 if X < XT(IX1) * 0 if XT(IX1) .le. X .le. XT(IX2) * +1 if XT(IX2) < X * ------------------------------------------------------------------ *--S replaces "?": ?SFIND, ?ERV1 * Both use IERM1 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ ilo = max( ix1, min( *lefti, ix2 - 1 ) ); if (x >= xt[ilo-(1)]) { if (x < xt[ilo + 1-(1)]) { *mode = 0; *lefti = ilo; return; } else { /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * procedure( SEARCH TO RIGHT ) */ step = 1; ihi = ilo + 1; /* do while(.true.) ! Using Fortran 90 "DO" syntax. */ L_5: ; if (ihi == ix2) { /* Here X .ge. XT(IX2) */ if (x == xt[ix2-(1)]) { *mode = 0; } else { *mode = 1; } ilo = ix2 - 1; *lefti = ilo; if (xt[ilo-(1)] >= xt[ix2-(1)]) { ierm1( "SSFIND", 2, 2, "Require T(IX2-1) < T(IX2)" , "IX2", ix2, ',' ); serv1( "T(IX2-1)", xt[ix2 - 1-(1)], ',' ); serv1( "T(IX2)", xt[ix2-(1)], '.' ); } return; } ilo = ihi; ihi = min( ihi + step, ix2 ); if (x < xt[ihi-(1)]) goto L_10; step *= 2.0e0; goto L_5; /* end do !while */ L_10: ; /* end proc !( SEARCH TO RIGHT ) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ } } else { /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * procedure( SEARCH TO LEFT ) */ step = 1; /* do while(.true.) */ L_15: ; if (ilo == ix1) { *mode = -1; *lefti = ix1; if (xt[ix1-(1)] >= xt[ix1 + 1-(1)]) { ierm1( "SSFIND", 1, 2, "Require T(IX1) < T(IX1+1)" , "IX1", ix1, ',' ); serv1( "T(IX1)", xt[ix1-(1)], ',' ); serv1( "T(IX1+1)", xt[ix1 + 1-(1)], '.' ); } return; } ihi = ilo; ilo = max( ilo - step, ix1 ); if (x >= xt[ilo-(1)]) goto L_20; step *= 2.0e0; goto L_15; /* end do !while */ L_20: ; /* end proc !( SEARCH TO LEFT ) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * procedure( BINARY SEARCH ) */ /* Here XT(ILO) .le. X .lt. XT(IHI) * * do while(.true.) */ L_25: ; middle = (ilo + ihi)/2; if (middle == ilo) goto L_30; if (x >= xt[middle-(1)]) { ilo = middle; } else { ihi = middle; } goto L_25; /* end do !while */ L_30: ; *mode = 0; *lefti = ilo; /* end proc !( BINARY SEARCH ) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ return; } /* end of function */