source: PSPA/parmelaPSPA/trunk/spln1.f @ 465

Last change on this file since 465 was 12, checked in by lemeur, 12 years ago

parmela pspa initial

File size: 6.9 KB
Line 
1      SUBROUTINE SPLN1 (N,X,Y,J,D,C,W)                                  SPLN1   
2c-----------------------------------------------------------------------
3c
4      include 'param_sz.h'
5c
6      dimension x(nspl),y(nspl),d(2),c(nsplc),w(nsplc)
7C     ------------------------------------------------------------------SPLN1   
8C             OVER THE INTERVAL X(I) TO X(I+1), THE INTERPOLATING       SPLN1   
9C             POLYNOMIAL                                                SPLN1   
10C                  Y=Y(I)+A(I)*Z+B(I)*Z**2+E(I)*Z**3                    SPLN1   
11C             WHERE Z=(X-X(I))/(X(I+1)-X(I))                            SPLN1   
12C             IS USED. THE COEFFICIENTS A(I),B(I) AND E(I) ARE COMPUTED SPLN1   
13C             BY SPLN1 AND STORED IN LOCATIONS C(3*I-2),C(3*I-1) AND    SPLN1   
14C             C(3*I) RESPECTIVELY.                                      SPLN1   
15C             WHILE WORKING IN THE ITH INTERVAL,THE VARIABLE Q WILL     SPLN1   
16C             REPRESENT Q=X(I+1) - X(I), AND Y(I) WILL REPRESENT        SPLN1   
17C             Y(I+1)-Y(I)                                               SPLN1   
18C     ------------------------------------------------------------------SPLN1   
19C                                                                       SPLN1   
20      Q=X(2) - X(1)                                                     SPLN1   
21      YI =Y(2) - Y(1)                                                   SPLN1   
22      IF (J.EQ.2) GO TO 100                                             SPLN1   
23C     ------------------------------------------------------------------SPLN1   
24C             IF THE FIRST DERIVATIVE AT THE END POINTS IS GIVEN,       SPLN1   
25C             A(1) IS KNOWN, AND THE SECOND EQUATION BECOMES            SPLN1   
26C             MERELY B(1)+E(1)=YI - Q*D(1).                             SPLN1   
27C     ------------------------------------------------------------------SPLN1   
28      C(1)=Q*D(1)                                                       SPLN1   
29      C(2)=1.0                                                          SPLN1   
30      W(2)=YI-C(1)                                                      SPLN1   
31      GO TO 200                                                         SPLN1   
32C     ------------------------------------------------------------------SPLN1   
33C             IF THE SECOND DERIVATIVE AT THE END POINTS IS GIVEN       SPLN1   
34C             B(1) IS KNOWN, THE SECOND EQUATION BECOMES                SPLN1   
35C             A(1)+E(1)=YI-0.5*Q*Q*D(1). DURING THE SOLUTION OF         SPLN1   
36C             THE 3N-4 EQUATIONS,A1 WILL BE KEPT IN CELL C(2)           SPLN1   
37C             INSTEAD OF C(1) TO RETAIN THE TRIDIAGONAL FORM OF THE     SPLN1   
38C             COEFFICIENT MATRIX.                                       SPLN1   
39C     ------------------------------------------------------------------SPLN1   
40100   C(2)=0.0                                                          SPLN1   
41      W(2)=0.5*Q*Q*D(1)                                                 SPLN1   
42200   M=N-2                                                             SPLN1   
43      IF(M.LE.0) GO TO 350                                              SPLN1   
44C     ------------------------------------------------------------------SPLN1   
45C             UPPER TRIANGULARIZATION OF THE TRIDIAGONAL SYSTEM OF      SPLN1   
46C             EQUATIONS FOR THE COEFFICIENT MATRIX FOLLOWS--            SPLN1   
47C     ------------------------------------------------------------------SPLN1   
48      DO 300 I=1,M                                                      SPLN1   
49      AI=Q                                                              SPLN1   
50      Q=X(I+2)- X(I+1)                                                  SPLN1   
51      H=AI/Q                                                            SPLN1   
52      C(3*I)=-H/(2.0-C(3*I-1))                                          SPLN1   
53      W(3*I)=(-YI-W(3*I-1))/(2.0 - C(3*I-1))                            SPLN1   
54      C(3*I+1)=-H*H/(H-C(3*I))                                          SPLN1   
55      W(3*I+1)=(YI-W(3*I))/(H-C(3*I))                                   SPLN1   
56      YI=Y(I+2)- Y(I+1)                                                 SPLN1   
57      C(3*I+2)=1.0/(1.0-C(3*I+1))                                       SPLN1   
58300   W(3*I+2)=(YI-W(3*I+1))/(1.0-C(3*I+1))                             SPLN1   
59C     ------------------------------------------------------------------SPLN1   
60C             E(N-1) IS DETERMINED DIRECTLY FROM THE LAST EQUATION      SPLN1   
61C             OBTAINED ABOVE, AND THE FIRST OR SECOND DERIVATIVE        SPLN1   
62C             VALUE GIVEN AT THE END POINT.                             SPLN1   
63C     ------------------------------------------------------------------SPLN1   
64350   IF(J.EQ.1) GO TO 400                                              SPLN1   
65      C(3*N-3)=(Q*Q*D(2)/2.0-W(3*N-4))/(3.0- C(3*N-4))                  SPLN1   
66      GO TO 500                                                         SPLN1   
67400   C(3*N-3)=(Q*D(2)-YI-W(3*N-4))/(2.0-C(3*N-4))                      SPLN1   
68500   M=3*N-6                                                           SPLN1   
69      IF(M.LE.0) GO TO 700                                              SPLN1   
70C     ------------------------------------------------------------------SPLN1   
71C             BACK SOLUTION FOR ALL COEFFICENTS EXCEPT                  SPLN1   
72C             A(1) AND B(1) FOLLOWS--                                   SPLN1   
73C     ------------------------------------------------------------------SPLN1   
74      DO 600 II=1,M                                                     SPLN1   
75      I=M-II+3                                                          SPLN1   
76600   C(I)=W(I)-C(I)*C(I+1)                                             SPLN1   
77700   IF(J.EQ.1) GO TO 800                                              SPLN1   
78C     ------------------------------------------------------------------SPLN1   
79C             IF THE SECOND DERIVATIVE IS GIVEN AT THE END POINTS,      SPLN1   
80C             A(1) CAN NOW BE COMPUTED FROM THE KNOWN VALUES OF         SPLN1   
81C             B(1) AND E(1). THEN A(1) AND B(1) ARE PUT INTO THEIR      SPLN1   
82C             PROPER PLACES IN THE C ARRAY.                             SPLN1   
83C     ------------------------------------------------------------------SPLN1   
84      C(1)=Y(2) - Y(1)-W(2)-C(3)                                        SPLN1   
85      C(2)=W(2)                                                         SPLN1   
86      RETURN                                                            SPLN1   
87800   C(2)=W(2)-C(3)                                                    SPLN1   
88      RETURN                                                            SPLN1   
89      END                                                               SPLN1   
90c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Note: See TracBrowser for help on using the repository browser.