]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MINICERN/mathlib/gen/e/dspps2.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / dspps2.F
CommitLineData
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