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