/*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 <math.h>
#include "fcrt.h"
#include "silup.h"
#include <float.h>
#include <stdlib.h>
#include <string.h>
/* 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 0<L3<NDEG+1.
* 4 (Argument K4) The absolute and relative errors expected in
* YT entries are specified in EOPT(K4) and EOPT(K4+1)
* respectively. The values provided here are used in
* estimating the error in the interpolation which is stored
* in EOPT(1).
* 5 (Argument K5, L5) Do the interpolation to the accuracy
* requested by the absolute error tolerance specified in
* EOPT(K5) and the relative error tolerance in EOPT(K5+1)
* respectively. An attempt is made to keep the final error
* < EOPT(K5) + EOPT(K5+1) * (abs(YT(i1) +abs(YT(i2)), where
* i1 and i2 are indices for table values close to X.
* Interpolation is done as for the default case, except that
* in this case NDEG gives the maximal degree polynomial to use
* in the interpolation, and polynomial interpolation of even
* degree can happen. (The form that gives the C1 continuity
* is not available in this case, and in fact continuity of the
* interpolant itself is not to be expected. The actual degree
* used in doing the interpolation is stored in the space for
* the argument L5. An error estimate is returned in EOPT(1).
* 6 (Argument K6) Do not use point k in the interpolation if
* YT(k) = EOPT(K6).
* 7 (Argument K7) YT(K7+i) gives the first derivative
* corresponding to the function value in YT(i). These
* derivatives are to be used in doing the interpolation. One
* gets a continuous interpolant only for NDEG = 3, 7, 11, and
* 15. The interpolating polynomial satisfies p(XT(i)) = YT(i),
* p'(XT(i)) = YT(K7+i) for values of i that give values of XT
* close to X. (Subject to keeping within one of the same
* number of XT's on either side of X, the maximum of |XT(i)-X|
* is minimized.) If NDEG is even a value of YT(i) is used
* without using the corresponding value of YT(K7+i).
* >253 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 */