1      SUBROUTINE ASMINT (NPT,X,Y,NVAL,XVAL,YVAL)
2C
3C***  NON-LINEAR INTERPOLATION ROUTINE FOR AIRFOIL MODULE
4C
5      INTEGER RP
6      DIMENSION X(1),Y(1),XVAL(1),YVAL(1),YP(2)
7C
8C   METHOD - SLOPES ARE DEFINED AT EACH DATA POINT AS A VALUE HALF WAY
9C            BETWEEN THE LINEAR SLOPES TO THE RIGHT AND LEFT OF THE
10C            POINT. A THIRD ORDER POLYNOMIAL IS DETERMINED FOR EACH
11C            INTERVAL BETWEEN DATA POINTS BY USING THE TWO END POINTS
12C            AND TWO END POINT SLOPES. SECOND ORDER POLYNOMIALS ARE
13C            DETERMINED FOR THE LEFT AND RIGHT END INTERVALS (AND FOR
14C            LEFT AND RIGHT EXTRAPOLATION) USING THE TWO END POINTS
15C            AT EACH END OF THE DATA AND THE SLOPES AT THE NEXT-TO-THE
16C            LAST DATA POINTS AT EACH END OF THE DATA. HENCE, THE OVER-
17C            ALL CURVE CONSISTS OF A LEFT END PARABOLA JOINED TO A
18C            SERIES OF JOINED CUBIC CURVES FINALLY JOINED TO A RIGHT
19C            END PARABOLA. THE OVERALL CURVE IS CONTINUOUS EVERYWHERE,
20C            HAS CONTINUOUS DERIVATIVES EVERYWHERE, AND CAN BE EXTRAP-
21C            OLATED TO THE RIGHT AND LEFT INDEFINITELY.
22C
23C   NPT      NUMBER OF POINTS IN X AND Y ARRAYS
24C   X        ABSCISSAS OF INPUT DATA
25C   Y        ORDINATES OF INPUT DATA
26C
27C   NVAL     NUMBER OF POINTS TO BE INTERPOLATED
28C   XVAL     ABSCISSA VALUES AT WHICH INTERPOLATION TO BE PERFORMED
29C   YVAL     INTERPOLATED ORDINATE VALUES
30C
31      NP1= NPT-1
32      NP2= NPT-2
33      DO 1080 I=1,NVAL
34         IF(XVAL(I).LT.X(2)) GO TO 1020
35         IF(XVAL(I).GT.X(NP1)) GO TO 1030
36            DO 1010 J=2,NP2
37            IF(XVAL(I).NE.X(J)) GO TO 1000
38            YVAL(I)= Y(J)
39            GO TO 1080
40 1000       K= J+1
41            IF(XVAL(I).GE.X(K))GO TO 1010
42            LOCATE= 2
43            GO TO 1040
44 1010    CONTINUE
45         YVAL(I)= Y(K)
46         GO TO 1080
47 1020    J=2
48         K=3
49         LOCATE= 1
50         GO TO 1040
51 1030    J=NP1
52         K=NPT
53         LOCATE= 3
54         GO TO 1040
55 1040    DO 1050 N=1,2
56            LP= J+N-2
57            MP= J+N-1
58            RP= K+N-1
59            SL= (Y(LP)-Y(MP))/(X(LP)-X(MP))
60            SR= (Y(MP)-Y(RP))/(X(MP)-X(RP))
61            ANGL= ATAN(SL)
62            ANGR= ATAN(SR)
63            ANGAV=(ANGL+ANGR)/2.0
64            YP(N)= SIN(ANGAV)/COS(ANGAV)
65            GO TO (1070,1050,1070),LOCATE
66 1050    CONTINUE
67C
68         X1S = X(J)*X(J)
69         X2S = X(K)*X(K)
70         X12F = X(J)-X(K)
71         Y12F = Y(J)-Y(K)
72         X12S = X1S-X2S
73         X12C = X1S*X(J)-X2S*X(K)
74         Y12P = YP(1)-YP(2)
75C
76         TW = 2.0
77         TH = 3.0
78         RED = TW*X(K)*X12F - X12S
79         GRN = TH*X2S*X12F - X12C
80         YEL = YP(2)*X12F-Y12F
81         E = TH*X12S*RED - TW*X12F*GRN
82         A = Y12P*RED - TW*X12F*YEL
83         A = A/E
84         B = TH*X12S*YEL - Y12P*GRN
85         B = B/E
86         C = (Y12F-A*X12C-B*X12S)/X12F
87         D = Y(K)-A*X2S*X(K)-B*X2S-C*X(K)
88C
89 1060    YVAL(I)=(((A*XVAL(I) +B)*XVAL(I))+C)*XVAL(I)+D
90         GO TO 1080
91C
92 1070    J=1
93         IF(LOCATE.EQ.3) J=NP1
94         K=J+1
95         L=2
96         IF(LOCATE.EQ.3) L=NP1
97         Z=(Y(J)-Y(K))/(X(J)-X(K))
98         A=0.0
99         B=(YP(1)-Z)/(2.0*X(L)-X(J)-X(K))
100         C=YP(1)-2.0*B*X(L)
101         D=Y(J)-((B*X(J)+C)*X(J))
102         GO TO 1060
103 1080 CONTINUE
104      RETURN
105      END
106