]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:26 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | FUNCTION DSPPS1(K,M,NDER,X,T,C,NERR) | |
11 | ||
12 | #include "gen/imp64.inc" | |
13 | DIMENSION T(*),C(*),B(27) | |
14 | CHARACTER NAME*(*) | |
15 | CHARACTER*80 ERRTXT | |
16 | PARAMETER (NAME = 'DSPPS1') | |
17 | ||
18 | ************************************************************************ | |
19 | * NORBAS, VERSION: 15.03.1993 | |
20 | ************************************************************************ | |
21 | * | |
22 | * DSPPS1 COMPUTES FUNCTION VALUES, VALUES OF DERIVATIVES, AND THE | |
23 | * VALUE OF INTEGRAL, RESPECTIVELY, OF A POLYNOMIAL SPLINE S(X) IN | |
24 | * B-SPLINE REPRESENTATION | |
25 | * | |
26 | * S(X) = SUMME(I=1,...,M-K-1) C(I) * B(I,K)(X) . | |
27 | * | |
28 | * THE FUNCTIONS B(I,K)(X) ARE NORMALIZED B-SPLINES OF DEGREE K | |
29 | * (0<= K <= 25) WITH INDEX I (1 <= I <= M-K-1) OVER A SET OF | |
30 | * SPLINE-KNOTS | |
31 | * T(1),T(2), ... ,T(M) ( M >= 2*K+2 ) | |
32 | * (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER | |
33 | * THAN K+1). | |
34 | * | |
35 | * THE FUNCTION VALUE OF THE NORMALIZED B-SPLINE B(I,K)(X) IS | |
36 | * IDENTICALLY ZERO OUTSIDE THE INTERVAL T(I) <= X < T(I+K+1). | |
37 | * | |
38 | * THE NORMALIZATION OF N(X) = B(I,K)(X) IS SUCH THAT THE INTGRAL OF | |
39 | * N(X) OVER THE WHOLE X-RANGE EQUALS | |
40 | * ( T(I+K+1) - T(I) ) / (K+1) . | |
41 | * | |
42 | * C(1),...,C(M-K-NDER-1) MUST CONTAIN THE COEFFICIENTS OF THE | |
43 | * POLYNOMIAL SPLINE S(X) OR ITS DERIVATIVE, REPRESENTED BY | |
44 | * NORMALIZED B-SPLINES OF DEGREE K-NDER . | |
45 | * | |
46 | * FOR TRANSFORMATION THE COEFFICIENTS OF THE POLYNOMIAL SPLINE S(X) | |
47 | * TO THE CORRESPONDING COEFFICIENTS OF THE NDER-TH DERIVATIVE OF | |
48 | * S(X) THE ROUTINE DSPCD1 MAY BE USED. | |
49 | * | |
50 | * ESPECIALLY FOR COMPUTING THE COEFFICIENTS C(1),...,C(M-K-1) OF | |
51 | * THE POLYNOMIAL VARIATION DIMINISHING SPLINE APPROXIMATION OF A USER | |
52 | * SUPPLIED FUNCTION F(X) THE ROUTINE DSPVD1 MAY BE USED. | |
53 | * | |
54 | * PARAMETERS: | |
55 | * | |
56 | * K (INTEGER) DEGREE (= ORDER - 1) OF B-SPLINES. | |
57 | * M (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES. | |
58 | * NDER (INTEGER) ON ENTRY, NDER MUST CONTAIN AN INTEGER VALUE .GE. -1 | |
59 | * = -1: DSPPS1 COMPUTES THE INTEGRAL OF S(TAU) OVER THE | |
60 | * RANGE TAU <= X. | |
61 | * = 0: DSPPS1 COMPUTES THE FUNCTION VALUE OF THE POLYNOMIAL | |
62 | * SPLINE S(X) AT X. | |
63 | * >= 1: DSPPS1 COMPUTES THE VALUE OF THE NDER-TH DERIVATIVE OF | |
64 | * S(X) AT X (IF NDER > K ZERO RETURNS). | |
65 | * X (DOUBLE PRECISION) ON ENTRY, X MUST CONTAIN THE VALUE OF THE | |
66 | * INDEPENDENT VARIABLE X OF S(X). | |
67 | * T (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M CONTAINING THE | |
68 | * KNOTS, ON ENTRY. | |
69 | * C (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M-K-NDER-1, | |
70 | * ON ENTRY C(J) CONTAINS THE J-TH COEFFICIENT OF THE B-SPLINE | |
71 | * REPRESENTATION OF S(X) OR OF ITS DERIVATIVE. | |
72 | * NERR (INTEGER) ERROR INDICATOR. ON EXIT: | |
73 | * = 0: NO ERROR DETECTED | |
74 | * = 1: AT LEAST ONE OF THE CONSTANTS K , M , NDER IS ILLEGAL | |
75 | * | |
76 | * USAGE: | |
77 | * | |
78 | * THE FUNCTION-CALL | |
79 | * Y = DSPPS1(K,M,NDER,X,T,C,NERR) | |
80 | * RETURNS | |
81 | * - THE VALUE OF THE INTEGRAL (NDER = -1) OR | |
82 | * - THE FUNCTION VALUE (NDER = 0 ) OR | |
83 | * - THE VALUE OF THE NDER-TH DERIVATIVE (NDER > 0 ) | |
84 | * OF THE POLYNOMIAL SPLINE (IN B-SPLINE REPRESENTATION) | |
85 | * S(X) = SUMME(I=1,...,M-K-1) C(I)*B(I,K)(X) | |
86 | * OF DEGREE K WITH THE SET OF KNOTS T(1),...,T(M) AT X. | |
87 | * | |
88 | * ERROR MESSAGES: | |
89 | * | |
90 | * IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT- | |
91 | * PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED: | |
92 | * K < 0 OR K > 25 OR | |
93 | * M < 2*K+2 OR | |
94 | * NDER < -1 . | |
95 | * | |
96 | ************************************************************************ | |
97 | ||
98 | PARAMETER (Z0 = 0 , Z1 = 1) | |
99 | ||
100 | NERR=1 | |
101 | IF(K .LT. 0 .OR. K .GT. 25) THEN | |
102 | WRITE(ERRTXT,101) 'K',K | |
103 | CALL MTLPRT(NAME,'E210.1',ERRTXT) | |
104 | ELSEIF(M .LT. 2*K+2) THEN | |
105 | WRITE(ERRTXT,101) 'M',M | |
106 | CALL MTLPRT(NAME,'E210.2',ERRTXT) | |
107 | ELSEIF(NDER .LT. -1) THEN | |
108 | WRITE(ERRTXT,101) 'NDER',NDER | |
109 | CALL MTLPRT(NAME,'E210.5',ERRTXT) | |
110 | ELSE | |
111 | ||
112 | NERR=0 | |
113 | DSPPS1=Z0 | |
114 | IF( X .LT. T(K+1) | |
115 | + .OR. X .GT. T(M-K) .AND. NDER .GE. 0 | |
116 | + .OR. K .LT. NDER ) RETURN | |
117 | ||
118 | IF(NDER .EQ. -1) THEN | |
119 | IF(X .GE. T(M-K)) THEN | |
120 | R=0 | |
121 | DO 10 I=1,M-K-1 | |
122 | 10 R=R+C(I)*(T(I+K+1)-T(I)) | |
123 | R=R/(K+1) | |
124 | ELSE | |
125 | KK=LKKSPL(X,T(K+1),M-2*K-1)+K | |
126 | R=0 | |
127 | DO 20 I=1,KK-K-2 | |
128 | 20 R=R+C(I)*(T(I+K+1)-T(I)) | |
129 | R=R/(K+1) | |
130 | IF(K .EQ. 0) THEN | |
131 | K1=MAX(1,KK-1) | |
132 | R=R+(X-T(K1))*C(K1) | |
133 | ELSE | |
134 | DO 50 I=MAX(1,KK-K-1),KK-1 | |
135 | CALL DVSET(K+1,Z0,B(1),B(2)) | |
136 | B(KK-I)=1/(T(KK)-T(KK-1)) | |
137 | ||
138 | DO 30 L=1,K | |
139 | DO 30 J=MAX(1,KK-I-L),MIN(K+1-L,KK-I) | |
140 | DIF=T(I+J+L)-T(I+J-1) | |
141 | B0=Z0 | |
142 | IF(DIF .NE. 0) B0=((X-T(I+J-1))*B(J)+(T(I+J+L)-X)*B(J+1))/DIF | |
143 | 30 B(J)=B0 | |
144 | S=Z0 | |
145 | DO 40 L=1,KK-I | |
146 | 40 S=S+(X-T(I+L-1))*B(L) | |
147 | S=S*(T(I+K+1)-T(I))/(K+1) | |
148 | 50 R=R+C(I)*S | |
149 | ENDIF | |
150 | ENDIF | |
151 | ELSE | |
152 | ||
153 | KK=LKKSPL(X,T(K+1),M-2*K-1)+K | |
154 | E1=X-T(KK-1) | |
155 | B(1)=Z1 | |
156 | DO 60 J=2,K-NDER+1 | |
157 | 60 B(J)=E1*B(J-1)/(T(KK-2+J)-T(KK-1)) | |
158 | IF(KK .NE. 0 .OR. NDER .NE. 0) THEN | |
159 | E2=T(KK)-X | |
160 | DO 70 J=1,K-NDER | |
161 | E3=X-T(KK-1-J) | |
162 | B(1)=E2*B(1)/(T(KK)-T(KK-J)) | |
163 | DO 70 L=2,K-NDER+1-J | |
164 | 70 B(L)=E3*B(L-1)/(T(KK-2+L)-T(KK-1-J))+ | |
165 | + (T(KK-1+L)-X)*B(L)/(T(KK-1+L)-T(KK-J)) | |
166 | ENDIF | |
167 | R=DVMPY(K-NDER+1,C(KK-1-K),C(KK-K),B(1),B(2)) | |
168 | ENDIF | |
169 | DSPPS1=R | |
170 | ENDIF | |
171 | RETURN | |
172 | ||
173 | 101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE') | |
174 | END | |
175 | ||
176 | ||
177 |