1      double precision function DHINT (X,NDERIV,NTAB,XTAB,YTAB,YPTAB)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5C>> 1994-11-11 DHINT  Krogh   Declared all vars.
6C>> 1994-10-20 DHINT  Krogh  Changes to use M77CON
7C>> 1987-12-09 DHINT  Lawson  Initial code.
8c--D replaces "?": ?HINT
9c
10c          This subprogram does look-up and Hermite cubic interpolation
11c     to evaluate a function at X defined by the data
12c     NTAB, XTAB(), YTAB(), & YPTAB().
13c     The values in XTAB() must be either strictly increasing or
14c     strictly decreasing.
15c     Optionally this subprogram will evaluate either the function or
16c     its first, second, or third derivative.
17c
18c     The look-up method is a linear search beginning with the index,
19c     LOOK, saved internally from the previous reference to this
20c     subprogram.  LOOK is reset if it exceeds the current NTAB.
21c     Evaluation at an interior tabular point uses data from the
22c     interval in the direction of decreasing abcissa values.
23c     Evaluation outside the table (extrapolation) uses data from the
24c     nearest tabular interval.
25c     ------------------------------------------------------------------
26c     Based on subprograms FXDP and DYDXDP by Lawson, Hanson, Lang, and
27c     Campbell, JPL, 1968-1974.
28c     Modified July 1984 and July 1987 by C. L. Lawson, JPL, for
29c     inclusion in the MATH 77 library.
30c     ------------------------------------------------------------------
31c     X  [in]  Argument at which interpolation is to be done.
32c           For best results X should be between XTAB(1) and XTAB(NTAB),
33c           however this subprogram will give an extrapolated
34c           value when X is outside this range.
35c
36c     NDERIV  [in]   = 0 to compute function value.
37c                    = 1 to compute value of first derivative.
38c                    = 2 to compute value of second derivative.
39c                    = 3 to compute value of third derivative.
40c
41c     NTAB  [in]  No. of points in XTAB(), YTAB(), and YPTAB().
42c                 Require NTAB .ge. 2.
43c
44c     XTAB()  [in]  Array of abcissas.  Values must either be strictly
45c                   increasing or strictly decreasing.
46c
47c     YTAB()  [in]  Array of function values.
48c                   YTAB(i) = f(XTAB(i))
49c
50c     YPTAB()  [in]  Array of first derivative values.
51c                    YPTAB(i) = f'(XTAB(i))
52c     ------------------------------------------------------------------
53      integer NDERIV, NTAB, LOOK, LOOK1
54      double precision XTAB(NTAB), YTAB(NTAB), YPTAB(NTAB)
55      double precision ONE, TWO, THREE, SIX, Q0, Q1, Q2, Q3, Q4, Q5,R,X
56      parameter(ONE = 1.0D0, TWO = 2.0D0, THREE = 3.0D0, SIX = 6.0D0)
57      save LOOK
58      data  LOOK/1/
59c     ------------------------------------------------------------------
60C
61C
62      if(XTAB(1) .lt. XTAB(2)) then
63c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64c                             Here when values in XTAB() are increasing.
65      IF( LOOK .GE. NTAB ) LOOK = max(1,NTAB/2)
66C                                      Here  1 .le. LOOK .le. NTAB-1
67      IF(X .gt. XTAB(LOOK+1)) GO TO 30
68C                 Search toward decreasing abcissae, decreasing indices.
69   25 IF(LOOK .eq. 1 ) GO TO 40
70      IF(X .gt. XTAB(LOOK)) GO TO 40
71      LOOK=LOOK-1
72      GO TO 25
73C                 Search toward increasing abcissae, increasing indices.
74   30 IF(LOOK .eq. NTAB-1) GO TO 40
75   35 LOOK=LOOK+1
76      IF(LOOK .eq. NTAB-1) GO TO 40
77      IF(X .gt. XTAB(LOOK+1)) GO TO 35
78   40 LOOK1 = LOOK+1
79      else
80c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81c                             Here when values in XTAB() are decreasing.
82      IF( LOOK .gt. NTAB .or. LOOK .lt. 2 ) LOOK = max(2,NTAB/2)
83C
84C                         Here  2 .le. LOOK .le. NTAB
85C
86      IF(X .gt. XTAB(LOOK-1)) GO TO 50
87C                 Search toward decreasing abcissae, increasing indices.
88   45 IF(LOOK .eq. NTAB ) GO TO 60
89      IF(X .gt. XTAB(LOOK)) GO TO 60
90      LOOK=LOOK+1
91      GO TO 45
92C                 Search toward increasing abcissae, decreasing indices.
93   50 IF(LOOK .eq. 2) GO TO 60
94   55 LOOK=LOOK-1
95      IF(LOOK .eq. 2) GO TO 60
96      IF(X .gt. XTAB(LOOK-1)) GO TO 55
97   60 LOOK1 = LOOK-1
98      endif
99c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100c              Here XTAB(LOOK) < XTAB(LOOK1) and X is in this interval
101c              unless we are extrapolating.
102c              Evaluate the Hermite interpolation formulas.
103C
104      R=XTAB(LOOK1)-XTAB(LOOK)
105      Q0=(X-XTAB(LOOK))/R
106      Q1=Q0-ONE
107C
108      go to (100, 101, 102, 103), NDERIV+1
109C
110  100 Q4=TWO*Q0
111      Q2=ONE+Q4
112      Q3=THREE-Q4
113      Q4=Q1**2
114      Q5=Q0**2
115      DHINT  =Q4*(Q2*YTAB(LOOK)+R*Q0*YPTAB(LOOK))+
116     *        Q5*(Q3*YTAB(LOOK1)+R*Q1*YPTAB(LOOK1))
117      return
118C
119  101 Q2=THREE*Q0
120      Q3=SIX*Q0/R
121      DHINT  =Q1*(Q3*(YTAB(LOOK)-YTAB(LOOK1))+
122     *        (Q2-ONE)*YPTAB(LOOK)) + Q0*(Q2-TWO)*YPTAB(LOOK1)
123      return
124c
125  102 Q2 = Q1 + Q0
126      DHINT = (TWO/R) * ( (THREE/R) * (YTAB(LOOK)-YTAB(LOOK1)) * Q2 +
127     *        YPTAB(LOOK)*(Q2+Q1) + YPTAB(LOOK1)*(Q2+Q0) )
128
129      return
130c
131  103 DHINT = (SIX/R**2) * ( (TWO/R) * (YTAB(LOOK)-YTAB(LOOK1)) +
132     *                       YPTAB(LOOK) + YPTAB(LOOK1) )
133      return
134C
135      end
136