SUBROUTINE SPLN1 (N,X,Y,J,D,C,W) SPLN1 c----------------------------------------------------------------------- c include 'param_sz.h' c dimension x(nspl),y(nspl),d(2),c(nsplc),w(nsplc) C ------------------------------------------------------------------SPLN1 C OVER THE INTERVAL X(I) TO X(I+1), THE INTERPOLATING SPLN1 C POLYNOMIAL SPLN1 C Y=Y(I)+A(I)*Z+B(I)*Z**2+E(I)*Z**3 SPLN1 C WHERE Z=(X-X(I))/(X(I+1)-X(I)) SPLN1 C IS USED. THE COEFFICIENTS A(I),B(I) AND E(I) ARE COMPUTED SPLN1 C BY SPLN1 AND STORED IN LOCATIONS C(3*I-2),C(3*I-1) AND SPLN1 C C(3*I) RESPECTIVELY. SPLN1 C WHILE WORKING IN THE ITH INTERVAL,THE VARIABLE Q WILL SPLN1 C REPRESENT Q=X(I+1) - X(I), AND Y(I) WILL REPRESENT SPLN1 C Y(I+1)-Y(I) SPLN1 C ------------------------------------------------------------------SPLN1 C SPLN1 Q=X(2) - X(1) SPLN1 YI =Y(2) - Y(1) SPLN1 IF (J.EQ.2) GO TO 100 SPLN1 C ------------------------------------------------------------------SPLN1 C IF THE FIRST DERIVATIVE AT THE END POINTS IS GIVEN, SPLN1 C A(1) IS KNOWN, AND THE SECOND EQUATION BECOMES SPLN1 C MERELY B(1)+E(1)=YI - Q*D(1). SPLN1 C ------------------------------------------------------------------SPLN1 C(1)=Q*D(1) SPLN1 C(2)=1.0 SPLN1 W(2)=YI-C(1) SPLN1 GO TO 200 SPLN1 C ------------------------------------------------------------------SPLN1 C IF THE SECOND DERIVATIVE AT THE END POINTS IS GIVEN SPLN1 C B(1) IS KNOWN, THE SECOND EQUATION BECOMES SPLN1 C A(1)+E(1)=YI-0.5*Q*Q*D(1). DURING THE SOLUTION OF SPLN1 C THE 3N-4 EQUATIONS,A1 WILL BE KEPT IN CELL C(2) SPLN1 C INSTEAD OF C(1) TO RETAIN THE TRIDIAGONAL FORM OF THE SPLN1 C COEFFICIENT MATRIX. SPLN1 C ------------------------------------------------------------------SPLN1 100 C(2)=0.0 SPLN1 W(2)=0.5*Q*Q*D(1) SPLN1 200 M=N-2 SPLN1 IF(M.LE.0) GO TO 350 SPLN1 C ------------------------------------------------------------------SPLN1 C UPPER TRIANGULARIZATION OF THE TRIDIAGONAL SYSTEM OF SPLN1 C EQUATIONS FOR THE COEFFICIENT MATRIX FOLLOWS-- SPLN1 C ------------------------------------------------------------------SPLN1 DO 300 I=1,M SPLN1 AI=Q SPLN1 Q=X(I+2)- X(I+1) SPLN1 H=AI/Q SPLN1 C(3*I)=-H/(2.0-C(3*I-1)) SPLN1 W(3*I)=(-YI-W(3*I-1))/(2.0 - C(3*I-1)) SPLN1 C(3*I+1)=-H*H/(H-C(3*I)) SPLN1 W(3*I+1)=(YI-W(3*I))/(H-C(3*I)) SPLN1 YI=Y(I+2)- Y(I+1) SPLN1 C(3*I+2)=1.0/(1.0-C(3*I+1)) SPLN1 300 W(3*I+2)=(YI-W(3*I+1))/(1.0-C(3*I+1)) SPLN1 C ------------------------------------------------------------------SPLN1 C E(N-1) IS DETERMINED DIRECTLY FROM THE LAST EQUATION SPLN1 C OBTAINED ABOVE, AND THE FIRST OR SECOND DERIVATIVE SPLN1 C VALUE GIVEN AT THE END POINT. SPLN1 C ------------------------------------------------------------------SPLN1 350 IF(J.EQ.1) GO TO 400 SPLN1 C(3*N-3)=(Q*Q*D(2)/2.0-W(3*N-4))/(3.0- C(3*N-4)) SPLN1 GO TO 500 SPLN1 400 C(3*N-3)=(Q*D(2)-YI-W(3*N-4))/(2.0-C(3*N-4)) SPLN1 500 M=3*N-6 SPLN1 IF(M.LE.0) GO TO 700 SPLN1 C ------------------------------------------------------------------SPLN1 C BACK SOLUTION FOR ALL COEFFICENTS EXCEPT SPLN1 C A(1) AND B(1) FOLLOWS-- SPLN1 C ------------------------------------------------------------------SPLN1 DO 600 II=1,M SPLN1 I=M-II+3 SPLN1 600 C(I)=W(I)-C(I)*C(I+1) SPLN1 700 IF(J.EQ.1) GO TO 800 SPLN1 C ------------------------------------------------------------------SPLN1 C IF THE SECOND DERIVATIVE IS GIVEN AT THE END POINTS, SPLN1 C A(1) CAN NOW BE COMPUTED FROM THE KNOWN VALUES OF SPLN1 C B(1) AND E(1). THEN A(1) AND B(1) ARE PUT INTO THEIR SPLN1 C PROPER PLACES IN THE C ARRAY. SPLN1 C ------------------------------------------------------------------SPLN1 C(1)=Y(2) - Y(1)-W(2)-C(3) SPLN1 C(2)=W(2) SPLN1 RETURN SPLN1 800 C(2)=W(2)-C(3) SPLN1 RETURN SPLN1 END SPLN1 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++