5 * Revision 1.1.1.1 1996/04/01 15:02:25 mclareni
10 SUBROUTINE DSPAP2(KX,KY,MX,MY,NX,NY,XI,YI,ZI,NDIMZ,
11 + KNOT,TX,TY,C,NDIMC,W,NW,NERR)
13 #include "gen/imp64.inc"
14 DIMENSION XI(*),YI(*),ZI(NDIMZ,*),TX(*),TY(*),W(*),C(NDIMC,*)
17 PARAMETER (NAME = 'DSPAP2')
19 ************************************************************************
20 * NORBAS, VERSION: 15.03.1993
21 ************************************************************************
23 * DSPAP2 COMPUTES THE COEFFICIENTS
24 * C(I,J) (I=1,...,NCX , J=1,...,NCY)
25 * OF A TWO-DIMENSIONAL POLYNOMIAL APPROXIMATION SPLINE Z = S(X,Y) IN
26 * REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL B-SPLINES B(I,J)(X,Y)
28 * S(X,Y) = SUMME(I=1,...,NCX)
29 * SUMME(J=1,...,NCY) C(I,J) * B(I,J)(X,Y)
31 * TO A USER SUPPLIED DATA SET
33 * (XI(I),YI(J),ZI(I,J)) (I=1,...,NX , J=1,...,NY)
35 * OF A FUNCTION Z = F(X,Y) , I.E.
37 * S(XI(I),YI(J)) = Z(I,J) (I=1,...,NX , J=1,...,NY) .
39 * THE TWO-DIMENSIONAL B-SPLINES B(I,J)(X,Y) ARE THE PRODUCT OF TWO
40 * ONE-DIMENSIONAL B-SPLINES BX , BY
41 * B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y)
42 * OF DEGREE KX AND KY ( 0 <= KX , KY <= 25 ) WITH INDICES I , J
43 * ( 1 <= I <= NX , 1 <= J <= NY ) OVER TWO SETS OF SPLINE-KNOTS
44 * TX(1),TX(2),...,TX(MX) ( MX = NCX+KX+1 )
45 * TY(1),TY(2),...,TY(MY) ( MY = NCY+KY+1 ) ,
47 * FOR FURTHER DETAILS TO THE ONE- AND TWO-DIMENSIONAL NORMALIZED
48 * B-SPLINES SEE THE COMMENTS TO DSPNB1 AND DSPNB2.
52 * NX (INTEGER) NUMBER OF APPROXIMATION POINTS IN X-DIRECTION :
53 * XI(I) , I=1,...,NX .
54 * NY (INTEGER) NUMBER OF APPROXIMATION POINTS IN Y-DIRECTION :
55 * YI(J) , J=1,...,NY .
56 * MX (INTEGER) NUMBER OF KNOTS IN X-DIRECTION .
57 * MY (INTEGER) NUMBER OF KNOTS IN Y-DIRECTION .
58 * KX (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION
59 * OVER THE SET OF KNOTS TX.
60 * KY (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION
61 * OVER THE SET OF KNOTS TY.
62 * NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY C IN THE
63 * CALLING PROGRAM, WITH NDIMC >= (MX-KX-1) .
64 * NDIMZ (INTEGER) DECLARED FIRST DIMENSION OF ARRAY ZI IN THE
65 * CALLING PROGRAM, WITH NDIMZ >= NX.
66 * XI (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NX .
67 * XI MUST CONTAIN THE APPROXIMATION POINTS IN X-DIRECTION IN
68 * ASCENDING ORDER, ON ENTRY.
69 * YI (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NY .
70 * YI MUST CONTAIN THE APPROXIMATION POINTS IN Y-DIRECTION IN
71 * ASCENDING ORDER, ON ENTRY.
72 * ZI (DOUBLE PRECISION) ARRAY OF ORDER (NDIMZ , >= NY) .
73 * ON ENTRY ZI MUST CONTAIN THE GIVEN FUNCTION VALUES
74 * Z(I,J) AT THE APPROXIMATION POINTS (X(I),Y(J))
75 * ( I=1,...,NX , J=1,...,NY ).
76 * KNOT (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS.
78 * = 1 : KNOTS ARE COMPUTED BY DSPAP2 IN THE FOLLOWING WAY:
79 * TX(J) = XI(1) , J = 1,...,KX+1
80 * TX(J) = XI(1)+(J-KX-1)*(XI(NX)-XI(1))/(NCX-KX) ,
82 * TX(NCX+J) = XI(NX) , J = 1,...,KX+1
83 * = 2 : KNOTS ARE COMPUTED BY DSPAP2 IN THE FOLLOWING WAY:
84 * TX(J) = XI(1) , J = 1,...,KX+1
85 * TX(J) = (XI(NX*(J-KX-2)/NCX+1)+XI(NX*J/NCX))/2 ,
87 * TX(NCX+J) = XI(NX) , J = 1,...,KX+1
88 * OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
89 * TX(1) <= ... <= TX(KX+1) <= XI(1)
90 * XI(1) < TX(KX+2) < ... < TX(NX) < XI(NX)
91 * XI(NX) <= TX(NX+1) <= ... <= TX(NX+KX+1)
92 * IN ALL CASES THE SAME CHOICE IS USED FOR KNOTS TY IN
94 * TX (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX.
95 * IF THE INPUT VALUE OF THE PARAMETER KNOT IS 1 OR 2 THE
96 * KNOTS ARE COMPUTED BY DSPAP2 AND THEY ARE GIVEN IN THE
98 * IN THE OTHER CASES THE ARRAY TX MUST CONTAIN THE USER
99 * SUPPLIED KNOTS IN X-DIRECTION, ON ENTRY.
100 * TY (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY.
101 * IF THE INPUT VALUE OF THE PARAMETER KNOT IS 1 OR 2 THE
102 * KNOTS ARE COMPUTED BY DSPAP2 AND THEY ARE GIVEN IN THE
104 * IN THE OTHER CASES THE ARRAY TY MUST CONTAIN THE USER
105 * SUPPLIED KNOTS IN Y-DIRECTION, ON ENTRY.
106 * W (DOUBLE PRECISION) WORKING ARRAY OF AT LEAST ORDER NW .
107 * NW (INTEGER) ORDER OF WORKING ARRAY W .
108 * NW >= N*(NC+6)+NC*(NC+1),
109 * WITH N=NX*NY , NC=NCX*NCY=(MX-KX-1)*(MY-KY-1)
110 * FOR GOOD PERFORMANCE, NW SHOULD GENERALLY BE LARGER.
111 * C (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= MY-KY-1).
112 * ON EXIT C(I,J) CONTAINS THE (I,J)-TH COEFFICIENT OF THE
113 * TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF S(X,Y) .
114 * NERR (INTEGER) ERROR INDICATOR. ON EXIT:
115 * = 0: NO ERROR DETECTED
116 * = 1: AT LEAST ONE OF THE CONSTANTS KX , KY , NX , NY , MX ,MY
121 * IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
122 * PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
123 * KX < 0 OR KX > 25 OR
124 * KY < 0 OR KY > 25 OR
125 * MX < 2*KX+2 OR MY < 2*KY+2 OR
126 * NX < MX-KX-1 OR NY < MY-KY-1 .
128 ************************************************************************
131 IF(KX .LT. 0 .OR. KX .GT. 25) THEN
132 WRITE(ERRTXT,101) 'KX',KX
133 CALL MTLPRT(NAME,'E210.1',ERRTXT)
134 ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN
135 WRITE(ERRTXT,101) 'KY',KY
136 CALL MTLPRT(NAME,'E210.1',ERRTXT)
137 ELSEIF(MX .LT. 2*KX+2) THEN
138 WRITE(ERRTXT,101) 'MX',MX
139 CALL MTLPRT(NAME,'E210.2',ERRTXT)
140 ELSEIF(MY .LT. 2*KY+2) THEN
141 WRITE(ERRTXT,101) 'MY',MY
142 CALL MTLPRT(NAME,'E210.2',ERRTXT)
143 ELSEIF(NX .LT. MX-KX-1) THEN
144 WRITE(ERRTXT,101) 'NX',NX
145 CALL MTLPRT(NAME,'E210.4',ERRTXT)
146 ELSEIF(NY .LT. MY-KY-1) THEN
147 WRITE(ERRTXT,101) 'NY',NY
148 CALL MTLPRT(NAME,'E210.4',ERRTXT)
162 CALL SPLAS2(N,NC,NCX,NCY,NX,NY,MX,MY,KX,KY,NDIMC,NDIMZ,XI,YI,ZI,
163 + KNOT,TX,TY,W(M1),W(M2),W(M3),W(M4),W(M5),LW,C,NERR)
169 101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE')