/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:15 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dilup.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_cdilup { double badpt, dx[MAXDEG + 1-(0)+1], ebnd, ebndr; long int kaos, kextrp, kgoc, lderiv, lexerr, lexit, ndegec, idx[MAXDEG + 1-(0)+1]; LOGICAL32 geterr; } cdilup; /* end of COMMON translations */ void /*FUNCTION*/ dilup( double x, double *y, long ntab, double xt[], double yt[], long ndeg, long *lup, long iopt[], double 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; double 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]={"DILUP$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 double epsr = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Eopt = &eopt[0] - 1; double *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; double *const Xt = &xt[0] - 1; double *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 DILUP Krogh Dimension of IDAT increased to 4. *>> 2010-03-05 DILUP Krogh Added more checks for NDEG too big. *>> 2010-03-03 DILUP Krogh Fixed bug in flagging extrpolation. *>> 2006-05-11 DILUP Krogh No flag of extrapolation on equality. *>> 2006-04-01 DILUP Krogh Fixed minor bug just introduced. *>> 2006-03-31 DILUP Krogh Don't flag extrap. when on end point. *>> 2003-04-17 DILUP Krogh Fixed bugs in derivatives with vectors. *>> 2000-12-03 DILUP Krogh Fixed so no reference to YT when LEXIT=0. *>> 2000-12-01 DILUP Krogh Removed unused parameter MENTXT. *>> 1998-01-26 DILUP Krogh Extrap. due to too many bad pts. flagged. *>> 1997-09-30 DILUP Krogh Decremented NDEGI when LUP=4 & GETERR. *>> 1996-04-08 DILUP Krogh With LUP=4 & GETERR, NDEGI was 1 too small. *>> 1996-03-30 DILUP Krogh Added external statement. *>> 1995-12-01 DILUP Krogh Fixed bugs connected with option 6. *>> 1995-11-10 DILUP Krogh Fixed so char. data at col. 72 is not ' '. *>> 1995-03-03 DILUP Krogh Bug in last change made DILUPM fail. *>> 1994-12-22 DILUP Krogh Added Hermite interpolation. *>> 1994-11-11 DILUP Krogh Declared all vars. *>> 1994-10-20 DILUP Krogh Changes to use M77CON *>> 1994-09-12 DILUP Krogh Added CHGTYP code. *>> 1993-04-28 DILUP Krogh Additions for Conversion to C. *>> 1992-05-27 DILUP Krogh Fixed bug in error estimate. *>> 1992-04-08 DILUP Krogh Removed unused labels 510 and 2080. *>> 1991-10-17 DILUP Krogh Initial Code. * *--D 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 * DILUPM, 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 DILUP.) 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 DILUPM, 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 *********************** * * D1MACH Returns system parameters. * DMESS 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 CDILUP. 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. * D1MACH Function to get parameters of floating point arithmetic. * DMESS 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 (D1MACH(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 CDILUP. 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 CDILUP. Keeps track of state on variable order. * = 0 No variable order -- this value only set in DILUPM. * = 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 CDILUP. 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 CDILUP. 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 DILUP. 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 CDILUP. Loc. where first deriv. is stored in EOPT(). * LEXERR In common CDILUP. Loc. where absolute and relative error * information on Y is stored in EOPT. * LEXIT In common CDILUP. Defines action after finding location in XT. * 0 Don't compute anything for Y, just return. (Used for DILUPM.) * 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 DMESS 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 DMESS 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 CDILUP. 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 * DMESS. * */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA DILUP$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; cdilup.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 (cdilup.kaos <= 0) goto L_50; lndeg = 0; kgo = 3; cdilup.kaos = 1; ndegi = min( MAXDEG, Iopt[3] + 1 ); ndege = -ndegi; goto L_60; } cdilup.geterr = FALSE; cdilup.kextrp = -1; cdilup.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 (cdilup.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 (cdilup.kextrp < 0) { cdilup.kextrp = max( ndege - 1, min( ndege, 2 ) ); if (cdilup.kextrp < 0) cdilup.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 (ARITHIF(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 (fabs( Xt[ili - 1] - xi ) < fabs( 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 (fabs( Xt[ili + 1] - xi ) < fabs( Xt[ili] - xi )) { ili += 1; n += 1; } L_240: if (*lup <= 0) *lup = -ili; L_250: cdilup.dx[0] = xi - Xt[ili]; /* Get bounding indices and interpolate */ L_260: iui = ili; k = ili; l = -1; if (cdilup.lexit == 0) { ndegi = 0; if (cdilup.geterr) { if (kgo != 3) { if (kgo == 1) ndege += 1; cdilup.kextrp = min( cdilup.kextrp + 1, ndege ); } else { ndege = -ndege; } } } /* Just got index for next XT */ L_270: if (cdilup.lexit == 0) { l += 1; cdilup.idx[l] = k; goto L_350; } tp2 = Yt[k]; if (lbadpt) { /* Check if point should be discarded */ if (tp2 == cdilup.badpt) { if (iui != ili) { if (labs( kk ) == 1) { n -= 1; } else { n += 1; } } goto L_1020; } } L_280: l += 1; cdilup.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])/(cdilup.dx[i - 1] - cdilup.dx[l]); tp3 = (tp3 - tp2)/(cdilup.dx[i - 1] - cdilup.dx[l]); } pi[l] = pi[l - 1]*cdilup.dx[l - 1]; tp1 += pi[l]*tp2; } } else { cdilup.dx[l] = cdilup.dx[l - 1]; pi[l] = pi[l - 1]*cdilup.dx[l]; tp2 = tp3; tp1 += pi[l]*tp2; } } else if (l == 0) { tp1 = tp2; } else { pi[l] = pi[l - 1]*cdilup.dx[l - 1]; if (l <= ndegi) { /* Get divided differences & interpolate. */ for (i = 1; i <= l; i++) { tp2 = (tp2 - cy[i - 1])/(cdilup.dx[i - 1] - cdilup.dx[l]); } tp1 += pi[l]*tp2; } } cy[l] = tp2; L_350: if (l < ndege) goto L_1000; L_360: if (cdilup.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 = cdilup.idx[l]; tp2 = Yt[k]; } goto L_290; /* Got simple interpolated value. */ L_500: *y = tp1; if (!cdilup.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*(fabs( tp1 - *y ) + .03125e0*fabs( 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])/(cdilup.dx[i - 1] - cdilup.dx[l]); } tp2 = (tp2 - cy[l - 1])/(cdilup.dx[0] - cdilup.dx[1]); cy[l] = tp2; pi[l] = pi[l - 1]*cdilup.dx[0]; *y = tp1 + pi[l]*tp2; if (cdilup.geterr) errest = 1.5e0*fabs( pi[l - 1] )*fabs( ((cdilup.dx[l - 1]* (cdilup.dx[0] - cdilup.dx[1])/(cdilup.dx[l - 1] - cdilup.dx[l]) - cdilup.dx[0])*tp2) + fabs( .03125e0*cy[l - 1] ) ); goto L_2000; /* Variable order, check estimated error and convergence. */ L_900: ; if (cdilup.kaos >= 3) goto L_930; if (cdilup.kaos == 2) goto L_920; /* First time */ if (((kder != 0) && (cdilup.lexit > 1)) && (l == 0)) goto L_970; cdilup.kaos = 2; e2l = fabs( tp1 ); e2 = cdilup.ebnd + 1.e30; goto L_970; /* Second time */ L_920: cdilup.kaos = 4; ebndi = .66666e0*(cdilup.ebnd + cdilup.ebndr*(fabs( cy[0] ) + fabs( Yt[cdilup.idx[1]] ))); ef = cdilup.dx[0]; if (cdilup.lexit > 1) { ef = cdilup.dx[l] - cdilup.dx[0]; adjsav = ef; } e2 = fabs( ef*tp2 ); em = .75e0; goto L_950; /* Usual case */ L_930: e2l = e2; e1 = e2l*(5.e0*em/(double)( l )); ef *= cdilup.dx[l - 1]; e2 = fabs( ef*tp2 ); em = 0.5e0*em + e2/(e2l + e2 + 1.e-20); if (e2 >= e1) { if (cdilup.kaos == 3) { /* Apparently diverging, so quit. */ l -= 1; goto L_960; } else { cdilup.kaos = 3; } } else { cdilup.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 (cdilup.lexit > 1) { if (kder == 0) tp2 = 1.5e0*fabs( cdilup.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; cdilup.dx[l + 1] = xl; goto L_270; } cdilup.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; cdilup.dx[l + 1] = xu; goto L_270; } cdilup.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 + (double)( 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 = cdilup.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 - (double)( ili ))*h; xl = xu; cdilup.dx[0] = xu; n = ili + ili; if (tp1 > (double)( ili )) n += 1; if ((xi - Xt[1])*(Xt[1] + h*(double)( ntabi - 1 ) - xi) >= 0.e0) goto L_260; n = 6*(ili - 1); goto L_1250; /* Already set up */ L_1400: kgo = cdilup.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 = cdilup.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; cdilup.kextrp = min( Iopt[i], ndege ); goto L_10; /* Get derivatives of interpolant */ L_1600: ; i += 2; cdilup.lderiv = Iopt[i - 1]; cdilup.lexit = Iopt[i] + 1; goto L_10; /* Get expected errors in YT. */ L_1900: ; i += 1; cdilup.lexerr = Iopt[i]; /* Set to get error estimate. */ L_1930: ; cdilup.geterr = TRUE; goto L_10; /* Set for automatic order selection. */ L_1950: ; io += 2; lndeg = io; cdilup.ebnd = fmax( 0.e0, Eopt[Iopt[io - 1]] ); cdilup.ebndr = fmax( 0.e0, Eopt[Iopt[io - 1] + 1] ); kgo = 3; cdilup.kaos = 1; ndegi = min( MAXDEG, ndege + 1 ); ndege = -ndegi; goto L_1930; /* Set to ignore special points */ L_1970: ; lbadpt = TRUE; io += 1; cdilup.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 (cdilup.lexit < 2) goto L_2100; /* Compute derivatives of Y */ saved = (kgo%10) == 4; if (saved) { dxs = cdilup.dx[l - 1]; cdilup.dx[l - 1] = cdilup.dx[0]; } n = cdilup.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 + cdilup.lderiv - 1] = 0.e0; } n = l; } tp1 = cy[1]; pid[0] = 1.e0; for (i = 1; i <= (l - 1); i++) { pid[i] = pi[i] + cdilup.dx[i]*pid[i - 1]; tp1 += pid[i]*cy[i + 1]; } Eopt[cdilup.lderiv] = tp1; for (k = 2; k <= n; k++) { tp1 = cy[k]; for (i = 1; i <= (l - k); i++) { pid[i] += cdilup.dx[i + k - 1]*pid[i - 1]; tp1 += pid[i]*cy[i + k]; } Eopt[k + cdilup.lderiv - 1] = tp1; } if (saved) cdilup.dx[l - 1] = dxs; /* Save info. and return */ L_2100: ; if (cdilup.geterr) { if (epsr == 0.e0) { epsr = DBL_EPSILON; } tp1 = epsr; tp2 = 0.e0; if (cdilup.lexerr != 0) { tp2 = fmax( 0.e0, Eopt[cdilup.lexerr] ); tp1 = fmax( tp1, Eopt[cdilup.lexerr + 1] ); } if (iiflg == 2) errest *= 32.e0; Eopt[1] = errest + tp2 + tp1*(fabs( cy[0] ) + fabs( cdilup.dx[0]* cy[1] )); if ((kgo == 3) || (kgo == 13)) { if (ebndi != 0.e0) { if (Eopt[1] > cdilup.ebnd) { Errdat[1] = Eopt[1]; Errdat[2] = cdilup.ebnd; iiflg = -1; } } } } L_2110: cdilup.kgoc = -kgo; cdilup.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) cdilup.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]; dmess( 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 -- DILUP */ } /* end of function */