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