]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:25 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE DSPAP2(KX,KY,MX,MY,NX,NY,XI,YI,ZI,NDIMZ, | |
11 | + KNOT,TX,TY,C,NDIMC,W,NW,NERR) | |
12 | ||
13 | #include "gen/imp64.inc" | |
14 | DIMENSION XI(*),YI(*),ZI(NDIMZ,*),TX(*),TY(*),W(*),C(NDIMC,*) | |
15 | CHARACTER NAME*(*) | |
16 | CHARACTER*80 ERRTXT | |
17 | PARAMETER (NAME = 'DSPAP2') | |
18 | ||
19 | ************************************************************************ | |
20 | * NORBAS, VERSION: 15.03.1993 | |
21 | ************************************************************************ | |
22 | * | |
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) | |
27 | * | |
28 | * S(X,Y) = SUMME(I=1,...,NCX) | |
29 | * SUMME(J=1,...,NCY) C(I,J) * B(I,J)(X,Y) | |
30 | * | |
31 | * TO A USER SUPPLIED DATA SET | |
32 | * | |
33 | * (XI(I),YI(J),ZI(I,J)) (I=1,...,NX , J=1,...,NY) | |
34 | * | |
35 | * OF A FUNCTION Z = F(X,Y) , I.E. | |
36 | * | |
37 | * S(XI(I),YI(J)) = Z(I,J) (I=1,...,NX , J=1,...,NY) . | |
38 | * | |
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 ) , | |
46 | * RESPECTIVELY. | |
47 | * FOR FURTHER DETAILS TO THE ONE- AND TWO-DIMENSIONAL NORMALIZED | |
48 | * B-SPLINES SEE THE COMMENTS TO DSPNB1 AND DSPNB2. | |
49 | * | |
50 | * PARAMETERS: | |
51 | * | |
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. | |
77 | * ON ENTRY: | |
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) , | |
81 | * J = KX+2,...,NCX | |
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 , | |
86 | * J = KX+2,...,NCX | |
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 | |
93 | * Y-DIRECTION . | |
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 | |
97 | * ARRAY TX, ON EXIT. | |
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 | |
103 | * ARRAY TY, ON EXIT. | |
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 | |
117 | * IS ILLEGAL | |
118 | * | |
119 | * ERROR MESSAGES: | |
120 | * | |
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 . | |
127 | * | |
128 | ************************************************************************ | |
129 | ||
130 | NERR=1 | |
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) | |
149 | ELSE | |
150 | ||
151 | NCX=MX-KX-1 | |
152 | NCY=MY-KY-1 | |
153 | NC=NCX*NCY | |
154 | N=NX*NY | |
155 | M1=1 | |
156 | M2=M1+N*NC | |
157 | M3=M2+NC | |
158 | M4=M3+NC*NC | |
159 | M5=M4+N | |
160 | LW=NW-M5+1 | |
161 | ||
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) | |
164 | ||
165 | ENDIF | |
166 | ||
167 | RETURN | |
168 | ||
169 | 101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE') | |
170 | END | |
171 | ||
172 | ||
173 |