]>
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 DSPPS2(KX,KY,MX,MY,NDERX,NDERY,X,Y,TX,TY,C,NDIMC,W,NERR) | |
11 | ||
12 | #include "gen/imp64.inc" | |
13 | DIMENSION TX(*),TY(*),C(NDIMC,*),W(*),BX(27),BY(27) | |
14 | CHARACTER NAME*(*) | |
15 | CHARACTER*80 ERRTXT | |
16 | PARAMETER (NAME = 'DSPPS2') | |
17 | ||
18 | ************************************************************************ | |
19 | * NORBAS, VERSION: 15.03.1993 | |
20 | ************************************************************************ | |
21 | * | |
22 | * DSPPS2 COMPUTES FUNCTION VALUES, VALUES OF DERIVATIVES, AND THE | |
23 | * VALUE OF INTEGRAL, RESPECTIVELY, OF A TWO-DIMENSIONAL POLYNOMIAL | |
24 | * SPLINE S(X,Y) IN REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL | |
25 | * B-SPLINES B(I,J)(X,Y) | |
26 | * | |
27 | * S(X,Y) = SUMME(I=1,...,MX-KX-1) | |
28 | * SUMME(J=1,...,MY-KY-1) C(I,J) * B(I,J)(X,Y) . | |
29 | * | |
30 | * THE TWO-DIMENSIONAL B-SPLINES B(I,J)(X,Y) ARE THE PRODUCT OF TWO | |
31 | * ONE-DIMENSIONAL B-SPLINES BX , BY | |
32 | * B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y) | |
33 | * OF DEGREE KX AND KY ( 0 <= KX , KY <= 25 ) WITH INDICES I , J | |
34 | * ( 1 <= I <= MX-KX-1 , 1 <= J <= MY-KY-1 ) OVER TWO SETS OF SPLINE- | |
35 | * KNOTS | |
36 | * TX(1),TX(2),...,TX(MX) ( MX >= 2*KX+2 ) | |
37 | * TY(1),TY(2),...,TY(MY) ( MY >= 2*KY+2 ) , | |
38 | * RESPECTIVELY. | |
39 | * FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE | |
40 | * THE COMMENTS TO DSPNB1. | |
41 | * | |
42 | * C(I,J) (I=1,...,MX-KX-NDERX-1 , J=1,...,MY-KY-NDERY-1) MUST CONTAIN | |
43 | * THE (I,J)-TH C COEFFICIENT OF THE POLYNOMIAL SPLINE S(X,Y) OR | |
44 | * OF ONE OF ITS PARTIAL DERIVATIVE, REPRESENTED BY NORMALIZED | |
45 | * TWO-DIMENSIONAL B-SPLINES OF DEGREE (KX-NDERX) AND (KY-NDERY), | |
46 | * RESPECTIVELY. | |
47 | * | |
48 | * FOR TRANSFORMATION THE COEFFICIENTS OF THE POLYNOMIAL SPLINE S(X,Y) | |
49 | * TO THE CORRESPONDING COEFFICIENTS OF THE NDERX-TH AND NDERY-TH | |
50 | * PARTIAL DERIVATIVE OF S(X,Y) THE ROUTINE DSPCD2 MAY BE USED. | |
51 | * | |
52 | * ESPECIALLY FOR COMPUTING THE COEFFICIENTS C(I,J) OF THE TWO- | |
53 | * DIMENSIONAL POLYNOMIAL VARIATION DIMINISHING SPLINE APPREOXIMATION | |
54 | * OF A USER SUPPLIED FUNCTION Z = F(X,Y) THE ROUTINE DSPVD2 MAY BE | |
55 | * USED. | |
56 | * | |
57 | * PARAMETERS: | |
58 | * | |
59 | * KX (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION | |
60 | * OVER THE SET OF KNOTS TX. | |
61 | * KY (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION | |
62 | * OVER THE SET OF KNOTS TY. | |
63 | * MX (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN X-DIRECTION. | |
64 | * MY (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN Y-DIRECTION. | |
65 | * NDERX (INTEGER) ON ENTRY, NDERX MUST CONTAIN AN INTEGER VALUE >= -1. | |
66 | * = -1: DSPPS2 COMPUTES THE INTEGRAL OF S(TAU,Y) OVER THE | |
67 | * RANGE TAU <= X. | |
68 | * = 0: DSPPS2 COMPUTES THE FUNCTION VALUE S(X,Y) FOR | |
69 | * FOR THE SPECIFIED VALUES OF X,Y. | |
70 | * >= 1: DSPNB2 COMPUTES THE VALUE OF THE NDERX-TH PARTIAL | |
71 | * DERIVATIVE OF S(X,Y) WITH RESPECT TO X | |
72 | * FOR THE SPECIFIED VALUES OF X,Y. | |
73 | * (IF NDERX > KX ZERO RETURNS). | |
74 | * NDERY (INTEGER) ON ENTRY, NDERY MUST CONTAIN AN INTEGER VALUE >= -1. | |
75 | * THE MEANING OF NDERY IS THE SAME AS THAT OF THE PARAMETER | |
76 | * NDERX WITH RESPECT TO Y-DIRECTION INSTEAD OF X-DIRECTION. | |
77 | * NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY C IN THE | |
78 | * CALLING PROGRAM, WITH NDIMC >= MX-KX-NDERX-1 . | |
79 | * X (DOUBLE PRECISION) ON ENTRY, X MUST CONTAIN THE VALUE OF THE | |
80 | * INDEPENDENT VARIABLE X OF S(X,Y) | |
81 | * Y (DOUBLE PRECISION) ON ENTRY, Y MUST CONTAIN THE VALUE OF THE | |
82 | * INDEPENDENT VARIABLE Y OF S(X,Y) | |
83 | * TX (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX CONTAINING THE | |
84 | * KNOTS IN X-DIRECTION, ON ENTRY. | |
85 | * TY (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY CONTAINING THE | |
86 | * KNOTS IN Y-DIRECTION, ON ENTRY. | |
87 | * C (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= MY-KY-NDERY-1), | |
88 | * CONTAINING THE COEFFICIENTS OF THE TWO-DIMENSIONAL B-SPLINE | |
89 | * REPRESENTATION OF S(X,Y) , ON ENTRY. | |
90 | * W (DOUBLE PRECISION) WORKING ARRAY OF AT LEAST ORDER MY-KY-1. | |
91 | * NERR (INTEGER) ERROR INDICATOR. ON EXIT: | |
92 | * = 0: NO ERROR DETECTED | |
93 | * = 1: AT LEAST ONE OF THE CONSTANTS KX , KY , MX , MY , | |
94 | * NDERX , NDERY IS ILLEGAL | |
95 | * | |
96 | * USAGE: | |
97 | * | |
98 | * THE FUNCTION-CALL | |
99 | * Z = DSPPS2(KX,KY,MX,MY,NDERX,NDERY,X,Y,TX,TY,C,NDIMC,W,NERR) | |
100 | * RETURNS | |
101 | * - THE VALUE OF THE INTEGRAL (NDERX =-1 .AND. NDERY =-1) OR | |
102 | * - THE FUNCTION VALUE (NDERX = 0 .AND. NDERY = 0) OR | |
103 | * - THE VALUE OF THE NDERX-TH AND NDERY-TH PARTIAL DERIVATIVE | |
104 | * (NDERX >= 0 .AND. NDERY >= 0 .AND. NDERX + NDERY > 0) | |
105 | * OF THE POLYNOMIAL SPLINE S(X,Y) AT (X,Y). | |
106 | * | |
107 | * ERROR MESSAGES: | |
108 | * | |
109 | * IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT- | |
110 | * PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED: | |
111 | * KX < 0 OR KX > 25 OR KY < 0 OR KY > 25 OR | |
112 | * MX < 2*KX+2 OR MY < 2*KY+2 OR | |
113 | * NDERX < -1 OR NDERY < -1 . | |
114 | * | |
115 | ************************************************************************ | |
116 | ||
117 | PARAMETER (Z0 = 0 , Z1 = 1) | |
118 | ||
119 | NERR=1 | |
120 | IF(KX .LT. 0 .OR. KX .GT. 25) THEN | |
121 | WRITE(ERRTXT,101) 'KX',KX | |
122 | CALL MTLPRT(NAME,'E210.1',ERRTXT) | |
123 | ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN | |
124 | WRITE(ERRTXT,101) 'KY',KY | |
125 | CALL MTLPRT(NAME,'E210.1',ERRTXT) | |
126 | ELSEIF(MX .LT. 2*KX+2) THEN | |
127 | WRITE(ERRTXT,101) 'MX',MX | |
128 | CALL MTLPRT(NAME,'E210.2',ERRTXT) | |
129 | ELSEIF(MY .LT. 2*KY+2) THEN | |
130 | WRITE(ERRTXT,101) 'MY',MY | |
131 | CALL MTLPRT(NAME,'E210.2',ERRTXT) | |
132 | ELSEIF(NDERX .LT. -1) THEN | |
133 | WRITE(ERRTXT,101) 'NDERX',NDERX | |
134 | CALL MTLPRT(NAME,'E210.5',ERRTXT) | |
135 | ELSEIF(NDERY .LT. -1) THEN | |
136 | WRITE(ERRTXT,101) 'NDERY',NDERY | |
137 | CALL MTLPRT(NAME,'E210.5',ERRTXT) | |
138 | ELSEIF(NDERX .EQ. -1 .AND. NDERY .NE. -1 .OR. | |
139 | + NDERX .NE. -1 .AND. NDERY .EQ. -1 ) THEN | |
140 | WRITE(ERRTXT,102) 'NDERX',NDERX,'NDERY',NDERY | |
141 | CALL MTLPRT(NAME,'E210.6',ERRTXT) | |
142 | ELSE | |
143 | ||
144 | NERR=0 | |
145 | IF(NDERX .EQ. -1 .AND. NDERY .EQ. -1) THEN | |
146 | DO 60 JJ=1,MY-KY-1 | |
147 | IF(X .GE. TX(MX-KX)) THEN | |
148 | R=Z0 | |
149 | DO 10 I=1,MX-KX-1 | |
150 | 10 R=R+C(I,JJ)*(TX(I+KX+1)-TX(I)) | |
151 | R=R/(KX+1) | |
152 | ELSE | |
153 | KK=LKKSPL(X,TX(KX+1),MX-2*KX-1)+KX | |
154 | R=Z0 | |
155 | DO 20 I=1,KK-KX-2 | |
156 | 20 R=R+C(I,JJ)*(TX(I+KX+1)-TX(I)) | |
157 | R=R/(KX+1) | |
158 | IF(KX .EQ. 0) THEN | |
159 | K1=MAX(1,KK-1) | |
160 | R=R+(X-TX(K1))*C(K1,JJ) | |
161 | ELSE | |
162 | DO 50 I=MAX(1,KK-KX-1),KK-1 | |
163 | CALL DVSET(KX+1,Z0,BX(1),BX(2)) | |
164 | BX(KK-I)=1/(TX(KK)-TX(KK-1)) | |
165 | DO 30 L=1,KX | |
166 | DO 30 J=MAX(1,KK-I-L),MIN(KX+1-L,KK-I) | |
167 | DIF=TX(I+J+L)-TX(I+J-1) | |
168 | B0=Z0 | |
169 | IF(DIF .NE. 0) | |
170 | + B0=((X-TX(I+J-1))*BX(J)+(TX(I+J+L)-X)*BX(J+1))/DIF | |
171 | 30 BX(J)=B0 | |
172 | S=Z0 | |
173 | DO 40 L=1,KK-I | |
174 | 40 S=S+(X-TX(I+L-1))*BX(L) | |
175 | S=S*(TX(I+KX+1)-TX(I))/(KX+1) | |
176 | 50 R=R+C(I,JJ)*S | |
177 | ENDIF | |
178 | ENDIF | |
179 | 60 W(JJ)=R | |
180 | IF(Y .GE. TY(MY-KY)) THEN | |
181 | R=Z0 | |
182 | DO 70 I=1,MY-KY-1 | |
183 | 70 R=R+W(I)*(TY(I+KY+1)-TY(I)) | |
184 | R=R/(KY+1) | |
185 | ELSE | |
186 | KK=LKKSPL(Y,TY(KY+1),MY-2*KY-1)+KY | |
187 | R=Z0 | |
188 | DO 80 I=1,KK-KY-2 | |
189 | 80 R=R+W(I)*(TY(I+KY+1)-TY(I)) | |
190 | R=R/(KY+1) | |
191 | IF(KY .EQ. 0) THEN | |
192 | K1=MAX(1,KK-1) | |
193 | R=R+(Y-TY(K1))*W(K1) | |
194 | ELSE | |
195 | DO 110 I=MAX(1,KK-KY-1),KK-1 | |
196 | CALL DVSET(KY+1,Z0,BY(1),BY(2)) | |
197 | BY(KK-I)=1/(TY(KK)-TY(KK-1)) | |
198 | DO 90 L=1,KY | |
199 | DO 90 J=MAX(1,KK-I-L),MIN(KY+1-L,KK-I) | |
200 | DIF=TY(I+J+L)-TY(I+J-1) | |
201 | B0=Z0 | |
202 | IF(DIF .NE. 0) | |
203 | + B0=((Y-TY(I+J-1))*BY(J)+(TY(I+J+L)-Y)*BY(J+1))/DIF | |
204 | 90 BY(J)=B0 | |
205 | S=Z0 | |
206 | DO 100 L=1,KK-I | |
207 | 100 S=S+(Y-TY(I+L-1))*BY(L) | |
208 | S=S*(TY(I+KY+1)-TY(I))/(KY+1) | |
209 | 110 R=R+W(I)*S | |
210 | ENDIF | |
211 | ENDIF | |
212 | DSPPS2=R | |
213 | RETURN | |
214 | ENDIF | |
215 | ||
216 | DSPPS2=Z0 | |
217 | IF(X .LT. TX(KX+1) .OR. X .GT. TX(MX-KX) .OR. | |
218 | + Y .LT. TY(KY+1) .OR. Y .GT. TY(MY-KY) .OR. | |
219 | + KX .LT. NDERX .OR. KY .LT. NDERY ) RETURN | |
220 | ||
221 | KKX=LKKSPL(X,TX(KX+1),MX-2*KX-1)+KX | |
222 | KKY=LKKSPL(Y,TY(KY+1),MY-2*KY-1)+KY | |
223 | E1=X-TX(KKX-1) | |
224 | BX(1)=Z1 | |
225 | DO 120 J=2,KX-NDERX+1 | |
226 | 120 BX(J)=E1*BX(J-1)/(TX(KKX-2+J)-TX(KKX-1)) | |
227 | E2=TX(KKX)-X | |
228 | DO 130 J=1,KX-NDERX | |
229 | E3=X-TX(KKX-1-J) | |
230 | BX(1)=E2*BX(1)/(TX(KKX)-TX(KKX-J)) | |
231 | DO 130 L=2,KX-NDERX+1-J | |
232 | 130 BX(L)=E3*BX(L-1)/(TX(KKX-2+L)-TX(KKX-1-J))+ | |
233 | + (TX(KKX-1+L)-X)*BX(L)/(TX(KKX-1+L)-TX(KKX-J)) | |
234 | E1=Y-TY(KKY-1) | |
235 | BY(1)=Z1 | |
236 | DO 140 J=2,KY-NDERY+1 | |
237 | 140 BY(J)=E1*BY(J-1)/(TY(KKY-2+J)-TY(KKY-1)) | |
238 | E2=TY(KKY)-Y | |
239 | DO 150 J=1,KY-NDERY | |
240 | E3=Y-TY(KKY-1-J) | |
241 | BY(1)=E2*BY(1)/(TY(KKY)-TY(KKY-J)) | |
242 | DO 150 L=2,KY-NDERY+1-J | |
243 | 150 BY(L)=E3*BY(L-1)/(TY(KKY-2+L)-TY(KKY-1-J))+ | |
244 | + (TY(KKY-1+L)-Y)*BY(L)/(TY(KKY-1+L)-TY(KKY-J)) | |
245 | R=Z0 | |
246 | DO 160 I=1,KX-NDERX+1 | |
247 | DO 160 J=1,KY-NDERY+1 | |
248 | 160 R=R+C(KKX-2-KX+I,KKY-2-KY+J)*BX(I)*BY(J) | |
249 | DSPPS2=R | |
250 | ENDIF | |
251 | RETURN | |
252 | ||
253 | 101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE') | |
254 | 102 FORMAT(1X,A5,' =',I6,A7,' =',I6,' INCONSISTENT') | |
255 | END | |
256 | ||
257 | ||
258 |