C ALGORITHM 697 , COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 3, SEPTEMBER, 1991, PP. 367. PROGRAM TTUVIP C C Test Program for the UVIP3P Subroutine C C Hiroshi Akima C U.S. Department of Commerce, NTIA/ITS C Version of 89/07/04 C C This program tests the UVIP3P subroutine by calling UVIP3P, in C a sequence, with two values of the NP argument, which is the C degree of the polynomial for the interval between each pair of C data points. With each NP value, this program calls UVIP3P in C three ways: C (1) in repeated calls with interpolation at a point in C each call, C (2) in one call with interpolation at all points, and C (3) in one call as in (2) but with inverted data set. C The program then compares the results with the expected values, C precalculated and stored in the program, and writes the C differences in two pages. If all entries in the last three C columns on each page are all zeros, UVIP3P is considered to be C working as expected. Otherwise, something is wrong in UVIP3P C and/or the test program itself. C C Dummy arguments of the UVIP3P subroutine are described in the C comments in the subroutine. C C Variables and arrays used in this program are C NP = degree of the polynomial for the interval between C each pair of data points, C NPO = integer array of dimension 2 containing the two C NP values, C ND = number of data points, C NI = number of interpolated points, C XD1, YD1 C = arrays of dimension ND containing the abscissas C and ordinates of the data points, C XI1 = array of dimension NI containing the abscissas of C the desired points, C YIE = doubly-dimensioned array of dimension (NI,2) C containing the expected values of the interpolated C results for the two NP values in two columns, C YI1 = array of dimension NI, where the ordinate values C interpolated from the XD1, YD1, and XI1 values in C repeated calls are to be stored, C YI2 = array of dimension NI, where the ordinate values C interpolated from the XD1, YD1, and XI1 values in C one call are to be stored, C XD3, YD3 C = arrays of dimension ND, where the abscissas and C ordinates of the inverted data points are to be C stored, C XI3 = array of dimension NI, where the abscissas of the C desired points for the inverted data points are to C be stored, C YI3 = array of dimension NI, where the ordinate values C interpolated from the XD3, YD3, and XI3 values in C one call are to be stored, C DYI1, DYI2, DYI3 C = arrays of dimension NI, where the differences C between the YI1, YI2, and YI3 values and their C respective expected values are to be stored. C C C Specification statements DIMENSION NPO(2) PARAMETER (ND=10, NI=31) DIMENSION XD1(ND),YD1(ND),YIE(NI,2),XI1(NI) DIMENSION XD3(ND),YD3(ND),XI3(NI) DIMENSION YI1(NI),YI2(NI),YI3(NI) DIMENSION DYI1(NI),DYI2(NI),DYI3(NI) C Data statements DATA NPO/ 3, 6/ DATA XD1/ 1 1.0, 2.0, 4.0, 6.5, 8.0,10.0,10.5,11.0,13.0,14.0/ DATA YD1/ 1 0.0, 0.0, 0.0, 0.0, 0.1, 1.0, 4.5, 8.0,10.0,15.0/ DATA XI1/ 1 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 2 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 3 8.0, 8.5, 9.0, 9.5,10.0,10.5,11.0,11.5, 4 12.0,12.5,13.0,13.5,14.0,14.5,15.0/ DATA YIE/ 1 8*0.0, 2 6*0.0, 0.015, 0.052, 3 0.100, 0.036,-0.045, 0.172, 1.000, 4.500, 8.000,10.075, 4 10.705,10.483,10.000,11.204,15.000,19.767,24.533, 1 8*0.0, 2 6*0.0, 0.020, 0.057, 3 0.100, 0.134, 0.166, 0.314, 1.000, 4.500, 8.000, 9.689, 4 10.101,10.180,10.000,11.663,15.000,19.767,24.533/ C Opens the file (for the UNIX operating system). C 10 OPEN (6,FILE='TTUOUT',FORM='PRINT') C Calculation C Do-loop with respect to the two NP values 20 DO 79 INP=1,2 NP=NPO(INP) C Interpolation in repeated calls at a point in each call 30 DO 31 II=1,NI CALL UVIP3P(NP,ND,XD1,YD1,1,XI1(II),YI1(II)) 31 CONTINUE C Interpolation in one call at all points CALL UVIP3P(NP,ND,XD1,YD1,NI,XI1,YI2) C Interpolation in one call with inverted data 40 DO 41 ID=1,ND IDI=ND+1-ID XD3(IDI)=XD1(1)+XD1(ND)-XD1(ID) YD3(IDI)=YD1(ID) 41 CONTINUE DO 42 II=1,NI XI3(II)=XD1(1)+XD1(ND)-XI1(II) 42 CONTINUE CALL UVIP3P(NP,ND,XD3,YD3,NI,XI3,YI3) C Calculation of the differences between the calculated and C expected values 50 DO 51 II=1,NI DYI1(II)=YI1(II)-YIE(II,INP) DYI2(II)=YI2(II)-YIE(II,INP) DYI3(II)=YI3(II)-YIE(II,INP) 51 CONTINUE C Writing the results 60 IF (INP.EQ.1) THEN WRITE (*,6060) NP ELSE WRITE (*,6061) NP END IF WRITE (*,6062) DO 69 II=1,NI IF (MOD(II,4).EQ.2) WRITE (*,6063) IF (II.LE.ND) THEN ID=II WRITE (*,6064) XD1(ID),YD1(ID), 1 XI1(II),YIE(II,INP),YI1(II), 2 DYI1(II),DYI2(II),DYI3(II) ELSE WRITE (*,6065) XI1(II),YIE(II,INP),YI1(II), 1 DYI1(II),DYI2(II),DYI3(II) END IF 69 CONTINUE 79 CONTINUE STOP C Format statements 6060 FORMAT (1H , 1 'TTUVIP',14X,'Test Program for UVIP3P',29X,'NP =',I2) 6061 FORMAT (1H1, 1 'TTUVIP',14X,'Test Program for UVIP3P',29X,'NP =',I2) 6062 FORMAT (1X/// 1 6X,'Data Points',9X,'Interpolated Points ', 2 8X,' Differences '// 3 6X,'XD1 YD1',8X,'XI1 YIE YI1', 4 8X,'DYI1 DYI2 DYI3'/) 6063 FORMAT (1X) 6064 FORMAT (1X,2F9.3,2X,3F9.3,2X,3F9.3) 6065 FORMAT (19X,2X,3F9.3,2X,3F9.3) END SUBROUTINE UVIP3P(NP,ND,XD,YD,NI,XI, YI) C C Univariate Interpolation (Improved Akima Method) C C Hiroshi Akima C U.S. Department of Commerce, NTIA/ITS C Version of 89/07/04 C C This subroutine performs univariate interpolation. It is based C on the improved A method developed by Hiroshi Akima, 'A method C of univariate interpolation that has the accuracy of a third- C degree polynomial,' ACM TOMS, vol. xx, pp. xxx-xxx, 19xx. (The C equation numbers referred to in the comments below are those in C the paper.) C C In this method, the interpolating function is a piecewise C function composed of a set of polynomials applicable to C successive intervals of the given data points. This method C uses third-degree polynomials as the default, but the user has C an option to use higher-degree polynomial to reduce undulations C in resulting curves. C C This method has the accuracy of a third-degree polynomial if C the degree of the polynomials for the interpolating function is C set to three. C C The input arguments are C NP = degree of the polynomials for the interpolating C function, C ND = number of input data points C (must be equal to 2 or greater), C XD = array of dimension ND, containing the abscissas of C the input data points C (must be in a monotonic increasing order), C YD = array of dimension ND, containing the ordinates of C the input data points, C NI = number of points for which interpolation is desired C (must be equal to 1 or greater), C XI = array of dimension NI, containing the abscissas of C the desired points. C C The output argument is C YI = array of dimension NI, where the ordinates of the C desired points are to be stored. C C If an integer value smaller than 3 is given to the NP argument, C this subroutine assumes NP = 3. C C The XI array elements need not be monotonic, but this C subroutine interpolates faster if the XI array elements are C given in a monotonic order. C C If the XI array element is less than XD(1) or greater than C XD(ND), this subroutine linearly interpolates the YI value. C C C Specification statement DIMENSION XD(*),YD(*),XI(*), YI(*) C Error check 10 IF (ND.LE.1) GO TO 90 IF (NI.LE.0) GO TO 91 DO 11 ID=2,ND IF (XD(ID).LE.XD(ID-1)) GO TO 92 11 CONTINUE C Branches off special cases. IF (ND.LE.4) GO TO 50 C General case -- Five data points of more C Calculates some local variables. 20 NP0=MAX(3,NP) NPM1=NP0-1 RENPM1=NPM1 RENNM2=NP0*(NP0-2) C Main calculation for the general case C First (outermost) DO-loop with respect to the desired points 30 DO 39 II=1,NI IF (II.EQ.1) IINTPV=-1 XII=XI(II) C Locates the interval that includes the desired point by binary C search. IF (XII.LE.XD(1)) THEN IINT=0 ELSE IF (XII.LT.XD(ND)) THEN IDMN=1 IDMX=ND IDMD=(IDMN+IDMX)/2 31 IF (XII.GE.XD(IDMD)) THEN IDMN=IDMD ELSE IDMX=IDMD END IF IDMD=(IDMN+IDMX)/2 IF (IDMD.GT.IDMN) GO TO 31 IINT=IDMD ELSE IINT=ND END IF C End of locating the interval of interest C Interpolation or extrapolation in one of the three subcases IF (IINT.LE.0) THEN C Subcase 1 -- Linear extrapolation when the abscissa of the C desired point is equal to that of the first data C point or less. C Estimates the first derivative when the interval is not the C same as the one for the previous desired point. -- C cf. Equation (8) IF (IINT.NE.IINTPV) THEN IINTPV=IINT X0=XD(1) X1=XD(2)-X0 X2=XD(3)-X0 X3=XD(4)-X0 Y0=YD(1) Y1=YD(2)-Y0 Y2=YD(3)-Y0 Y3=YD(4)-Y0 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) A1=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT END IF C Evaluates the YI value. YI(II)=Y0+A1*(XII-X0) C End of Subcase 1 ELSE IF (IINT.GE.ND) THEN C Subcase 2 -- Linear extrapolation when the abscissa of the C desired point is equal to that of the last data C point or greater. C Estimates the first derivative when the interval is not the C same as the one for the previous desired point. -- C cf. Equation (8) IF (IINT.NE.IINTPV) THEN IINTPV=IINT X0=XD(ND) X1=XD(ND-1)-X0 X2=XD(ND-2)-X0 X3=XD(ND-3)-X0 Y0=YD(ND) Y1=YD(ND-1)-Y0 Y2=YD(ND-2)-Y0 Y3=YD(ND-3)-Y0 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) A1=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT END IF C Evaluates the YI value. YI(II)=Y0+A1*(XII-X0) C End of Subcase 2 ELSE C Subcase 3 -- Interpolation when the abscissa of the desired C point is between those of the first and last C data points. C Calculates the coefficients of the third-degree polynomial (for C NP.LE.3) or the factors for the higher-degree polynomials (for C NP.GT.3), when the interval is not the same as the one for the C previous desired point. IF (IINT.NE.IINTPV) THEN IINTPV=IINT C The second DO-loop with respect to the two endpoints of the C interval DO 37 IEPT=1,2 C Calculates the estimate of the first derivative at an endpoint. C Initial setting for calculation ID0=IINT+IEPT-1 X0=XD(ID0) Y0=YD(ID0) SMPEF=0.0 SMWTF=0.0 SMPEI=0.0 SMWTI=0.0 C The third (innermost) DO-loop with respect to the four primary C estimate of the first derivative DO 36 IPE=1,4 C Selects point numbers of four consecutive data points for C calculating the primary estimate of the first derivative. IF (IPE.EQ.1) THEN ID1=ID0-3 ID2=ID0-2 ID3=ID0-1 ELSE IF (IPE.EQ.2) THEN ID1=ID0+1 ELSE IF (IPE.EQ.3) THEN ID2=ID0+2 ELSE ID3=ID0+3 END IF C Checks if any point number falls outside the legitimate range C (between 1 and ND). Skips calculation of the primary estimate C if any does. IF (ID1.LT.1.OR.ID2.LT.1.OR.ID3.LT.1.OR. 1 ID1.GT.ND.OR.ID2.GT.ND.OR.ID3.GT.ND) 2 GO TO 36 C Calculates the primary estimate of the first derivative -- C cf. Equation (8) X1=XD(ID1)-X0 X2=XD(ID2)-X0 X3=XD(ID3)-X0 Y1=YD(ID1)-Y0 Y2=YD(ID2)-Y0 Y3=YD(ID3)-Y0 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) PE=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT C Calculates the volatility factor, VOL, and distance factor, C SXX, for the primary estimate. -- cf. Equations (9) and (11) SX=X1+X2+X3 SY=Y1+Y2+Y3 SXX=X1*X1+X2*X2+X3*X3 SXY=X1*Y1+X2*Y2+X3*Y3 DNM=4.0*SXX-SX*SX B0=(SXX*SY-SX*SXY)/DNM B1=(4.0*SXY-SX*SY)/DNM DY0=-B0 DY1=Y1-(B0+B1*X1) DY2=Y2-(B0+B1*X2) DY3=Y3-(B0+B1*X3) VOL=DY0*DY0+DY1*DY1+DY2*DY2+DY3*DY3 C Calculates the EPSLN value, which is used to decide whether or C not the volatility factor, VOL, is essentially zero. EPSLN=(YD(ID0)**2+YD(ID1)**2 1 +YD(ID2)**2+YD(ID3)**2)*1.0E-12 C Accumulates the weighted primary estimates. -- C cf. Equations (13) and (14) IF (VOL.GT.EPSLN) THEN C - For finite weight. WT=1.0/(VOL*SXX) SMPEF=SMPEF+PE*WT SMWTF=SMWTF+WT ELSE C - For infinite weight. SMPEI=SMPEI+PE SMWTI=SMWTI+1.0 END IF 36 CONTINUE C End of the third DO-loop C Calculates the final estimate of the first derivative. -- C cf. Equation (14) IF (SMWTI.LT.0.5) THEN C - When no infinite weights exist. YP=SMPEF/SMWTF ELSE C - When infinite weights exist. YP=SMPEI/SMWTI END IF IF (IEPT.EQ.1) THEN YP0=YP ELSE YP1=YP END IF C End of the calculation of the estimate of the first derivative C at an endpoint 37 CONTINUE C End of the second DO-loop IF (NP0.LE.3) THEN C Calculates the coefficients of the third-degree polynomial C (when NP.LE.3). -- cf. Equation (4) DX=XD(IINT+1)-XD(IINT) DY=YD(IINT+1)-YD(IINT) A0=YD(IINT) A1=YP0 YP1=YP1-YP0 YP0=YP0-DY/DX A2=-(3.0*YP0+YP1)/DX A3= (2.0*YP0+YP1)/(DX*DX) ELSE C Calculates the factors for the higher-degree polynomials C (when NP.GT.3). -- cf. Equation (20) DX=XD(IINT+1)-XD(IINT) DY=YD(IINT+1)-YD(IINT) T0=YP0*DX-DY T1=YP1*DX-DY AA0= (T0+RENPM1*T1)/RENNM2 AA1=-(RENPM1*T0+T1)/RENNM2 END IF END IF C End of the calculation of the coefficients of the third-degree C polynomial (when NP.LE.3) or the factors for the higher-degree C polynomials (when NP.GT.3), when the interval is not the same C as the one for the previous desired point. C Evaluates the YI value. IF (NP0.LE.3) THEN C - With a third-degree polynomial. -- cf. Equation (3) XX=XII-XD(IINT) YI(II)=A0+XX*(A1+XX*(A2+XX*A3)) ELSE C - With a higher-degree polynomial. -- cf. Equation (19) U=(XII-XD(IINT))/DX UC=1.0-U V=AA0*((U**NP0)-U)+AA1*((UC**NP0)-UC) YI(II)=YD(IINT)+DY*U+V END IF C End of Subcase 3 END IF 39 CONTINUE C End of the first DO-loop C End of general case RETURN C Special cases -- Four data points or less C Preliminary processing for the special cases 50 X0=XD(1) Y0=YD(1) X1=XD(2)-X0 Y1=YD(2)-Y0 IF (ND.EQ.2) GO TO 60 X2=XD(3)-X0 Y2=YD(3)-Y0 IF (ND.EQ.3) GO TO 70 X3=XD(4)-X0 Y3=YD(4)-Y0 GO TO 80 C Special Case 1 -- Two data points C (Linear interpolation and extrapolation) 60 A1=Y1/X1 DO 61 II=1,NI YI(II)=Y0+A1*(XI(II)-X0) 61 CONTINUE C End of Special Case 1 RETURN C Special Case 2 -- Three data points C (Quadratic interpolation and linear extrapolation) 70 DLT=X1*X2*(X2-X1) A1=(X2*X2*Y1-X1*X1*Y2)/DLT A2=(X1*Y2-X2*Y1)/DLT A12=2.0*A2*X2+A1 DO 71 II=1,NI XX=XI(II)-X0 IF (XX.LE.0.0) THEN YI(II)=Y0+A1*XX ELSE IF (XX.LT.X2) THEN YI(II)=Y0+XX*(A1+XX*A2) ELSE YI(II)=Y0+Y2+A12*(XX-X2) END IF 71 CONTINUE C End of Special Case 2 RETURN C Special Case 3 -- Four data points C (Cubic interpolation and linear extrapolation) 80 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) A1=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT A2=(X2*X3*(X2*X2-X3*X3)*Y1 1 +X3*X1*(X3*X3-X1*X1)*Y2 2 +X1*X2*(X1*X1-X2*X2)*Y3)/DLT A3=(X2*X3*(X3-X2)*Y1 1 +X3*X1*(X1-X3)*Y2 2 +X1*X2*(X2-X1)*Y3)/DLT A13=(3.0*A3*X3+2.0*A2)*X3+A1 DO 81 II=1,NI XX=XI(II)-X0 IF (XX.LE.0.0) THEN YI(II)=Y0+A1*XX ELSE IF (XX.LT.X3) THEN YI(II)=Y0+XX*(A1+XX*(A2+XX*A3)) ELSE YI(II)=Y0+Y3+A13*(XX-X3) END IF 81 CONTINUE C End of Special Case 3 RETURN C Error exit 90 WRITE (*,99090) ND GO TO 99 91 WRITE (*,99091) NI GO TO 99 92 WRITE (*,99092) ID,XD(ID-1),XD(ID) 99 WRITE (*,99099) RETURN C Format statements for error messages 99090 FORMAT (1X/ ' *** Insufficient data points.' 1 7X,'ND =',I3) 99091 FORMAT (1X/ ' *** No desired points.' 1 7X,'NI =',I3) 99092 FORMAT (1X/ ' *** Two data points identical or out of ', 1 'sequence.'/ 2 7X,'ID, XD(ID-1), XD(ID) =',I5,2F10.3) 99099 FORMAT (' Error detected in the UVIP3P subroutine'/) END PROGRAM TTUVIP C C Test Program for the UVIP3P Subroutine C C Hiroshi Akima C U.S. Department of Commerce, NTIA/ITS C Version of 89/07/04 C C This program tests the UVIP3P subroutine by calling UVIP3P, in C a sequence, with two values of the NP argument, which is the C degree of the polynomial for the interval between each pair of C data points. With each NP value, this program calls UVIP3P in C three ways: C (1) in repeated calls with interpolation at a point in C each call, C (2) in one call with interpolation at all points, and C (3) in one call as in (2) but with inverted data set. C The program then compares the results with the expected values, C precalculated and stored in the program, and writes the C differences in two pages. If all entries in the last three C columns on each page are all zeros, UVIP3P is considered to be C working as expected. Otherwise, something is wrong in UVIP3P C and/or the test program itself. C C Dummy arguments of the UVIP3P subroutine are described in the C comments in the subroutine. C C Variables and arrays used in this program are C NP = degree of the polynomial for the interval between C each pair of data points, C NPO = integer array of dimension 2 containing the two C NP values, C ND = number of data points, C NI = number of interpolated points, C XD1, YD1 C = arrays of dimension ND containing the abscissas C and ordinates of the data points, C XI1 = array of dimension NI containing the abscissas of C the desired points, C YIE = doubly-dimensioned array of dimension (NI,2) C containing the expected values of the interpolated C results for the two NP values in two columns, C YI1 = array of dimension NI, where the ordinate values C interpolated from the XD1, YD1, and XI1 values in C repeated calls are to be stored, C YI2 = array of dimension NI, where the ordinate values C interpolated from the XD1, YD1, and XI1 values in C one call are to be stored, C XD3, YD3 C = arrays of dimension ND, where the abscissas and C ordinates of the inverted data points are to be C stored, C XI3 = array of dimension NI, where the abscissas of the C desired points for the inverted data points are to C be stored, C YI3 = array of dimension NI, where the ordinate values C interpolated from the XD3, YD3, and XI3 values in C one call are to be stored, C DYI1, DYI2, DYI3 C = arrays of dimension NI, where the differences C between the YI1, YI2, and YI3 values and their C respective expected values are to be stored. C C C Specification statements DIMENSION NPO(2) PARAMETER (ND=10, NI=31) DIMENSION XD1(ND),YD1(ND),YIE(NI,2),XI1(NI) DIMENSION XD3(ND),YD3(ND),XI3(NI) DIMENSION YI1(NI),YI2(NI),YI3(NI) DIMENSION DYI1(NI),DYI2(NI),DYI3(NI) C Data statements DATA NPO/ 3, 6/ DATA XD1/ 1 1.0, 2.0, 4.0, 6.5, 8.0,10.0,10.5,11.0,13.0,14.0/ DATA YD1/ 1 0.0, 0.0, 0.0, 0.0, 0.1, 1.0, 4.5, 8.0,10.0,15.0/ DATA XI1/ 1 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 2 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 3 8.0, 8.5, 9.0, 9.5,10.0,10.5,11.0,11.5, 4 12.0,12.5,13.0,13.5,14.0,14.5,15.0/ DATA YIE/ 1 8*0.0, 2 6*0.0, 0.015, 0.052, 3 0.100, 0.036,-0.045, 0.172, 1.000, 4.500, 8.000,10.075, 4 10.705,10.483,10.000,11.204,15.000,19.767,24.533, 1 8*0.0, 2 6*0.0, 0.020, 0.057, 3 0.100, 0.134, 0.166, 0.314, 1.000, 4.500, 8.000, 9.689, 4 10.101,10.180,10.000,11.663,15.000,19.767,24.533/ C Opens the file (for the UNIX operating system). 10 OPEN (6,FILE='TTUOUT',FORM='PRINT') C Calculation C Do-loop with respect to the two NP values 20 DO 79 INP=1,2 NP=NPO(INP) C Interpolation in repeated calls at a point in each call 30 DO 31 II=1,NI CALL UVIP3P(NP,ND,XD1,YD1,1,XI1(II),YI1(II)) 31 CONTINUE C Interpolation in one call at all points CALL UVIP3P(NP,ND,XD1,YD1,NI,XI1,YI2) C Interpolation in one call with inverted data 40 DO 41 ID=1,ND IDI=ND+1-ID XD3(IDI)=XD1(1)+XD1(ND)-XD1(ID) YD3(IDI)=YD1(ID) 41 CONTINUE DO 42 II=1,NI XI3(II)=XD1(1)+XD1(ND)-XI1(II) 42 CONTINUE CALL UVIP3P(NP,ND,XD3,YD3,NI,XI3,YI3) C Calculation of the differences between the calculated and C expected values 50 DO 51 II=1,NI DYI1(II)=YI1(II)-YIE(II,INP) DYI2(II)=YI2(II)-YIE(II,INP) DYI3(II)=YI3(II)-YIE(II,INP) 51 CONTINUE C Writing the results 60 IF (INP.EQ.1) THEN WRITE (*,6060) NP ELSE WRITE (*,6061) NP END IF WRITE (*,6062) DO 69 II=1,NI IF (MOD(II,4).EQ.2) WRITE (*,6063) IF (II.LE.ND) THEN ID=II WRITE (*,6064) XD1(ID),YD1(ID), 1 XI1(II),YIE(II,INP),YI1(II), 2 DYI1(II),DYI2(II),DYI3(II) ELSE WRITE (*,6065) XI1(II),YIE(II,INP),YI1(II), 1 DYI1(II),DYI2(II),DYI3(II) END IF 69 CONTINUE 79 CONTINUE STOP C Format statements 6060 FORMAT (1H , 1 'TTUVIP',14X,'Test Program for UVIP3P',29X,'NP =',I2) 6061 FORMAT (1H1, 1 'TTUVIP',14X,'Test Program for UVIP3P',29X,'NP =',I2) 6062 FORMAT (1X/// 1 6X,'Data Points',9X,'Interpolated Points ', 2 8X,' Differences '// 3 6X,'XD1 YD1',8X,'XI1 YIE YI1', 4 8X,'DYI1 DYI2 DYI3'/) 6063 FORMAT (1X) 6064 FORMAT (1X,2F9.3,2X,3F9.3,2X,3F9.3) 6065 FORMAT (19X,2X,3F9.3,2X,3F9.3) END SUBROUTINE UVIP3P(NP,ND,XD,YD,NI,XI, YI) C C Univariate Interpolation (Improved Akima Method) C C Hiroshi Akima C U.S. Department of Commerce, NTIA/ITS C Version of 89/07/04 C C This subroutine performs univariate interpolation. It is based C on the improved A method developed by Hiroshi Akima, 'A method C of univariate interpolation that has the accuracy of a third- C degree polynomial,' ACM TOMS, vol. xx, pp. xxx-xxx, 19xx. (The C equation numbers referred to in the comments below are those in C the paper.) C C In this method, the interpolating function is a piecewise C function composed of a set of polynomials applicable to C successive intervals of the given data points. This method C uses third-degree polynomials as the default, but the user has C an option to use higher-degree polynomial to reduce undulations C in resulting curves. C C This method has the accuracy of a third-degree polynomial if C the degree of the polynomials for the interpolating function is C set to three. C C The input arguments are C NP = degree of the polynomials for the interpolating C function, C ND = number of input data points C (must be equal to 2 or greater), C XD = array of dimension ND, containing the abscissas of C the input data points C (must be in a monotonic increasing order), C YD = array of dimension ND, containing the ordinates of C the input data points, C NI = number of points for which interpolation is desired C (must be equal to 1 or greater), C XI = array of dimension NI, containing the abscissas of C the desired points. C C The output argument is C YI = array of dimension NI, where the ordinates of the C desired points are to be stored. C C If an integer value smaller than 3 is given to the NP argument, C this subroutine assumes NP = 3. C C The XI array elements need not be monotonic, but this C subroutine interpolates faster if the XI array elements are C given in a monotonic order. C C If the XI array element is less than XD(1) or greater than C XD(ND), this subroutine linearly interpolates the YI value. C C C Specification statement DIMENSION XD(*),YD(*),XI(*), YI(*) C Error check 10 IF (ND.LE.1) GO TO 90 IF (NI.LE.0) GO TO 91 DO 11 ID=2,ND IF (XD(ID).LE.XD(ID-1)) GO TO 92 11 CONTINUE C Branches off special cases. IF (ND.LE.4) GO TO 50 C General case -- Five data points of more C Calculates some local variables. 20 NP0=MAX(3,NP) NPM1=NP0-1 RENPM1=NPM1 RENNM2=NP0*(NP0-2) C Main calculation for the general case C First (outermost) DO-loop with respect to the desired points 30 DO 39 II=1,NI IF (II.EQ.1) IINTPV=-1 XII=XI(II) C Locates the interval that includes the desired point by binary C search. IF (XII.LE.XD(1)) THEN IINT=0 ELSE IF (XII.LT.XD(ND)) THEN IDMN=1 IDMX=ND IDMD=(IDMN+IDMX)/2 31 IF (XII.GE.XD(IDMD)) THEN IDMN=IDMD ELSE IDMX=IDMD END IF IDMD=(IDMN+IDMX)/2 IF (IDMD.GT.IDMN) GO TO 31 IINT=IDMD ELSE IINT=ND END IF C End of locating the interval of interest C Interpolation or extrapolation in one of the three subcases IF (IINT.LE.0) THEN C Subcase 1 -- Linear extrapolation when the abscissa of the C desired point is equal to that of the first data C point or less. C Estimates the first derivative when the interval is not the C same as the one for the previous desired point. -- C cf. Equation (8) IF (IINT.NE.IINTPV) THEN IINTPV=IINT X0=XD(1) X1=XD(2)-X0 X2=XD(3)-X0 X3=XD(4)-X0 Y0=YD(1) Y1=YD(2)-Y0 Y2=YD(3)-Y0 Y3=YD(4)-Y0 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) A1=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT END IF C Evaluates the YI value. YI(II)=Y0+A1*(XII-X0) C End of Subcase 1 ELSE IF (IINT.GE.ND) THEN C Subcase 2 -- Linear extrapolation when the abscissa of the C desired point is equal to that of the last data C point or greater. C Estimates the first derivative when the interval is not the C same as the one for the previous desired point. -- C cf. Equation (8) IF (IINT.NE.IINTPV) THEN IINTPV=IINT X0=XD(ND) X1=XD(ND-1)-X0 X2=XD(ND-2)-X0 X3=XD(ND-3)-X0 Y0=YD(ND) Y1=YD(ND-1)-Y0 Y2=YD(ND-2)-Y0 Y3=YD(ND-3)-Y0 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) A1=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT END IF C Evaluates the YI value. YI(II)=Y0+A1*(XII-X0) C End of Subcase 2 ELSE C Subcase 3 -- Interpolation when the abscissa of the desired C point is between those of the first and last C data points. C Calculates the coefficients of the third-degree polynomial (for C NP.LE.3) or the factors for the higher-degree polynomials (for C NP.GT.3), when the interval is not the same as the one for the C previous desired point. IF (IINT.NE.IINTPV) THEN IINTPV=IINT C The second DO-loop with respect to the two endpoints of the C interval DO 37 IEPT=1,2 C Calculates the estimate of the first derivative at an endpoint. C Initial setting for calculation ID0=IINT+IEPT-1 X0=XD(ID0) Y0=YD(ID0) SMPEF=0.0 SMWTF=0.0 SMPEI=0.0 SMWTI=0.0 C The third (innermost) DO-loop with respect to the four primary C estimate of the first derivative DO 36 IPE=1,4 C Selects point numbers of four consecutive data points for C calculating the primary estimate of the first derivative. IF (IPE.EQ.1) THEN ID1=ID0-3 ID2=ID0-2 ID3=ID0-1 ELSE IF (IPE.EQ.2) THEN ID1=ID0+1 ELSE IF (IPE.EQ.3) THEN ID2=ID0+2 ELSE ID3=ID0+3 END IF C Checks if any point number falls outside the legitimate range C (between 1 and ND). Skips calculation of the primary estimate C if any does. IF (ID1.LT.1.OR.ID2.LT.1.OR.ID3.LT.1.OR. 1 ID1.GT.ND.OR.ID2.GT.ND.OR.ID3.GT.ND) 2 GO TO 36 C Calculates the primary estimate of the first derivative -- C cf. Equation (8) X1=XD(ID1)-X0 X2=XD(ID2)-X0 X3=XD(ID3)-X0 Y1=YD(ID1)-Y0 Y2=YD(ID2)-Y0 Y3=YD(ID3)-Y0 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) PE=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT C Calculates the volatility factor, VOL, and distance factor, C SXX, for the primary estimate. -- cf. Equations (9) and (11) SX=X1+X2+X3 SY=Y1+Y2+Y3 SXX=X1*X1+X2*X2+X3*X3 SXY=X1*Y1+X2*Y2+X3*Y3 DNM=4.0*SXX-SX*SX B0=(SXX*SY-SX*SXY)/DNM B1=(4.0*SXY-SX*SY)/DNM DY0=-B0 DY1=Y1-(B0+B1*X1) DY2=Y2-(B0+B1*X2) DY3=Y3-(B0+B1*X3) VOL=DY0*DY0+DY1*DY1+DY2*DY2+DY3*DY3 C Calculates the EPSLN value, which is used to decide whether or C not the volatility factor, VOL, is essentially zero. EPSLN=(YD(ID0)**2+YD(ID1)**2 1 +YD(ID2)**2+YD(ID3)**2)*1.0E-12 C Accumulates the weighted primary estimates. -- C cf. Equations (13) and (14) IF (VOL.GT.EPSLN) THEN C - For finite weight. WT=1.0/(VOL*SXX) SMPEF=SMPEF+PE*WT SMWTF=SMWTF+WT ELSE C - For infinite weight. SMPEI=SMPEI+PE SMWTI=SMWTI+1.0 END IF 36 CONTINUE C End of the third DO-loop C Calculates the final estimate of the first derivative. -- C cf. Equation (14) IF (SMWTI.LT.0.5) THEN C - When no infinite weights exist. YP=SMPEF/SMWTF ELSE C - When infinite weights exist. YP=SMPEI/SMWTI END IF IF (IEPT.EQ.1) THEN YP0=YP ELSE YP1=YP END IF C End of the calculation of the estimate of the first derivative C at an endpoint 37 CONTINUE C End of the second DO-loop IF (NP0.LE.3) THEN C Calculates the coefficients of the third-degree polynomial C (when NP.LE.3). -- cf. Equation (4) DX=XD(IINT+1)-XD(IINT) DY=YD(IINT+1)-YD(IINT) A0=YD(IINT) A1=YP0 YP1=YP1-YP0 YP0=YP0-DY/DX A2=-(3.0*YP0+YP1)/DX A3= (2.0*YP0+YP1)/(DX*DX) ELSE C Calculates the factors for the higher-degree polynomials C (when NP.GT.3). -- cf. Equation (20) DX=XD(IINT+1)-XD(IINT) DY=YD(IINT+1)-YD(IINT) T0=YP0*DX-DY T1=YP1*DX-DY AA0= (T0+RENPM1*T1)/RENNM2 AA1=-(RENPM1*T0+T1)/RENNM2 END IF END IF C End of the calculation of the coefficients of the third-degree C polynomial (when NP.LE.3) or the factors for the higher-degree C polynomials (when NP.GT.3), when the interval is not the same C as the one for the previous desired point. C Evaluates the YI value. IF (NP0.LE.3) THEN C - With a third-degree polynomial. -- cf. Equation (3) XX=XII-XD(IINT) YI(II)=A0+XX*(A1+XX*(A2+XX*A3)) ELSE C - With a higher-degree polynomial. -- cf. Equation (19) U=(XII-XD(IINT))/DX UC=1.0-U V=AA0*((U**NP0)-U)+AA1*((UC**NP0)-UC) YI(II)=YD(IINT)+DY*U+V END IF C End of Subcase 3 END IF 39 CONTINUE C End of the first DO-loop C End of general case RETURN C Special cases -- Four data points or less C Preliminary processing for the special cases 50 X0=XD(1) Y0=YD(1) X1=XD(2)-X0 Y1=YD(2)-Y0 IF (ND.EQ.2) GO TO 60 X2=XD(3)-X0 Y2=YD(3)-Y0 IF (ND.EQ.3) GO TO 70 X3=XD(4)-X0 Y3=YD(4)-Y0 GO TO 80 C Special Case 1 -- Two data points C (Linear interpolation and extrapolation) 60 A1=Y1/X1 DO 61 II=1,NI YI(II)=Y0+A1*(XI(II)-X0) 61 CONTINUE C End of Special Case 1 RETURN C Special Case 2 -- Three data points C (Quadratic interpolation and linear extrapolation) 70 DLT=X1*X2*(X2-X1) A1=(X2*X2*Y1-X1*X1*Y2)/DLT A2=(X1*Y2-X2*Y1)/DLT A12=2.0*A2*X2+A1 DO 71 II=1,NI XX=XI(II)-X0 IF (XX.LE.0.0) THEN YI(II)=Y0+A1*XX ELSE IF (XX.LT.X2) THEN YI(II)=Y0+XX*(A1+XX*A2) ELSE YI(II)=Y0+Y2+A12*(XX-X2) END IF 71 CONTINUE C End of Special Case 2 RETURN C Special Case 3 -- Four data points C (Cubic interpolation and linear extrapolation) 80 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1) A1=(((X2*X3)**2)*(X3-X2)*Y1 1 +((X3*X1)**2)*(X1-X3)*Y2 2 +((X1*X2)**2)*(X2-X1)*Y3)/DLT A2=(X2*X3*(X2*X2-X3*X3)*Y1 1 +X3*X1*(X3*X3-X1*X1)*Y2 2 +X1*X2*(X1*X1-X2*X2)*Y3)/DLT A3=(X2*X3*(X3-X2)*Y1 1 +X3*X1*(X1-X3)*Y2 2 +X1*X2*(X2-X1)*Y3)/DLT A13=(3.0*A3*X3+2.0*A2)*X3+A1 DO 81 II=1,NI XX=XI(II)-X0 IF (XX.LE.0.0) THEN YI(II)=Y0+A1*XX ELSE IF (XX.LT.X3) THEN YI(II)=Y0+XX*(A1+XX*(A2+XX*A3)) ELSE YI(II)=Y0+Y3+A13*(XX-X3) END IF 81 CONTINUE C End of Special Case 3 RETURN C Error exit 90 WRITE (*,99090) ND GO TO 99 91 WRITE (*,99091) NI GO TO 99 92 WRITE (*,99092) ID,XD(ID-1),XD(ID) 99 WRITE (*,99099) RETURN C Format statements for error messages 99090 FORMAT (1X/ ' *** Insufficient data points.' 1 7X,'ND =',I3) 99091 FORMAT (1X/ ' *** No desired points.' 1 7X,'NI =',I3) 99092 FORMAT (1X/ ' *** Two data points identical or out of ', 1 'sequence.'/ 2 7X,'ID, XD(ID-1), XD(ID) =',I5,2F10.3) 99099 FORMAT (' Error detected in the UVIP3P subroutine'/) END