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