1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
|
SUBROUTINE ALG01 (XDATA,YDATA,NDATA,XIN,YOUT,SLOPE,NXY,NTYPE,NWOT)
C
REAL M
C
DIMENSION XDATA(2),YDATA(2),XIN(1),YOUT(1),SLOPE(1),A(21),B(21),
1 D(21),M(21)
C
IF (NTYPE.EQ.1.OR.NDATA.LT.3) GO TO 210
A(1)=1.0
B(1)=0.0
D(1)=0.0
N=NDATA-1
DO 110 I=2,N
A(I)=(XDATA(I+1)-XDATA(I-1))/3.0-(XDATA(I)-XDATA(I-1))*B(I-1)/(6.0
1*A(I-1))
B(I)=(XDATA(I+1)-XDATA(I))/6.0
110 D(I)=(YDATA(I+1)-YDATA(I))/(XDATA(I+1)-XDATA(I))-(YDATA(I)-YDATA(I
1-1))/(XDATA(I)-XDATA(I-1))-(XDATA(I)-XDATA(I-1))*D(I-1)/(6.0*A(I-1
2))
A(NDATA)=0.0
B(NDATA)=1.0
D(NDATA)=0.0
M(NDATA)=A(NDATA)*D(N)/(A(NDATA)*B(N)-A(N)*B(NDATA))
DO 120 II=2,NDATA
I=NDATA+1-II
120 M(I)=(D(I)-B(I)*M(I+1))/A(I)
ASSIGN 150 TO IY
IF (NWOT.EQ.1) ASSIGN 160 TO IY
ASSIGN 160 TO ISLOPE
IF (NWOT.EQ.0) ASSIGN 200 TO ISLOPE
J=2
DO 200 I=1,NXY
IF (XIN(I).LT.XDATA(1)) GO TO 170
IF (XIN(I).GT.XDATA(NDATA)) GO TO 180
130 IF (XIN(I).LE.XDATA(J)) GO TO 140
J=J+1
GO TO 130
140 DX=XDATA(J)-XDATA(J-1)
GO TO IY, (150,160)
150 YOUT(I)=M(J-1)/(6.0*DX)*(XDATA(J)-XIN(I))**3+M(J)/(6.0*DX)*(XIN(I)
1-XDATA(J-1))**3+(XDATA(J)-XIN(I))*(YDATA(J-1)/DX-M(J-1)/6.0*DX)+(X
2IN(I)-XDATA(J-1))*(YDATA(J)/DX-M(J)/6.0*DX)
GO TO ISLOPE, (160,200)
160 SLOPE(I)=(-M(J-1)*(XDATA(J)-XIN(I))**2/2.0+M(J)*(XIN(I)-XDATA(J-1)
1)**2/2.0+YDATA(J)-YDATA(J-1))/DX-(M(J)-M(J-1))/6.0*DX
GO TO 200
170 JP=1
KP=2
GO TO 190
180 JP=NDATA
KP=N
190 YPRIME=(YDATA(KP)-YDATA(JP))/(XDATA(KP)-XDATA(JP))-M(KP)/6.0*(XDAT
1A(KP)-XDATA(JP))
IF (NWOT.NE.1) YOUT(I)=YDATA(JP)+(XIN(I)-XDATA(JP))*YPRIME
IF (NWOT.NE.0) SLOPE(I)=YPRIME
200 CONTINUE
RETURN
210 IF (NDATA.NE.1) GO TO 230
DO 220 I=1,NXY
220 YOUT(I)=YDATA(1)
RETURN
230 IF (NWOT.EQ.1) GO TO 254
J=2
DO 250 I=1,NXY
240 IF (XIN(I).LE.XDATA(J).OR.J.EQ.NDATA) GO TO 250
J=J+1
GO TO 240
250 YOUT(I)=YDATA(J-1)+(YDATA(J)-YDATA(J-1))/(XDATA(J)-XDATA(J-1))*(XI
1N(I)-XDATA(J-1))
IF (NWOT.NE.2) RETURN
254 YPRIME=(YDATA(2)-YDATA(1))/(XDATA(2)-XDATA(1))
DO 260 I=1,NXY
260 SLOPE(I)=YPRIME
RETURN
END
|