This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / dspvd2.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       SUBROUTINE DSPVD2(F,KX,KY,MX,MY,TX,TY,C,NDIMC,NERR)
11
12 #include "gen/imp64.inc"
13       DIMENSION TX(*),TY(*),C(NDIMC,*)
14       CHARACTER NAME*(*)
15       CHARACTER*80 ERRTXT
16       PARAMETER (NAME = 'DSPVD2')
17
18 ************************************************************************
19 *   NORBAS, VERSION: 15.03.1993
20 ************************************************************************
21 *
22 *   DSPVD2 COMPUTES THE COEFFICIENTS
23 *          C(I,J)   (I=1,...,MX-KX-1 , J=1,...,MY-KY-1)
24 *   OF A TWO-DIMENSIONAL POLYNOMIAL VARIATION DIMINISHING SPLINE
25 *   APPROXIMATION S(X,Y) IN REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL
26 *   B-SPLINES  B(I,J)(X,Y)
27 *
28 *     S(X,Y) = SUMME(I=1,...,MX-KX-1)
29 *              SUMME(J=1,...,MY-KY-1) C(I,J) * B(I,J)(X,Y) .
30 *
31 *   THE TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)  ARE THE PRODUCT OF TWO
32 *   ONE-DIMENSIONAL B-SPLINES  BX , BY
33 *          B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y)
34 *   OF DEGREE  KX AND KY ( 0 <= KX , KY <= 25 )   WITH INDICES  I , J
35 *   ( 1 <= I <= MX-KX-1 , 1 <= J <= MY-KY-1 ) OVER TWO SETS OF SPLINE-
36 *   KNOTS
37 *       TX(1),TX(2),...,TX(MX)   ( MX >= 2*KX+2 )
38 *       TY(1),TY(2),...,TY(MY)   ( MY >= 2*KY+2 ) ,
39 *   RESPECTIVELY.
40 *   FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE
41 *   THE COMMENTS TO  DSPNB1.
42 *
43 *   PARAMETERS:
44 *
45 *   F     (DOUBLE PRECISION) USER SUPPLIED FUNCTION  Z=F(X,Y)  FOR WHICH
46 *         THE CORRESPONDING SPLINE APPROXIMATION HAS TO BE COMPUTED.
47 *   KX    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION
48 *         OVER THE SET OF KNOTS  TX.
49 *   KY    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION
50 *         OVER THE SET OF KNOTS  TY.
51 *   MX    (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN X-DIRECTION.
52 *   MY    (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN Y-DIRECTION.
53 *   NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  C  IN THE
54 *         CALLING PROGRAM, WITH  NDIMC >= MX-KX-1 .
55 *   TX    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX CONTAINING THE
56 *         KNOTS IN X-DIRECTION, ON ENTRY.
57 *   TY    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY CONTAINING THE
58 *         KNOTS IN Y-DIRECTION, ON ENTRY.
59 *   C     (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= MY-KY-1).
60 *         ON EXIT  C(I,J)  CONTAINS THE (I,J)-TH COEFFICIENT OF THE
61 *         TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF  S(X,Y) .
62 *   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
63 *         = 0: NO ERROR DETECTED
64 *         = 1: AT LEAST ONE OF THE CONSTANTS KX, KY, MX, MY IS ILLEGAL
65 *
66 *   ERROR MESSAGES:
67 *
68 *   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
69 *   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
70 *     KX < 1       OR    KX > 25      OR    KY < 1   OR   KY > 25   OR
71 *     MX < 2*KX+2  OR    MY < 2*KY+2  .
72 *
73 ************************************************************************
74
75       PARAMETER (Z1 = 1)
76
77       NERR=1
78       IF(KX .LT. 1 .OR. KX .GT. 25) THEN
79        WRITE(ERRTXT,101) 'KX',KX
80        CALL MTLPRT(NAME,'E210.1',ERRTXT)
81       ELSEIF(KY .LT. 1 .OR. KY .GT. 25) THEN
82        WRITE(ERRTXT,101) 'KY',KY
83        CALL MTLPRT(NAME,'E210.1',ERRTXT)
84       ELSEIF(MX .LT. 2*KX+2) THEN
85        WRITE(ERRTXT,101) 'MX',MX
86        CALL MTLPRT(NAME,'E210.2',ERRTXT)
87       ELSEIF(MY .LT. 2*KY+2) THEN
88        WRITE(ERRTXT,101) 'MY',MY
89        CALL MTLPRT(NAME,'E210.2',ERRTXT)
90       ELSE
91
92       NERR=0
93       RX=Z1/KX
94       RY=Z1/KY
95       DO 10 I = 1,MX-KX-1
96       XI=RX*DVSUM(KX,TX(I+1),TX(I+2))
97       DO 10 J = 1,MY-KY-1
98    10 C(I,J)=F(XI,RY*DVSUM(KY,TY(J+1),TY(J+2)))
99       ENDIF
100
101       RETURN
102
103   101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
104       END
105
106
107