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