1 subroutine DILUP (X, Y, NTAB, XT, YT, NDEG, LUP, IOPT, EOPT) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 2010-08-26 DILUP Krogh Dimension of IDAT increased to 4. 6c>> 2010-03-05 DILUP Krogh Added more checks for NDEG too big. 7c>> 2010-03-03 DILUP Krogh Fixed bug in flagging extrpolation. 8c>> 2006-05-11 DILUP Krogh No flag of extrapolation on equality. 9c>> 2006-04-01 DILUP Krogh Fixed minor bug just introduced. 10c>> 2006-03-31 DILUP Krogh Don't flag extrap. when on end point. 11c>> 2003-04-17 DILUP Krogh Fixed bugs in derivatives with vectors. 12c>> 2000-12-03 DILUP Krogh Fixed so no reference to YT when LEXIT=0. 13c>> 2000-12-01 DILUP Krogh Removed unused parameter MENTXT. 14c>> 1998-01-26 DILUP Krogh Extrap. due to too many bad pts. flagged. 15c>> 1997-09-30 DILUP Krogh Decremented NDEGI when LUP=4 & GETERR. 16c>> 1996-04-08 DILUP Krogh With LUP=4 & GETERR, NDEGI was 1 too small. 17c>> 1996-03-30 DILUP Krogh Added external statement. 18C>> 1995-12-01 DILUP Krogh Fixed bugs connected with option 6. 19C>> 1995-11-10 DILUP Krogh Fixed so char. data at col. 72 is not ' '. 20C>> 1995-03-03 DILUP Krogh Bug in last change made DILUPM fail. 21C>> 1994-12-22 DILUP Krogh Added Hermite interpolation. 22C>> 1994-11-11 DILUP Krogh Declared all vars. 23c>> 1994-10-20 DILUP Krogh Changes to use M77CON 24c>> 1994-09-12 DILUP Krogh Added CHGTYP code. 25c>> 1993-04-28 DILUP Krogh Additions for Conversion to C. 26c>> 1992-05-27 DILUP Krogh Fixed bug in error estimate. 27c>> 1992-04-08 DILUP Krogh Removed unused labels 510 and 2080. 28c>> 1991-10-17 DILUP Krogh Initial Code. 29c 30c--D replaces "?": ?ILUP, ?ILUPM, C?ILUP, ?MESS 31c 32c Polynomial Interpolation with look up. 33c Design/Code by Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA 34c 35c In addition to doing standard 1 dimensional interpolation, this 36c subroutine supports efficient interpolation of several functions 37c defined at the same values of the independent variable and supports 38c DILUPM, a subroutine for doing multidimensional interpolation. Error 39c estimates, Hermite interpolation, and derivatives of the interpolant 40c can be obtained via options. 41c Algorithms used are described in "Efficient Algorithms for Polynomial 42c Interpolation and Numerical Differentiation", by Fred T. Krogh, Math. 43c of Comp. Vol. 24, #109 (Jan. 1970), pp. 185-190. 44c 45c ************* Formal Arguments *************************** 46c 47c X Independent variable where value of interpolant is desired. 48c Y Value of interpolant, computed by this subroutine. The 49c interpolant is always a piecewise polynomial. 50c NTAB Number of points in the table. 51c XT Array of independent variable values. Must be monotone 52c increasing or monotone decreasing. If XT(I) = XT(I+1) for some 53c I, then the XT(J)'s that are used in the interpolation will 54c either have all J's .le. I, or all J's .ge. I+1. If the XT's 55c are equally spaced an option allows one to provide only XT(1), 56c and the increment between the XT's. 57c YT Array of dependent variable values. Y(XT(I)) = YT(I). 58c NDEG If NDEG < 2 or odd, then it gives the degree of the polynomial 59c used. Else the polynomial used is a linear combination of two 60c polynomials of degree NDEG, such that the resulting polynomial 61c is of degree NDEG+1 and has a continuous first derivative. 62c Typically accuracy will improve as NDEG increases up to some 63c point where it will start to get worse because of either 64c rounding errors or the inherant instability of high degree 65c polynoimial interpolation. If more than MAXDEG-th degree is 66c desired, parameter MAXDEG must be changed below. It is 67c currently 15. If X is so close to the end of the table that 68c the same number of points can not be selected on both sides of 69c x the degree of the interpolant is NDEG exactly. When 70c extrapolating, the degree used is max(2, 2*(NDEG/2)), where the 71c divide is truncated to the nearest integer. 72c LUP Defines the type of look up method. (Changed if LUP < 1.) 73c < 0 Use a sequential search starting with an index of -LUP, and 74c set LUP to -k on exit where k minimizes abs(X-XT(k)). 75c = 0 As for < 0, except start with a binary search. 76c = 1 Use a binary search. 77c = 2 Start a sequential search with an index = [1.5 + 78c (NTAB-1) * (X-XT(1)) / (XT(NTAB)-XT(1))]. 79c = 3 YT(k) corresponds to XT(1) + (k-1)*XT(2), no search is needed 80c to do the look up and XT can have dimension 2. 81c = 4 Internal information connected with X-XT(k) values used in 82c in the last interpolation is reused. (Only use if there are 83c no intervening calls to DILUP.) No options should be 84c specified and only YT should be different. Intended for 85c interpolating components after the first of a vector valued 86c function. 87c IOPT IOPT(1) is used to return a status as follows. 88c -10 An option index is out of range. 89c -9 NTAB is outside of allowed limits. 90c -8 NDEG is outside of allowed limits. 91c -7 LUP > 3 (use old parameters), used when not ready for it. 92c -6 Option 3 (compute derivatives), has requested more than 93c MAXDEG derivatives. 94c -5 LUP = 3 and XT(2) = 0. 95c -4 XT(1) = XT(NTAB), and NTAB is not 1. 96c -3 < 2 points available because of bad points. 97c -2 Only one table entry available, req. err. est. not computed. 98c -1 The accuracy requested was not obtained. 99c 0 Nothing special to flag. 100c 1 X was outside the domain of the table, extrapolation used. 101c 2 NTAB is so small, it restricted the degree of the polynomial. 102c 103c Starting with IOPT(2) options are specified by integers in the 104c range 0-6, followed in some cases by integers providing 105c argument(s). Each option, together with its options its 106c arguments if any is followed in IOPT by the next option (or 0). 107c 0 No more options; this must be last in the option list. 108c 1 An error estimate is to be returned in EOPT(1). 109c 2 (Argument: K2) K2 gives the polynomial degree to use when 110c extrapolating. 111c 3 (Argument: K3, L3) Save (k-th derivative of interpolating 112c polynomial) / (k!) in EOPT(K3+K-1) for k = 1, 2, ..., L3. 113c These values are the coefficients of the polynomial in the 114c monomial basis expanded about X. One must have 0<L3<NDEG+1. 115c 4 (Argument K4) The absolute and relative errors expected in 116c YT entries are specified in EOPT(K4) and EOPT(K4+1) 117c respectively. The values provided here are used in 118c estimating the error in the interpolation which is stored 119c in EOPT(1). 120c 5 (Argument K5, L5) Do the interpolation to the accuracy 121c requested by the absolute error tolerance specified in 122c EOPT(K5) and the relative error tolerance in EOPT(K5+1) 123c respectively. An attempt is made to keep the final error 124c < EOPT(K5) + EOPT(K5+1) * (abs(YT(i1) +abs(YT(i2)), where 125c i1 and i2 are indices for table values close to X. 126c Interpolation is done as for the default case, except that 127c in this case NDEG gives the maximal degree polynomial to use 128c in the interpolation, and polynomial interpolation of even 129c degree can happen. (The form that gives the C1 continuity 130c is not available in this case, and in fact continuity of the 131c interpolant itself is not to be expected. The actual degree 132c used in doing the interpolation is stored in the space for 133c the argument L5. An error estimate is returned in EOPT(1). 134c 6 (Argument K6) Do not use point k in the interpolation if 135c YT(k) = EOPT(K6). 136c 7 (Argument K7) YT(K7+i) gives the first derivative 137c corresponding to the function value in YT(i). These 138c derivatives are to be used in doing the interpolation. One 139c gets a continuous interpolant only for NDEG = 3, 7, 11, and 140c 15. The interpolating polynomial satisfies p(XT(i)) = YT(i), 141c p'(XT(i)) = YT(K7+i) for values of i that give values of XT 142c close to X. (Subject to keeping within one of the same 143c number of XT's on either side of X, the maximum of |XT(i)-X| 144c is minimized.) If NDEG is even a value of YT(i) is used 145c without using the corresponding value of YT(K7+i). 146c >253 Used for calls from DILUPM, a number of variables in the 147c common block have been set. If this is used, it must be the 148c first option, and common variables must be set. The rest 149c of IOPT is not examined for options. 150c EOPT Array used to return an error estimate and also used for 151c options. 152c EOPT(1) if an error estimate is returned, this contains a crude 153c estimate of the error in the interpolation. 154c 155c ************ External Procedures *********************** 156c 157c D1MACH Returns system parameters. 158c DMESS Prints error messages. 159c 160c ************ Variables referenced *********************** 161c 162c ADJSAV Saved difference of XT values used to adjust error estimate on 163c variable order when derivatives are being computed. 164c BADPT In common CDILUP. YT(K) is not to be used if YT(K) = BADPT, 165c and LBADPT = .true. 166c CY Array giving the YT values used to construct the interpolant. 167c Also used to store divided differences. 168c D1MACH Function to get parameters of floating point arithmetic. 169c DMESS Program calling MESS to print error messages. 170c DX Array giving the differences X - XT(I), where I runs through 171c the points that are used in the interpolation. 172c DXS Saved value of DX when computing derivatives 173c E1 Error estimate from one iteration ago (for variable order) 174c E2 Error estimate from this iteration (for variable order) 175c E2L Used in computing error estimate for variable order. 176c EBND If estimated error is < EBND then order is suff. high. 177c EBNDI Internal value for bounding error. 178c EBNDR Used in getting part of EBND due to relative accuracy request. 179c EF Factor used in estimating errors for variable order. 180c EN Used in computing error estimate for variable order. 181c EOPT Formal argument, see above. 182c EPSR Relative error level (D1MACH(4)). 183c ERRDAT Holds floating point data for error messages. 184c ERREST Estimate of the error in the interpolation. 185c GETERR Logical variable that is .true. if we are getting an error est. 186c H In indexed lookup contains the difference between XT values. 187c I Temporary index. 188c IO Index used in processing options. 189c IDX In common CDILUP. Gives indices of points selected for the 190c interpolation in order selected. 191c IIFLG Internal value for IOPT(1). 192c ILI Lower (sometimes upper) bound on the XT indices that will be 193c used in the interpolation. 194c INDXED Logical variable that is .true. if we have an indexed XT. 195c IOPT Formal argument, see above. 196c IUI As for ILI, except an upper (sometimes lower) bound. 197c K Index usually into YT and XT, into CY for derivatives. 198c KAOS In common CDILUP. Keeps track of state on variable order. 199c = 0 No variable order -- this value only set in DILUPM. 200c = 1 First time 201c = 2 Second time 202c = 3 Just had an increase or a very weak decrease in the error. 203c = 4 Just had a strong decrease in the error. 204c On an error with IIFLG .lt. 0, this is set to -10 * stop level - 205c the print level. 206c KDER 0 if not doing Hermite interpolation, else is < 0 if have next 207c higher order difference computed (in TP3), > 0 if it isn't. 208c KEXTRP In common CDILUP. Gives degree to use when extrapolating. 209c KGO Indicates type of interpolation as follows: 210c 1 Standard polynomial interpolation. 211c 2 Get error estimate after 1. 212c 3 Compute an interpolant with some desired accuracy. 213c 4 Compute an interpolant with a continuous derivative. 214c >10 Same as if KGO were 10 smaller, except XT values have already 215c been selected. 216c KGOC In common CDILUP. Value saved for -KGO. If YT contains values 217c of Y computed elsewhere in the order specified in IDX, then 218c KGOC should be set to abs(KGOC) before calling DILUP. If -2, 219c an extra value was computed for an error estimate. 220c KK Used in selecting data points to use. 221c 1 decrease ILI only 222c 0 increase IUI 223c -1 decrease ILI 224c <-1, increase IUI only 225c L Temporary index used in searching the XT array. Also one less 226c than the number of points used in the current interpolant. 227c LBADPT Logical variable set = .true. if checking for bad points. 228c LDER Offset of y' values from the y values. 229c LDERIV In common CDILUP. Loc. where first deriv. is stored in EOPT(). 230c LEXERR In common CDILUP. Loc. where absolute and relative error 231c information on Y is stored in EOPT. 232c LEXIT In common CDILUP. Defines action after finding location in XT. 233c 0 Don't compute anything for Y, just return. (Used for DILUPM.) 234c 1 Usual case, just interpolate. 235c >1 Compute LEXIT-1 derivatives of the interpolant. 236c LINC Logical variable used in the sequential search. Set .true. if 237c the values in XT are increasing with I, and set .false. otherwise. 238c LNDEG Location in IOPT to save degree when variable order is used. 239c LOPT Index of the last option. 240c LUP Formal argument, see above. 241c MACT Array used for error message actions, see DMESS for details. 242c MACTV Array used for part of error message for IOPT. 243c MAXDEG Parameter giving the maximum degree polynomial interpolation 244c supported. MAXDEG must be odd and > 2. 245c MESS Message program called from DMESS to print error messages. 246c MEEMES Parameter giving value for printing an error message in MESS. 247c MEIVEC Parameter giving value for printing an integer vector in MESS. 248c MEMDA1 Parameter giving value for indicating value in MACT for MESS. 249c MEMDA2 Parameter giving value for indicating value in MACT for MESS. 250c MEMDAT Parameter giving value for specifying MACT index for MESS. 251c MERET Parameter giving value to indicate no more actions for MESS. 252c METEXT Parameter giving value to tell MESS to print from MTXTAA. 253c LTXTxx Parameter names of this form were generated by PMESS in making 254c up text for error messages and are used to locate various parts 255c of the message. 256c MLOC Array giving starting loctions for error message text. 257c MTXTAA Character array holding error message text for MESS. 258c N Used in logic for deciding bounds. If N = 0, ILI can not be 259c reduced further. If N = 3*NTABI, IUI can not be increased further. 260c NDEG Formal argument, see above. 261c NDEGE Degree of polynomial to be used. Starts out = NDEG. 262c NDEGEC In common CDILUP. Degree actually used, saved in common. 263c NDEGI Degree of polynomial up to which y & differences are computed. 264c NDEGQ Used when KGO > 10. In this case If L .ge. NDEGQ special 265c action is needed. (Usually means quitting.) 266c NTAB Formal argument, see above. 267c NTABI Internal value of NTAB, = number of points in XT and YT. 268c PI Array use to store interpolation coefficients 269c PID Temporary storage used when computing derivatives. 270c SAVED Logical variable set = .true. if a DX value needs to be 271c replaced and restored when computing derivatives. 272c TP1 Used for temporary accumulation of values and temp. storage. 273c TP2 Used for temporary storage. 274c TP3 Used for divided difference when doing Hermite interpolation. 275c X Formal argument, see above. 276c XI Internal value for X, the place where interpolation is desired. 277c XL Value of "left hand" X when selecting indexed points. 278c XT Formal argument, see above. 279c XU Value of "right hand" X when selecting indexed points. 280c Y Formal argument, see above. 281c YL Last value of Y when selecting the order. 282c YT Formal argument, see above. 283c YTORD Logical variable set .true. when YT values are ordered on entry 284c 285c ************* Formal Variable Declarations *************** 286c 287 external D1MACH 288 integer IOPT(*), LUP, NDEG, NTAB 289 double precision X, XT(*), Y, YT(*), EOPT(*) 290c 291c ************* Common Block and Parameter ***************** 292c 293c MAXDEG must be odd and > 2. 294 integer MAXDEG 295 parameter (MAXDEG = 15) 296c 297 logical GETERR 298 integer KAOS, KEXTRP, KGOC,LDERIV,LEXERR,LEXIT,NDEGEC, 299 1 IDX(0:MAXDEG+1) 300 double precision BADPT, DX(0:MAXDEG+1), EBND, EBNDR 301 common / CDILUP / BADPT, DX, EBND, EBNDR, KAOS, KEXTRP, KGOC, 302 1 LDERIV, LEXERR,LEXIT,NDEGEC,IDX,GETERR 303c 304c ************* Local Variables **************************** 305c 306 logical INDXED, LBADPT, LINC, YTORD, SAVED 307 integer I, IIFLG, ILI, IUI, IO, K, KDER, KGO, KK, L, 308 1 LDER, LNDEG, LOPT, N, NDEGE, NDEGI, NDEGQ, NTABI 309 double precision ADJSAV, CY(0:MAXDEG+1), DXS, E1, E2, E2L, EBNDI, 310 1 EF, EM, EPSR, ERREST, H, PI(0:MAXDEG+1), 311 2 PID(0:MAXDEG), TP1, TP2, TP3, XI, XL, XU, YL 312 double precision D1MACH 313 save EPSR, /CDILUP/ 314c 315c ************************ Error Message Stuff and Data **************** 316c 317c Parameter defined below are all defined in the error message program 318c DMESS. 319c 320 integer MECONT,MEMDA1,MERET,MEEMES,METEXT,MEIVEC 321 parameter (MEMDA1 =27) 322 parameter (MECONT = 50) 323 parameter (MERET =51) 324 parameter (MEEMES =52) 325 parameter (METEXT =53) 326 parameter (MEIVEC =57) 327c 328 double precision ERRDAT(2) 329 integer MLOC(9), IDAT(4), MACT(5), MACTV(6) 330 331c ********* Error message text *************** 332c[Last 2 letters of Param. name] [Text generating message.] 333cAA DILUP$B 334cAB Estimated error = $F, Requested error = $F.$E 335cAC NTAB = 1, no error estimate computed.$E 336cAD Too many bad points, only $I point(s) available.$E 337cAE XT(1) = XT(NTAB=$I) = $F ??$E 338cAF LUP = 3, and XT(2) = 0.$E 339cAG Order of requested derivative, $I, is > bound of $I.$E 340cAH LUP > 3, with unprepared common block.$E 341cAI NDEG = $I must be .le. $I or NTAB = $I must be at $C 342c least $I. You must either decrease NDEG, increase NTAB, $C 343c or if requesting an error estimate, turn off that request.$E 344cAJ NDEG = $I is not in the allowed interval of [0, $I].$E 345cAK NTAB = $I is not in the allowed interval of [1, $I].$E 346cAL IOPT($I) = $I is not a valid option.$E 347c $ 348cAM IOPT(1:$M): $E 349 integer LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG,LTXTAH, 350 * LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM 351 parameter (LTXTAA= 1,LTXTAB= 8,LTXTAC= 53,LTXTAD= 92,LTXTAE=142, 352 * LTXTAF=171,LTXTAG=196,LTXTAH=250,LTXTAI=290,LTXTAJ=458, 353 * LTXTAK=512,LTXTAL=566,LTXTAM= 1) 354 character MTXTAA(3) * (201) 355 character MTXTAB(1) * (14) 356 data MTXTAA/'DILUP$BEstimated error = $F, Requested error = $F.$EN 357 *TAB = 1, no error estimate computed.$EToo many bad points, only $I 358 * point(s) available.$EXT(1) = XT(NTAB=$I) = $F ??$ELUP = 3, and XT 359 *(2) = 0.$EOrder ','of requested derivative, $I, is > bound of $I.$ 360 *ELUP > 3, with unprepared common block.$ENDEG = $I must be .le. $I 361 * or NTAB = $I must be at least $I. You must either decrease NDEG 362 *, increase NTAB, or if',' requesting an error estimate, turn off t 363 *hat request.$ENDEG = $I is not in the allowed interval of [0, $I]. 364 *$ENTAB = $I is not in the allowed interval of [1, $I].$EIOPT($I) = 365 * $I is not a valid option.$E'/ 366 data MTXTAB/'IOPT(1:$M): $E'/ 367c **** End of automatically generated text 368c 369c 1 2 3 4 5 370 data MACT / MEEMES,0,0,0, MECONT / 371c 2 3 4 5 6 372 data MACTV /MEMDA1, 0, METEXT, MEIVEC, 0, MERET / 373 374c 375 data MLOC / LTXTAB, LTXTAC, LTXTAD, LTXTAE, LTXTAF, LTXTAG, 376 1 LTXTAH, LTXTAI, LTXTAJ / 377c 378 data EPSR / 0.D0 / 379c 380c 381c ************ Start of executable code ******************* 382c 383 IIFLG = 0 384 KDER = 0 385 PI(0) = 1.D0 386 LBADPT = .false. 387 ILI = -LUP 388 if (ILI .le. -4) go to 1400 389 LEXERR = 0 390 NTABI = NTAB 391 if ((NTABI .le. 0) .or. (NTABI .gt. 9999999)) then 392 IIFLG = -9 393 IDAT(1) = NTABI 394 IDAT(2) = 9999999 395 go to 2120 396 end if 397 NDEGI = NDEG 398 NDEGE = NDEGI 399 XI = X 400 KGO = 1 401 IO = 1 402 if (IOPT(2) .ge. 254) then 403 LBADPT = IOPT(2) .eq. 255 404 if (KAOS .le. 0) go to 50 405 LNDEG = 0 406 KGO = 3 407 KAOS = 1 408 NDEGI = min(MAXDEG, IOPT(3)+1) 409 NDEGE = -NDEGI 410 go to 60 411 end if 412 GETERR = .false. 413 KEXTRP = -1 414 LEXIT = 1 415c Loop to take care of options 416 10 IO = IO + 1 417 LOPT = IOPT(IO) 418 if (LOPT .ne. 0) then 419 go to (1930, 1500, 1600, 1900, 1950, 1970, 1980), LOPT 420 IIFLG = -10 421 IDAT(1) = IO 422 IDAT(2) = LOPT 423 go to 2120 424 end if 425 50 continue 426 K = NDEGI + 1 427 if (KDER .ne. 0) then 428 K = K / 2 429 else if (mod(NDEGI,2) .eq. 0) then 430 K = K + 1 431 end if 432 if (GETERR) K = K + 1 433 if ((NDEGI.lt.0) .or. (NDEGI.gt.MAXDEG) .or. (K.gt.NTAB)) then 434 IIFLG = -8 435 IDAT(1) = NDEGI 436 IDAT(2) = MAXDEG 437 IDAT(3) = NTAB 438 IDAT(4) = K 439 go to 2120 440 end if 441 if (KGO .eq. 1) then 442 if (NDEGI .ge. 2) then 443 if (mod(NDEGI, 2) .eq. 0) then 444 if (KDER .eq. 0) then 445 NDEGE = NDEGE + 1 446 KGO = 4 447 end if 448 end if 449 end if 450 end if 451 60 if (KEXTRP .lt. 0) then 452 KEXTRP = max(NDEGE-1, min(NDEGE, 2)) 453 if (KEXTRP .lt. 0) KEXTRP = NDEGI 454 end if 455 if (ILI .gt. 0) go to 200 456 if (ILI + 2) 1300, 1200, 100 457c 458c Binary search, then sequential 459 100 continue 460c In binary search XT(ILI) .le. XI .le. XT(IUI) 461c or we are extrapolating 462 ILI = 1 463 IUI = NTABI 464 if (XT(NTABI) - XT(1)) 110, 1220, 130 465 110 ILI = NTABI 466 L = 1 467 120 IUI = L 468 130 L = (IUI - ILI) / 2 469 if (L .eq. 0) go to 210 470 L = ILI + L 471 if (XT(L) .gt. XI) go to 120 472 ILI = L 473 go to 130 474c 475c Sequential search, then exit 476 200 continue 477 if (ILI .gt. NTABI) go to 100 478 if (XT(NTABI) .eq. XT(1)) go to 1220 479 210 INDXED = .false. 480 LINC = XT(NTABI) .gt. XT(1) 481 if ((XT(ILI) .gt. XI) .eqv. LINC) go to 230 482 220 if (ILI .eq. NTABI) go to 1230 483 ILI = ILI + 1 484 if ((XT(ILI) .lt. XI) .eqv. LINC) go to 220 485 N = 2*ILI 486 if (abs(XT(ILI-1)-XI) .lt. abs(XT(ILI) - XI)) then 487 ILI = ILI - 1 488 N = N - 1 489 end if 490 go to 240 491 230 if (ILI .eq. 1) go to 1230 492 ILI = ILI - 1 493 if ((XT(ILI) .gt. XI) .eqv. LINC) go to 230 494 N = 2*ILI + 1 495 if (abs(XT(ILI+1)-XI) .lt. abs(XT(ILI) - XI)) then 496 ILI = ILI + 1 497 N = N + 1 498 end if 499 240 if (LUP .le. 0) LUP = -ILI 500 250 DX(0) = XI - XT(ILI) 501c Get bounding indices and interpolate 502 260 IUI = ILI 503 K = ILI 504 L = -1 505 if (LEXIT .eq. 0) then 506 NDEGI = 0 507 if (GETERR) then 508 if (KGO .ne. 3) then 509 if (KGO .eq. 1) NDEGE = NDEGE + 1 510 KEXTRP = min(KEXTRP+1, NDEGE) 511 else 512 NDEGE = -NDEGE 513 end if 514 end if 515 end if 516c Just got index for next XT 517 270 if (LEXIT .eq. 0) then 518 L = L + 1 519 IDX(L) = K 520 go to 350 521 end if 522 TP2 = YT(K) 523 if (LBADPT) then 524c Check if point should be discarded 525 if (TP2 .eq. BADPT) then 526 if (IUI .ne. ILI) then 527 if (abs(KK) .eq. 1) then 528 N = N - 1 529 else 530 N = N + 1 531 end if 532 end if 533 go to 1020 534 end if 535 end if 536 280 L = L + 1 537 IDX(L) = K 538c 539 290 if (KDER .ne. 0) then 540c Hermite interpolation 541 KDER = -KDER 542 if (KDER .lt. 0) then 543 TP3 = YT(LDER + K) 544 if (L .eq. 0) then 545 TP1 = TP2 546 else 547 do 300 I = 1, L 548 TP2 = (TP2 - CY(I-1)) / (DX(I-1) - DX(L)) 549 TP3 = (TP3 - TP2) / (DX(I-1) - DX(L)) 550 300 continue 551 PI(L) = PI(L-1) * DX(L-1) 552 TP1 = TP1 + PI(L) * TP2 553 end if 554 else 555 DX(L) = DX(L-1) 556 PI(L) = PI(L-1) * DX(L) 557 TP2 = TP3 558 TP1 = TP1 + PI(L) * TP2 559 end if 560 else if (L .eq. 0) then 561 TP1 = TP2 562 else 563 PI(L) = PI(L-1) * DX(L-1) 564 if (L .le. NDEGI) then 565c Get divided differences & interpolate. 566 do 320 I = 1, L 567 TP2 = (TP2 - CY(I-1)) / (DX(I-1) - DX(L)) 568 320 continue 569 TP1 = TP1 + PI(L) * TP2 570 end if 571 end if 572 CY(L) = TP2 573 350 if (L .lt. NDEGE) go to 1000 574 360 if (LEXIT .eq. 0) go to 2110 575 go to (500, 600, 900, 800), KGO 576 if (L .ge. NDEGQ) go to (500,600,900,800), KGO-10 577c Already got the points selected. 578 400 L = L + 1 579 if (YTORD) then 580 TP2 = YT(L+1) 581 else 582 K = IDX(L) 583 TP2 = YT(K) 584 end if 585 go to 290 586c Got simple interpolated value. 587 500 Y = TP1 588 if (.not. GETERR) go to 2000 589 NDEGI = NDEGI + 1 590 KGO = KGO + 1 591 if (NDEGE .ge. 0) go to 1000 592 go to 400 593c Got info. for error estimate. 594 600 if (L .eq. 0) then 595 IIFLG = -2 596 else 597 ERREST=1.5D0*(abs(TP1-Y)+.03125D0*abs(PI(L-1)*CY(L-1))) 598 L = L - 1 599 end if 600 go to 2000 601c C1 interpolant. 602 800 do 810 I = 1, L-1 603 TP2 = (TP2 - CY(I-1)) / (DX(I-1) - DX(L)) 604 810 continue 605 TP2 = (TP2 - CY(L-1)) / (DX(0) - DX(1)) 606 CY(L) = TP2 607 PI(L) = PI(L-1) * DX(0) 608 Y = TP1 + PI(L) * TP2 609 if (GETERR) ERREST = 1.5D0 * abs(PI(L-1)) * abs(((DX(L-1) * 610 1 (DX(0) - DX(1)) / (DX(L-1) - DX(L)) - DX(0)) * TP2) + 611 2 abs(.03125D0*CY(L-1))) 612 go to 2000 613c Variable order, check estimated error and convergence. 614 900 continue 615 if (KAOS .ge. 3) go to 930 616 if (KAOS .eq. 2) go to 920 617c First time 618 if ((KDER.ne.0) .and. (LEXIT.gt.1) .and. (L.eq.0)) go to 970 619 KAOS = 2 620 E2L = abs(TP1) 621 E2 = EBND + 1.D30 622 go to 970 623c Second time 624 920 KAOS = 4 625 EBNDI = .66666D0*(EBND + EBNDR*(abs(CY(0))+abs(YT(IDX(1))))) 626 EF = DX(0) 627 if (LEXIT .gt. 1) then 628 EF = DX(L) - DX(0) 629 ADJSAV = EF 630 end if 631 E2 = abs(EF * TP2) 632 EM = .75D0 633 go to 950 634c Usual case 635 930 E2L = E2 636 E1 = E2L * (5.D0 * EM / dble(L)) 637 EF = EF * DX(L-1) 638 E2 = abs(EF*TP2) 639 EM = 0.5D0 * EM + E2 / (E2L + E2 + 1.D-20) 640 if (E2 .ge. E1) then 641 if (KAOS .eq. 3) then 642c Apparently diverging, so quit. 643 L = L - 1 644 go to 960 645 else 646 KAOS = 3 647 end if 648 else 649 KAOS = 4 650 end if 651 950 YL = TP1 652 if (E2L + E2 .gt. EBNDI) then 653 if ((L+NDEGE .lt. 0) .and. (IIFLG .ne. 2)) go to 970 654 if (KGO .eq. 13) then 655c May need early exit to get next point. 656 IOPT(2) = 0 657 if (L .lt. NDEG) go to 2110 658 end if 659 end if 660 960 TP2 = 1.5D0 661 if (LEXIT .gt. 1) then 662 if (KDER .eq. 0) TP2 = 1.5D0 * abs(DX(0) / ADJSAV) 663 end if 664 ERREST = TP2 * (E2 + .0625D0 * E2L) 665 if (LNDEG .ne. 0) IOPT(LNDEG) = L 666 Y = YL 667 go to 2000 668 970 if (KGO .gt. 10) go to 400 669c 670 1000 if (KDER .lt. 0) go to 280 671 1020 KK = min(N - IUI - ILI, 2) - 1 672c In this section of code, KK=: 1, decrease ILI only; 0, increase IUI; 673c -1, decrease ILI; and <-1, increase IUI only 674 if (abs(KK) .eq. 1) then 675 ILI = ILI - 1 676 K = ILI 677 if (ILI .ne. 0) then 678 if (INDXED) then 679 XL = XL + H 680 DX(L+1) = XL 681 go to 270 682 end if 683 DX(L+1) = XI - XT(K) 684 if (XT(ILI+1) .ne. XT(ILI)) go to 270 685 end if 686 if (KK .eq. 1) go to 1100 687 N = 0 688c If all points will be on one side, flag extrapolation. 689c if (L .le. 0) go to 1250 690 else 691 IUI = IUI + 1 692 K = IUI 693 if (IUI .le. NTABI) then 694 if (INDXED) then 695 XU = XU - H 696 DX(L+1) = XU 697 go to 270 698 end if 699 DX(L+1) = XI - XT(K) 700 if (XT(IUI-1) .ne. XT(IUI)) go to 270 701 end if 702 if (KK .ne. 0) go to 1100 703 N = 3*NTABI 704c If all points will be on one side, flag extarpolation. 705c if (L .le. 0) go to 1250 706 end if 707 if (KGO .lt. 4) go to 1020 708c If we can't get the same number of points on either side, the 709c continuous derivative case becomes standard interpolation. 710 KGO = 1 711 NDEGE = NDEGE - 1 712 if (L .lt. NDEGE) go to 1020 713 go to 360 714c No more data accessible. 715 1100 if (L .le. 0) then 716c Too many bad points, couldn't find points. 717 Y = TP1 718 IDAT(1) = L + 1 719 IIFLG = -3 720 go to 2110 721 end if 722 NDEGI = min(NDEGI, L) 723 NDEGE = 0 724 IIFLG = 2 725 go to 360 726c 727c Secant start, then use sequential 728 1200 continue 729 if (XT(1) .eq. XT(NTABI)) go to 1220 730 ILI = max(1, min(NTABI, int(1.5D0+dble(NTABI-1)*(XI-XT(1))/ 731 1 (XT(NTABI) - XT(1))))) 732 go to 210 733c 734c Special cases 735c 1 entry in XT 736 1220 ILI = 1 737 IUI = 1 738 KGO = 1 739 K = 1 740 if (NDEGE .ne. 0) IIFLG = 2 741 NDEGE = 0 742 if (NTABI .eq. 1) go to 250 743c Error -- XT(1) .eq. XT(NTAB), and NTAB .ne. 1 744 IIFLG = -4 745 IDAT(1) = NTABI 746 ERRDAT(1) = XT(1) 747 KGO = 0 748 go to 2120 749c Check if really extrapolating 750 1230 N = 6 * (ILI - 1) 751 if ((XI - XT(1)) * (XT(NTABI) - XI) .ge. 0.D0) go to 240 752c Extrapolating 753 1250 IIFLG = 1 754 if (KGO .ge. 4) KGO = 1 755 NDEGI = KEXTRP 756 NDEGE = sign(NDEGI, NDEGE) 757 if (INDXED) go to 260 758 go to 240 759c 760c Index search, then exit 761 1300 continue 762 INDXED = .true. 763 H = XT(2) 764 if (H .eq. 0) then 765 IIFLG = -5 766 go to 2120 767 end if 768 TP1 = 1.D0 + (XI - XT(1)) / H 769 ILI = min(max(1, nint(TP1)), NTABI) 770 XU = (TP1 - dble(ILI)) * H 771 XL = XU 772 DX(0) = XU 773 N = ILI + ILI 774 if (TP1 .gt. dble(ILI)) N = N + 1 775 if ((XI - XT(1)) * ((XT(1) + H*dble(NTABI-1) - XI)) .ge. 0.E0) 776 1 go to 260 777 N = 6 * (ILI - 1) 778 go to 1250 779c 780c Already set up 781 1400 KGO = KGOC 782 if (KGO .eq. 0) then 783 IIFLG = -7 784 go to 2120 785 end if 786 YTORD = KGO .gt. 0 787 KGO = abs(KGO) 788 if (KGO .lt. 10) KGO = KGO + 10 789 if (KGO .eq. 12) KGO = 11 790 NDEGI = NDEGEC 791 NDEGQ = NDEGI 792 if (KGO .eq. 13) NDEGQ = 0 793 NDEGE = -NDEGI 794 L = -1 795 if (KGO .eq. 14) NDEGI = NDEGI - 1 796 go to 400 797c 798c Set number of points for extrapolation 799 1500 continue 800 I = I + 1 801 KEXTRP = min(IOPT(I), NDEGE) 802 go to 10 803c 804c Get derivatives of interpolant 805 1600 continue 806 I = I + 2 807 LDERIV = IOPT(I-1) 808 LEXIT = IOPT(I) + 1 809 go to 10 810c 811c Get expected errors in YT. 8121900 continue 813 I = I + 1 814 LEXERR = IOPT(I) 815c 816c Set to get error estimate. 817 1930 continue 818 GETERR = .true. 819 go to 10 820c 821c Set for automatic order selection. 822 1950 continue 823 IO = IO + 2 824 LNDEG = IO 825 EBND = max(0.D0, EOPT(IOPT(IO-1))) 826 EBNDR = max(0.D0, EOPT(IOPT(IO-1)+1)) 827 KGO = 3 828 KAOS = 1 829 NDEGI = min(MAXDEG, NDEGE+1) 830 NDEGE = -NDEGI 831 go to 1930 832c 833c Set to ignore special points 834 1970 continue 835 LBADPT = .true. 836 IO = IO + 1 837 BADPT = EOPT(IOPT(IO)) 838 go to 10 839c 840c Set up for Hermite interpolation. 841 1980 continue 842 IO = IO + 1 843 LDER = IOPT(IO) 844 KDER = abs(LDER) 845 go to 10 846c 847c Put any new options just above here 848c End of the loop 849c 850 2000 if (LEXIT .lt. 2) go to 2100 851c Compute derivatives of Y 852 SAVED = mod(KGO, 10) .eq. 4 853 if (SAVED) then 854 DXS = DX(L-1) 855 DX(L-1) = DX(0) 856 end if 857 N = LEXIT - 1 858 if (N .gt. L) then 859 if (N .gt. MAXDEG) then 860 IDAT(1) = N 861 IDAT(2) = MAXDEG 862 N = MAXDEG 863 IIFLG = -6 864 end if 865 do 2020 I = L+1, N 866 EOPT(I+LDERIV-1) = 0.D0 867 2020 continue 868 N = L 869 end if 870 TP1 = CY(1) 871 PID(0) = 1.D0 872 do 2030 I = 1, L-1 873 PID(I) = PI(I) + DX(I) * PID(I-1) 874 TP1 = TP1 + PID(I) * CY(I+1) 875 2030 continue 876 EOPT(LDERIV) = TP1 877 do 2070 K = 2, N 878 TP1 = CY(K) 879 do 2060 I = 1, L-K 880 PID(I) = PID(I) + DX(I+K-1) * PID(I-1) 881 TP1 = TP1 + PID(I) * CY(I+K) 882 2060 continue 883 EOPT(K+LDERIV-1) = TP1 884 2070 continue 885 if (SAVED) DX(L-1) = DXS 886c Save info. and return 887 2100 continue 888 if (GETERR) then 889 if (EPSR .eq. 0.D0) then 890 EPSR = D1MACH(4) 891 end if 892 TP1 = EPSR 893 TP2 = 0.D0 894 if (LEXERR .ne. 0) then 895 TP2 = max(0.D0, EOPT(LEXERR)) 896 TP1 = max(TP1, EOPT(LEXERR+1)) 897 end if 898 if (IIFLG .eq. 2) ERREST = 32.D0 * ERREST 899 EOPT(1) = ERREST + TP2 + TP1*(abs(CY(0)) + abs(DX(0)*CY(1))) 900 if ((KGO .eq. 3) .or. (KGO .eq. 13)) then 901 if (EBNDI .ne. 0.D0) then 902 if (EOPT(1) .gt. EBND) then 903 ERRDAT(1) = EOPT(1) 904 ERRDAT(2) = EBND 905 IIFLG = -1 906 end if 907 end if 908 end if 909 end if 910 2110 KGOC = -KGO 911 NDEGEC = L 912 2120 IOPT(1) = IIFLG 913 if (IIFLG .ge. 0) return 914c Take care of error messages, IIFLG < -2 should stop. 915 MACT(2) = 88 916 if (IIFLG .lt. -1) KAOS = -MACT(2) 917 if (IIFLG .ge. -3) MACT(2) = min(26, 24 - IIFLG) 918 if (IOPT(2) .ge. 254) then 919 if (IIFLG .eq. -1) return 920 MACT(2) = 28 921 end if 922 MACT(3) = -IIFLG 923 MACT(4) = MLOC(-IIFLG) 924 call DMESS(MACT, MTXTAA, IDAT, ERRDAT) 925 MACTV(2) = IO 926 MACTV(5) = IO 927 K = 1 928 if (IO .eq. 1) K = 6 929 call MESS(MACTV(K), MTXTAB, IOPT) 930 return 931c End of Interpolation routine -- DILUP 932 end 933