/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:19 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "silup.h" #include #include #include /* PARAMETER translations */ #define LTXTAB 8 #define LTXTAC 53 #define LTXTAD 92 #define LTXTAE 142 #define LTXTAF 171 #define LTXTAG 196 #define LTXTAH 250 #define LTXTAI 290 #define LTXTAJ 458 #define MAXDEG 15 #define MECONT 50 #define MEEMES 52 #define MEIVEC 57 #define MEMDA1 27 #define MERET 51 #define METEXT 53 /* end of PARAMETER translations */ /* COMMON translations */ struct t_csilup { float badpt, dx[MAXDEG + 1-(0)+1], ebnd, ebndr; long int kaos, kextrp, kgoc, lderiv, lexerr, lexit, ndegec, idx[MAXDEG + 1-(0)+1]; LOGICAL32 geterr; } csilup; /* end of COMMON translations */ void /*FUNCTION*/ silup( float x, float *y, long ntab, float xt[], float yt[], long ndeg, long *lup, long iopt[], float eopt[]) { LOGICAL32 indxed, lbadpt, linc, saved, ytord; long int _l0, i, idat[4], iiflg, ili, io, iui, k, kder, kgo, kk, l, lder, lndeg, lopt, n, ndege, ndegi, ndegq, ntabi; float adjsav, cy[MAXDEG + 1-(0)+1], dxs, e1, e2, e2l, ebndi, ef, em, errdat[2], errest, h, pi[MAXDEG + 1-(0)+1], pid[MAXDEG-(0)+1], tp1, tp2, tp3, xi, xl, xu, yl; static char mtxtaa[3][202]={"SILUP$BEstimated error = $F, Requested error = $F.$ENTAB = 1, no error estimate computed.$EToo many bad points, only $I point(s) available.$EXT(1) = XT(NTAB=$I) = $F ??$ELUP = 3, and XT(2) = 0.$EOrder ", "of requested derivative, $I, is > bound of $I.$ELUP > 3, with unprepared common block.$ENDEG = $I must be .le. $I or NTAB = $I must be at least $I. You must either decrease NDEG, increase NTAB, or if", " requesting an error estimate, turn off that request.$ENDEG = $I is not in the allowed interval of [0, $I].$ENTAB = $I is not in the allowed interval of [1, $I].$EIOPT($I) = $I is not a valid option.$E"}; static char mtxtab[1][15]={"IOPT(1:$M): $E"}; static long mact[5]={MEEMES,0,0,0,MECONT}; static long mactv[6]={MEMDA1,0,METEXT,MEIVEC,0,MERET}; static long mloc[9]={LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG, LTXTAH,LTXTAI,LTXTAJ}; static float epsr = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Eopt = &eopt[0] - 1; float *const Errdat = &errdat[0] - 1; long *const Idat = &idat[0] - 1; long *const Iopt = &iopt[0] - 1; long *const Mact = &mact[0] - 1; long *const Mactv = &mactv[0] - 1; long *const Mloc = &mloc[0] - 1; float *const Xt = &xt[0] - 1; float *const Yt = &yt[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2010-08-26 SILUP Krogh Dimension of IDAT increased to 4. *>> 2010-03-05 SILUP Krogh Added more checks for NDEG too big. *>> 2010-03-03 SILUP Krogh Fixed bug in flagging extrpolation. *>> 2006-05-11 SILUP Krogh No flag of extrapolation on equality. *>> 2006-04-01 SILUP Krogh Fixed minor bug just introduced. *>> 2006-03-31 SILUP Krogh Don't flag extrap. when on end point. *>> 2003-04-17 SILUP Krogh Fixed bugs in derivatives with vectors. *>> 2000-12-03 SILUP Krogh Fixed so no reference to YT when LEXIT=0. *>> 2000-12-01 SILUP Krogh Removed unused parameter MENTXT. *>> 1998-01-26 SILUP Krogh Extrap. due to too many bad pts. flagged. *>> 1997-09-30 SILUP Krogh Decremented NDEGI when LUP=4 & GETERR. *>> 1996-04-08 SILUP Krogh With LUP=4 & GETERR, NDEGI was 1 too small. *>> 1996-03-30 SILUP Krogh Added external statement. *>> 1995-12-01 SILUP Krogh Fixed bugs connected with option 6. *>> 1995-11-10 SILUP Krogh Fixed so char. data at col. 72 is not ' '. *>> 1995-03-03 SILUP Krogh Bug in last change made SILUPM fail. *>> 1994-12-22 SILUP Krogh Added Hermite interpolation. *>> 1994-11-11 SILUP Krogh Declared all vars. *>> 1994-10-20 SILUP Krogh Changes to use M77CON *>> 1994-09-12 SILUP Krogh Added CHGTYP code. *>> 1993-04-28 SILUP Krogh Additions for Conversion to C. *>> 1992-05-27 SILUP Krogh Fixed bug in error estimate. *>> 1992-04-08 SILUP Krogh Removed unused labels 510 and 2080. *>> 1991-10-17 SILUP Krogh Initial Code. * *--S replaces "?": ?ILUP, ?ILUPM, C?ILUP, ?MESS * * Polynomial Interpolation with look up. * Design/Code by Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA * * In addition to doing standard 1 dimensional interpolation, this * subroutine supports efficient interpolation of several functions * defined at the same values of the independent variable and supports * SILUPM, a subroutine for doing multidimensional interpolation. Error * estimates, Hermite interpolation, and derivatives of the interpolant * can be obtained via options. * Algorithms used are described in "Efficient Algorithms for Polynomial * Interpolation and Numerical Differentiation", by Fred T. Krogh, Math. * of Comp. Vol. 24, #109 (Jan. 1970), pp. 185-190. * * ************* Formal Arguments *************************** * * X Independent variable where value of interpolant is desired. * Y Value of interpolant, computed by this subroutine. The * interpolant is always a piecewise polynomial. * NTAB Number of points in the table. * XT Array of independent variable values. Must be monotone * increasing or monotone decreasing. If XT(I) = XT(I+1) for some * I, then the XT(J)'s that are used in the interpolation will * either have all J's .le. I, or all J's .ge. I+1. If the XT's * are equally spaced an option allows one to provide only XT(1), * and the increment between the XT's. * YT Array of dependent variable values. Y(XT(I)) = YT(I). * NDEG If NDEG < 2 or odd, then it gives the degree of the polynomial * used. Else the polynomial used is a linear combination of two * polynomials of degree NDEG, such that the resulting polynomial * is of degree NDEG+1 and has a continuous first derivative. * Typically accuracy will improve as NDEG increases up to some * point where it will start to get worse because of either * rounding errors or the inherant instability of high degree * polynoimial interpolation. If more than MAXDEG-th degree is * desired, parameter MAXDEG must be changed below. It is * currently 15. If X is so close to the end of the table that * the same number of points can not be selected on both sides of * x the degree of the interpolant is NDEG exactly. When * extrapolating, the degree used is max(2, 2*(NDEG/2)), where the * divide is truncated to the nearest integer. * LUP Defines the type of look up method. (Changed if LUP < 1.) * < 0 Use a sequential search starting with an index of -LUP, and * set LUP to -k on exit where k minimizes abs(X-XT(k)). * = 0 As for < 0, except start with a binary search. * = 1 Use a binary search. * = 2 Start a sequential search with an index = [1.5 + * (NTAB-1) * (X-XT(1)) / (XT(NTAB)-XT(1))]. * = 3 YT(k) corresponds to XT(1) + (k-1)*XT(2), no search is needed * to do the look up and XT can have dimension 2. * = 4 Internal information connected with X-XT(k) values used in * in the last interpolation is reused. (Only use if there are * no intervening calls to SILUP.) No options should be * specified and only YT should be different. Intended for * interpolating components after the first of a vector valued * function. * IOPT IOPT(1) is used to return a status as follows. * -10 An option index is out of range. * -9 NTAB is outside of allowed limits. * -8 NDEG is outside of allowed limits. * -7 LUP > 3 (use old parameters), used when not ready for it. * -6 Option 3 (compute derivatives), has requested more than * MAXDEG derivatives. * -5 LUP = 3 and XT(2) = 0. * -4 XT(1) = XT(NTAB), and NTAB is not 1. * -3 < 2 points available because of bad points. * -2 Only one table entry available, req. err. est. not computed. * -1 The accuracy requested was not obtained. * 0 Nothing special to flag. * 1 X was outside the domain of the table, extrapolation used. * 2 NTAB is so small, it restricted the degree of the polynomial. */ /* Starting with IOPT(2) options are specified by integers in the * range 0-6, followed in some cases by integers providing * argument(s). Each option, together with its options its * arguments if any is followed in IOPT by the next option (or 0). * 0 No more options; this must be last in the option list. * 1 An error estimate is to be returned in EOPT(1). * 2 (Argument: K2) K2 gives the polynomial degree to use when * extrapolating. * 3 (Argument: K3, L3) Save (k-th derivative of interpolating * polynomial) / (k!) in EOPT(K3+K-1) for k = 1, 2, ..., L3. * These values are the coefficients of the polynomial in the * monomial basis expanded about X. One must have 0253 Used for calls from SILUPM, a number of variables in the * common block have been set. If this is used, it must be the * first option, and common variables must be set. The rest * of IOPT is not examined for options. * EOPT Array used to return an error estimate and also used for * options. * EOPT(1) if an error estimate is returned, this contains a crude * estimate of the error in the interpolation. * * ************ External Procedures *********************** * * R1MACH Returns system parameters. * SMESS Prints error messages. * * ************ Variables referenced *********************** * * ADJSAV Saved difference of XT values used to adjust error estimate on * variable order when derivatives are being computed. * BADPT In common CSILUP. YT(K) is not to be used if YT(K) = BADPT, * and LBADPT = .true. * CY Array giving the YT values used to construct the interpolant. * Also used to store divided differences. * R1MACH Function to get parameters of floating point arithmetic. * SMESS Program calling MESS to print error messages. * DX Array giving the differences X - XT(I), where I runs through * the points that are used in the interpolation. * DXS Saved value of DX when computing derivatives * E1 Error estimate from one iteration ago (for variable order) * E2 Error estimate from this iteration (for variable order) * E2L Used in computing error estimate for variable order. * EBND If estimated error is < EBND then order is suff. high. * EBNDI Internal value for bounding error. * EBNDR Used in getting part of EBND due to relative accuracy request. * EF Factor used in estimating errors for variable order. * EN Used in computing error estimate for variable order. * EOPT Formal argument, see above. * EPSR Relative error level (R1MACH(4)). * ERRDAT Holds floating point data for error messages. * ERREST Estimate of the error in the interpolation. * GETERR Logical variable that is .true. if we are getting an error est. * H In indexed lookup contains the difference between XT values. * I Temporary index. * IO Index used in processing options. * IDX In common CSILUP. Gives indices of points selected for the * interpolation in order selected. * IIFLG Internal value for IOPT(1). * ILI Lower (sometimes upper) bound on the XT indices that will be * used in the interpolation. * INDXED Logical variable that is .true. if we have an indexed XT. * IOPT Formal argument, see above. * IUI As for ILI, except an upper (sometimes lower) bound. * K Index usually into YT and XT, into CY for derivatives. * KAOS In common CSILUP. Keeps track of state on variable order. * = 0 No variable order -- this value only set in SILUPM. * = 1 First time * = 2 Second time */ /* = 3 Just had an increase or a very weak decrease in the error. * = 4 Just had a strong decrease in the error. * On an error with IIFLG .lt. 0, this is set to -10 * stop level - * the print level. * KDER 0 if not doing Hermite interpolation, else is < 0 if have next * higher order difference computed (in TP3), > 0 if it isn't. * KEXTRP In common CSILUP. Gives degree to use when extrapolating. * KGO Indicates type of interpolation as follows: * 1 Standard polynomial interpolation. * 2 Get error estimate after 1. * 3 Compute an interpolant with some desired accuracy. * 4 Compute an interpolant with a continuous derivative. * >10 Same as if KGO were 10 smaller, except XT values have already * been selected. * KGOC In common CSILUP. Value saved for -KGO. If YT contains values * of Y computed elsewhere in the order specified in IDX, then * KGOC should be set to abs(KGOC) before calling SILUP. If -2, * an extra value was computed for an error estimate. * KK Used in selecting data points to use. * 1 decrease ILI only * 0 increase IUI * -1 decrease ILI * <-1, increase IUI only * L Temporary index used in searching the XT array. Also one less * than the number of points used in the current interpolant. * LBADPT Logical variable set = .true. if checking for bad points. * LDER Offset of y' values from the y values. * LDERIV In common CSILUP. Loc. where first deriv. is stored in EOPT(). * LEXERR In common CSILUP. Loc. where absolute and relative error * information on Y is stored in EOPT. * LEXIT In common CSILUP. Defines action after finding location in XT. * 0 Don't compute anything for Y, just return. (Used for SILUPM.) * 1 Usual case, just interpolate. * >1 Compute LEXIT-1 derivatives of the interpolant. * LINC Logical variable used in the sequential search. Set .true. if * the values in XT are increasing with I, and set .false. otherwise. * LNDEG Location in IOPT to save degree when variable order is used. * LOPT Index of the last option. * LUP Formal argument, see above. * MACT Array used for error message actions, see SMESS for details. * MACTV Array used for part of error message for IOPT. * MAXDEG Parameter giving the maximum degree polynomial interpolation * supported. MAXDEG must be odd and > 2. * MESS Message program called from SMESS to print error messages. * MEEMES Parameter giving value for printing an error message in MESS. * MEIVEC Parameter giving value for printing an integer vector in MESS. * MEMDA1 Parameter giving value for indicating value in MACT for MESS. * MEMDA2 Parameter giving value for indicating value in MACT for MESS. * MEMDAT Parameter giving value for specifying MACT index for MESS. * MERET Parameter giving value to indicate no more actions for MESS. * METEXT Parameter giving value to tell MESS to print from MTXTAA. * LTXTxx Parameter names of this form were generated by PMESS in making * up text for error messages and are used to locate various parts * of the message. * MLOC Array giving starting loctions for error message text. * MTXTAA Character array holding error message text for MESS. * N Used in logic for deciding bounds. If N = 0, ILI can not be * reduced further. If N = 3*NTABI, IUI can not be increased further. * NDEG Formal argument, see above. * NDEGE Degree of polynomial to be used. Starts out = NDEG. * NDEGEC In common CSILUP. Degree actually used, saved in common. * NDEGI Degree of polynomial up to which y & differences are computed. * NDEGQ Used when KGO > 10. In this case If L .ge. NDEGQ special * action is needed. (Usually means quitting.) * NTAB Formal argument, see above. * NTABI Internal value of NTAB, = number of points in XT and YT. * PI Array use to store interpolation coefficients * PID Temporary storage used when computing derivatives. * SAVED Logical variable set = .true. if a DX value needs to be * replaced and restored when computing derivatives. * TP1 Used for temporary accumulation of values and temp. storage. * TP2 Used for temporary storage. * TP3 Used for divided difference when doing Hermite interpolation. * X Formal argument, see above. * XI Internal value for X, the place where interpolation is desired. * XL Value of "left hand" X when selecting indexed points. * XT Formal argument, see above. * XU Value of "right hand" X when selecting indexed points. * Y Formal argument, see above. * YL Last value of Y when selecting the order. * YT Formal argument, see above. * YTORD Logical variable set .true. when YT values are ordered on entry * * ************* Formal Variable Declarations *************** * */ /* ************* Common Block and Parameter ***************** * * MAXDEG must be odd and > 2. */ /* ************* Local Variables **************************** * */ /* ************************ Error Message Stuff and Data **************** * * Parameter defined below are all defined in the error message program * SMESS. * */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA SILUP$B *AB Estimated error = $F, Requested error = $F.$E *AC NTAB = 1, no error estimate computed.$E *AD Too many bad points, only $I point(s) available.$E *AE XT(1) = XT(NTAB=$I) = $F ??$E *AF LUP = 3, and XT(2) = 0.$E *AG Order of requested derivative, $I, is > bound of $I.$E *AH LUP > 3, with unprepared common block.$E *AI NDEG = $I must be .le. $I or NTAB = $I must be at $C * least $I. You must either decrease NDEG, increase NTAB, $C * or if requesting an error estimate, turn off that request.$E *AJ NDEG = $I is not in the allowed interval of [0, $I].$E *AK NTAB = $I is not in the allowed interval of [1, $I].$E *AL IOPT($I) = $I is not a valid option.$E * $ *AM IOPT(1:$M): $E */ /* **** End of automatically generated text * * 1 2 3 4 5 */ /* 2 3 4 5 6 */ /* ************ Start of executable code ******************* * */ iiflg = 0; kder = 0; pi[0] = 1.e0; lbadpt = FALSE; ili = -*lup; if (ili <= -4) goto L_1400; csilup.lexerr = 0; ntabi = ntab; if ((ntabi <= 0) || (ntabi > 9999999)) { iiflg = -9; Idat[1] = ntabi; Idat[2] = 9999999; goto L_2120; } ndegi = ndeg; ndege = ndegi; xi = x; kgo = 1; io = 1; if (Iopt[2] >= 254) { lbadpt = Iopt[2] == 255; if (csilup.kaos <= 0) goto L_50; lndeg = 0; kgo = 3; csilup.kaos = 1; ndegi = min( MAXDEG, Iopt[3] + 1 ); ndege = -ndegi; goto L_60; } csilup.geterr = FALSE; csilup.kextrp = -1; csilup.lexit = 1; /* Loop to take care of options */ L_10: io += 1; lopt = Iopt[io]; if (lopt != 0) { switch (lopt) { case 1: goto L_1930; case 2: goto L_1500; case 3: goto L_1600; case 4: goto L_1900; case 5: goto L_1950; case 6: goto L_1970; case 7: goto L_1980; } iiflg = -10; Idat[1] = io; Idat[2] = lopt; goto L_2120; } L_50: ; k = ndegi + 1; if (kder != 0) { k /= 2; } else if ((ndegi%2) == 0) { k += 1; } if (csilup.geterr) k += 1; if (((ndegi < 0) || (ndegi > MAXDEG)) || (k > ntab)) { iiflg = -8; Idat[1] = ndegi; Idat[2] = MAXDEG; Idat[3] = ntab; Idat[4] = k; goto L_2120; } if (kgo == 1) { if (ndegi >= 2) { if ((ndegi%2) == 0) { if (kder == 0) { ndege += 1; kgo = 4; } } } } L_60: if (csilup.kextrp < 0) { csilup.kextrp = max( ndege - 1, min( ndege, 2 ) ); if (csilup.kextrp < 0) csilup.kextrp = ndegi; } if (ili > 0) goto L_200; switch (IARITHIF(ili + 2)) { case -1: goto L_1300; case 0: goto L_1200; case 1: goto L_100; } /* Binary search, then sequential */ L_100: ; /* In binary search XT(ILI) .le. XI .le. XT(IUI) * or we are extrapolating */ ili = 1; iui = ntabi; switch (SARITHIF(Xt[ntabi] - Xt[1])) { case -1: goto L_110; case 0: goto L_1220; case 1: goto L_130; } L_110: ili = ntabi; l = 1; L_120: iui = l; L_130: l = (iui - ili)/2; if (l == 0) goto L_210; l += ili; if (Xt[l] > xi) goto L_120; ili = l; goto L_130; /* Sequential search, then exit */ L_200: ; if (ili > ntabi) goto L_100; if (Xt[ntabi] == Xt[1]) goto L_1220; L_210: indxed = FALSE; linc = Xt[ntabi] > Xt[1]; if ((Xt[ili] > xi) == linc) goto L_230; L_220: if (ili == ntabi) goto L_1230; ili += 1; if ((Xt[ili] < xi) == linc) goto L_220; n = 2*ili; if (fabsf( Xt[ili - 1] - xi ) < fabsf( Xt[ili] - xi )) { ili -= 1; n -= 1; } goto L_240; L_230: if (ili == 1) goto L_1230; ili -= 1; if ((Xt[ili] > xi) == linc) goto L_230; n = 2*ili + 1; if (fabsf( Xt[ili + 1] - xi ) < fabsf( Xt[ili] - xi )) { ili += 1; n += 1; } L_240: if (*lup <= 0) *lup = -ili; L_250: csilup.dx[0] = xi - Xt[ili]; /* Get bounding indices and interpolate */ L_260: iui = ili; k = ili; l = -1; if (csilup.lexit == 0) { ndegi = 0; if (csilup.geterr) { if (kgo != 3) { if (kgo == 1) ndege += 1; csilup.kextrp = min( csilup.kextrp + 1, ndege ); } else { ndege = -ndege; } } } /* Just got index for next XT */ L_270: if (csilup.lexit == 0) { l += 1; csilup.idx[l] = k; goto L_350; } tp2 = Yt[k]; if (lbadpt) { /* Check if point should be discarded */ if (tp2 == csilup.badpt) { if (iui != ili) { if (labs( kk ) == 1) { n -= 1; } else { n += 1; } } goto L_1020; } } L_280: l += 1; csilup.idx[l] = k; L_290: if (kder != 0) { /* Hermite interpolation */ kder = -kder; if (kder < 0) { tp3 = Yt[lder + k]; if (l == 0) { tp1 = tp2; } else { for (i = 1; i <= l; i++) { tp2 = (tp2 - cy[i - 1])/(csilup.dx[i - 1] - csilup.dx[l]); tp3 = (tp3 - tp2)/(csilup.dx[i - 1] - csilup.dx[l]); } pi[l] = pi[l - 1]*csilup.dx[l - 1]; tp1 += pi[l]*tp2; } } else { csilup.dx[l] = csilup.dx[l - 1]; pi[l] = pi[l - 1]*csilup.dx[l]; tp2 = tp3; tp1 += pi[l]*tp2; } } else if (l == 0) { tp1 = tp2; } else { pi[l] = pi[l - 1]*csilup.dx[l - 1]; if (l <= ndegi) { /* Get divided differences & interpolate. */ for (i = 1; i <= l; i++) { tp2 = (tp2 - cy[i - 1])/(csilup.dx[i - 1] - csilup.dx[l]); } tp1 += pi[l]*tp2; } } cy[l] = tp2; L_350: if (l < ndege) goto L_1000; L_360: if (csilup.lexit == 0) goto L_2110; switch (kgo) { case 1: goto L_500; case 2: goto L_600; case 3: goto L_900; case 4: goto L_800; } if (l >= ndegq) switch (kgo - 10) { case 1: goto L_500; case 2: goto L_600; case 3: goto L_900; case 4: goto L_800; } /* Already got the points selected. */ L_400: l += 1; if (ytord) { tp2 = Yt[l + 1]; } else { k = csilup.idx[l]; tp2 = Yt[k]; } goto L_290; /* Got simple interpolated value. */ L_500: *y = tp1; if (!csilup.geterr) goto L_2000; ndegi += 1; kgo += 1; if (ndege >= 0) goto L_1000; goto L_400; /* Got info. for error estimate. */ L_600: if (l == 0) { iiflg = -2; } else { errest = 1.5e0*(fabsf( tp1 - *y ) + .03125e0*fabsf( pi[l - 1]* cy[l - 1] )); l -= 1; } goto L_2000; /* C1 interpolant. */ L_800: for (i = 1; i <= (l - 1); i++) { tp2 = (tp2 - cy[i - 1])/(csilup.dx[i - 1] - csilup.dx[l]); } tp2 = (tp2 - cy[l - 1])/(csilup.dx[0] - csilup.dx[1]); cy[l] = tp2; pi[l] = pi[l - 1]*csilup.dx[0]; *y = tp1 + pi[l]*tp2; if (csilup.geterr) errest = 1.5e0*fabsf( pi[l - 1] )*fabsf( ((csilup.dx[l - 1]* (csilup.dx[0] - csilup.dx[1])/(csilup.dx[l - 1] - csilup.dx[l]) - csilup.dx[0])*tp2) + fabsf( .03125e0*cy[l - 1] ) ); goto L_2000; /* Variable order, check estimated error and convergence. */ L_900: ; if (csilup.kaos >= 3) goto L_930; if (csilup.kaos == 2) goto L_920; /* First time */ if (((kder != 0) && (csilup.lexit > 1)) && (l == 0)) goto L_970; csilup.kaos = 2; e2l = fabsf( tp1 ); e2 = csilup.ebnd + 1.e30; goto L_970; /* Second time */ L_920: csilup.kaos = 4; ebndi = .66666e0*(csilup.ebnd + csilup.ebndr*(fabsf( cy[0] ) + fabsf( Yt[csilup.idx[1]] ))); ef = csilup.dx[0]; if (csilup.lexit > 1) { ef = csilup.dx[l] - csilup.dx[0]; adjsav = ef; } e2 = fabsf( ef*tp2 ); em = .75e0; goto L_950; /* Usual case */ L_930: e2l = e2; e1 = e2l*(5.e0*em/(float)( l )); ef *= csilup.dx[l - 1]; e2 = fabsf( ef*tp2 ); em = 0.5e0*em + e2/(e2l + e2 + 1.e-20); if (e2 >= e1) { if (csilup.kaos == 3) { /* Apparently diverging, so quit. */ l -= 1; goto L_960; } else { csilup.kaos = 3; } } else { csilup.kaos = 4; } L_950: yl = tp1; if (e2l + e2 > ebndi) { if ((l + ndege < 0) && (iiflg != 2)) goto L_970; if (kgo == 13) { /* May need early exit to get next point. */ Iopt[2] = 0; if (l < ndeg) goto L_2110; } } L_960: tp2 = 1.5e0; if (csilup.lexit > 1) { if (kder == 0) tp2 = 1.5e0*fabsf( csilup.dx[0]/adjsav ); } errest = tp2*(e2 + .0625e0*e2l); if (lndeg != 0) Iopt[lndeg] = l; *y = yl; goto L_2000; L_970: if (kgo > 10) goto L_400; L_1000: if (kder < 0) goto L_280; L_1020: kk = min( n - iui - ili, 2 ) - 1; /* In this section of code, KK=: 1, decrease ILI only; 0, increase IUI; * -1, decrease ILI; and <-1, increase IUI only */ if (labs( kk ) == 1) { ili -= 1; k = ili; if (ili != 0) { if (indxed) { xl += h; csilup.dx[l + 1] = xl; goto L_270; } csilup.dx[l + 1] = xi - Xt[k]; if (Xt[ili + 1] != Xt[ili]) goto L_270; } if (kk == 1) goto L_1100; n = 0; /* If all points will be on one side, flag extrapolation. * if (L .le. 0) go to 1250 */ } else { iui += 1; k = iui; if (iui <= ntabi) { if (indxed) { xu -= h; csilup.dx[l + 1] = xu; goto L_270; } csilup.dx[l + 1] = xi - Xt[k]; if (Xt[iui - 1] != Xt[iui]) goto L_270; } if (kk != 0) goto L_1100; n = 3*ntabi; /* If all points will be on one side, flag extarpolation. * if (L .le. 0) go to 1250 */ } if (kgo < 4) goto L_1020; /* If we can't get the same number of points on either side, the * continuous derivative case becomes standard interpolation. */ kgo = 1; ndege -= 1; if (l < ndege) goto L_1020; goto L_360; /* No more data accessible. */ L_1100: if (l <= 0) { /* Too many bad points, couldn't find points. */ *y = tp1; Idat[1] = l + 1; iiflg = -3; goto L_2110; } ndegi = min( ndegi, l ); ndege = 0; iiflg = 2; goto L_360; /* Secant start, then use sequential */ L_1200: ; if (Xt[1] == Xt[ntabi]) goto L_1220; ili = max( 1, min( ntabi, (long)( 1.5e0 + (float)( ntabi - 1 )* (xi - Xt[1])/(Xt[ntabi] - Xt[1]) ) ) ); goto L_210; /* Special cases * 1 entry in XT */ L_1220: ili = 1; iui = 1; kgo = 1; k = 1; if (ndege != 0) iiflg = 2; ndege = 0; if (ntabi == 1) goto L_250; /* Error -- XT(1) .eq. XT(NTAB), and NTAB .ne. 1 */ iiflg = -4; Idat[1] = ntabi; Errdat[1] = Xt[1]; kgo = 0; goto L_2120; /* Check if really extrapolating */ L_1230: n = 6*(ili - 1); if ((xi - Xt[1])*(Xt[ntabi] - xi) >= 0.e0) goto L_240; /* Extrapolating */ L_1250: iiflg = 1; if (kgo >= 4) kgo = 1; ndegi = csilup.kextrp; ndege = isign( ndegi, ndege ); if (indxed) goto L_260; goto L_240; /* Index search, then exit */ L_1300: ; indxed = TRUE; h = Xt[2]; if (h == 0) { iiflg = -5; goto L_2120; } tp1 = 1.e0 + (xi - Xt[1])/h; ili = min( max( 1, nint( tp1 ) ), ntabi ); xu = (tp1 - (float)( ili ))*h; xl = xu; csilup.dx[0] = xu; n = ili + ili; if (tp1 > (float)( ili )) n += 1; if ((xi - Xt[1])*(Xt[1] + h*(float)( ntabi - 1 ) - xi) >= 0.e0) goto L_260; n = 6*(ili - 1); goto L_1250; /* Already set up */ L_1400: kgo = csilup.kgoc; if (kgo == 0) { iiflg = -7; goto L_2120; } ytord = kgo > 0; kgo = labs( kgo ); if (kgo < 10) kgo += 10; if (kgo == 12) kgo = 11; ndegi = csilup.ndegec; ndegq = ndegi; if (kgo == 13) ndegq = 0; ndege = -ndegi; l = -1; if (kgo == 14) ndegi -= 1; goto L_400; /* Set number of points for extrapolation */ L_1500: ; i += 1; csilup.kextrp = min( Iopt[i], ndege ); goto L_10; /* Get derivatives of interpolant */ L_1600: ; i += 2; csilup.lderiv = Iopt[i - 1]; csilup.lexit = Iopt[i] + 1; goto L_10; /* Get expected errors in YT. */ L_1900: ; i += 1; csilup.lexerr = Iopt[i]; /* Set to get error estimate. */ L_1930: ; csilup.geterr = TRUE; goto L_10; /* Set for automatic order selection. */ L_1950: ; io += 2; lndeg = io; csilup.ebnd = fmaxf( 0.e0, Eopt[Iopt[io - 1]] ); csilup.ebndr = fmaxf( 0.e0, Eopt[Iopt[io - 1] + 1] ); kgo = 3; csilup.kaos = 1; ndegi = min( MAXDEG, ndege + 1 ); ndege = -ndegi; goto L_1930; /* Set to ignore special points */ L_1970: ; lbadpt = TRUE; io += 1; csilup.badpt = Eopt[Iopt[io]]; goto L_10; /* Set up for Hermite interpolation. */ L_1980: ; io += 1; lder = Iopt[io]; kder = labs( lder ); goto L_10; /* Put any new options just above here * End of the loop * */ L_2000: if (csilup.lexit < 2) goto L_2100; /* Compute derivatives of Y */ saved = (kgo%10) == 4; if (saved) { dxs = csilup.dx[l - 1]; csilup.dx[l - 1] = csilup.dx[0]; } n = csilup.lexit - 1; if (n > l) { if (n > MAXDEG) { Idat[1] = n; Idat[2] = MAXDEG; n = MAXDEG; iiflg = -6; } for (i = l + 1; i <= n; i++) { Eopt[i + csilup.lderiv - 1] = 0.e0; } n = l; } tp1 = cy[1]; pid[0] = 1.e0; for (i = 1; i <= (l - 1); i++) { pid[i] = pi[i] + csilup.dx[i]*pid[i - 1]; tp1 += pid[i]*cy[i + 1]; } Eopt[csilup.lderiv] = tp1; for (k = 2; k <= n; k++) { tp1 = cy[k]; for (i = 1; i <= (l - k); i++) { pid[i] += csilup.dx[i + k - 1]*pid[i - 1]; tp1 += pid[i]*cy[i + k]; } Eopt[k + csilup.lderiv - 1] = tp1; } if (saved) csilup.dx[l - 1] = dxs; /* Save info. and return */ L_2100: ; if (csilup.geterr) { if (epsr == 0.e0) { epsr = FLT_EPSILON; } tp1 = epsr; tp2 = 0.e0; if (csilup.lexerr != 0) { tp2 = fmaxf( 0.e0, Eopt[csilup.lexerr] ); tp1 = fmaxf( tp1, Eopt[csilup.lexerr + 1] ); } if (iiflg == 2) errest *= 32.e0; Eopt[1] = errest + tp2 + tp1*(fabsf( cy[0] ) + fabsf( csilup.dx[0]* cy[1] )); if ((kgo == 3) || (kgo == 13)) { if (ebndi != 0.e0) { if (Eopt[1] > csilup.ebnd) { Errdat[1] = Eopt[1]; Errdat[2] = csilup.ebnd; iiflg = -1; } } } } L_2110: csilup.kgoc = -kgo; csilup.ndegec = l; L_2120: Iopt[1] = iiflg; if (iiflg >= 0) return; /* Take care of error messages, IIFLG < -2 should stop. */ Mact[2] = 88; if (iiflg < -1) csilup.kaos = -Mact[2]; if (iiflg >= -3) Mact[2] = min( 26, 24 - iiflg ); if (Iopt[2] >= 254) { if (iiflg == -1) return; Mact[2] = 28; } Mact[3] = -iiflg; Mact[4] = Mloc[-iiflg]; smess( mact, (char*)mtxtaa,202, idat, errdat ); Mactv[2] = io; Mactv[5] = io; k = 1; if (io == 1) k = 6; mess( &Mactv[k], (char*)mtxtab,15, iopt ); return; /* End of Interpolation routine -- SILUP */ } /* end of function */