1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1890.c,v 1.3 2005-02-28 09:04:49 afr Exp $
45 *
46 */
47
48
49 #define S1890
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1890(double oknots[],int oik,int oin,double * par[],int * der[],int * jstat)55 s1890 (double oknots[], int oik, int oin, double *par[], int *der[],
56 int *jstat)
57 #else
58 void
59 s1890 (oknots, oik, oin, par, der, jstat)
60 double oknots[];
61 int oik;
62 int oin;
63 double *par[];
64 int *der[];
65 int *jstat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 *********************************************************************
71 *
72 * PURPOSE : To produce a set of parameter values and their corresponding
73 * derivative indicatores.
74 * The parameters will later be used in an interpolation problem.
75 *
76 * INPUT : oknots - The original knot vector.
77 * oik - The order of the basis.
78 * oin - The number of degrees of freedom in the basis given
79 * by the knot vector.
80 *
81 * OUTPUT : par - The computed parameter values.
82 * der - The derivate indicators. (all 0.0)
83 * jstat - Status variable:
84 * > 0 : warning
85 * = 0 : ok
86 * < 0 : error
87 *
88 * METHOD : First one parameter value is inserted at the start parameter
89 * value, then ik+1 knots are added successively and divided by
90 * ik+1 to produce internal parameter values. Then the end
91 * parameter value is inserted. This procedure can result in
92 * parameter values outside the legal range, these are then
93 * corrected. All parameter values are assigned derivative
94 * indicator 0.
95 *
96 * REFERENCES : Fortran version:
97 * Tor Dokken, SI, 1981-10
98 *
99 * CALLS : s6err.
100 *
101 * WRITTEN BY : Christophe R. Birkeland & Trond Vidar Stensby, SI, 1991-06
102 *
103 *********************************************************************
104 */
105 {
106 int kpos = 0;
107 int count1, count2; /* Loop control variables */
108 int start, stop;
109 int numb; /* Number of wrong parameters */
110
111 double sum; /* Sum of knot values */
112 double pvl; /* Single parameter value */
113 double delta; /* Used for correcting wrong
114 * parameter values */
115
116 *jstat = 0;
117
118
119 /* Test if legal input. */
120
121 if (oik <= 1 || oin < oik)
122 goto err112;
123
124
125 /* Test if input knot vector degenerate. */
126
127 if (oknots[oik - 1] >= oknots[oin])
128 goto err112;
129
130
131 /* Allocate arrays par and der. */
132
133 *par = newarray (oin, DOUBLE);
134 if (*par == SISL_NULL)
135 goto err101;
136 *der = new0array (oin, INT);
137 if (*der == SISL_NULL)
138 goto err101;
139
140
141 /* P R O D U C E P A R A M E T E R V A L U E S.
142 * First we produce parameter values by a simple algorithm.
143 * The parameter values calculated in a wrong way are then corrected. */
144
145 (*par)[0] = oknots[oik - 1];
146 (*par)[oin - 1] = oknots[oin];
147
148 for (count1 = 2; count1 < oin; count1++)
149 {
150 stop = count1 + oik;
151 sum = (double) 0.0;
152 for (count2 = count1; count2 <= stop; count2++)
153 sum += oknots[count2 - 1];
154 (*par)[count1 - 1] = sum / (oik + 1);
155 }
156
157 /* Find second distinct knot value. */
158
159 pvl = oknots[oik - 1];
160 for (count1 = oik; oknots[count1] <= pvl; count1++) ;
161
162
163 /* Find number of parameter values with wrong value at start of curve. */
164
165 pvl = (oknots[oik - 1] + oknots[count1]) / (double)2.0;
166 for (numb = 0, start = 1; (*par)[start] <= pvl; start++, numb++) ;
167
168 if (numb > 0)
169 {
170 delta = (pvl - (*par)[0]) / (numb + 1);
171
172 /* Fill inn missing parameter values. */
173
174 pvl = (*par)[0] + delta;
175
176 for (count1 = 1; count1 <= numb; count1++)
177 {
178 (*par)[count1] = pvl;
179 pvl += delta;
180 }
181 }
182
183 /* Find last but one distinct knot value. */
184
185 pvl = oknots[oin];
186 for (count1 = oin - 1; oknots[count1] >= pvl; count1--) ;
187
188
189 /* Find end parameters in wrong interval. */
190
191 pvl = (oknots[count1] + oknots[oin + 1]) / (double) 2.0;
192 for (numb = 0, stop = oin - 2; (*par)[stop] >= pvl; stop--, numb++) ;
193
194 if (numb > 0)
195 {
196 delta = ((*par)[oin - 1] - pvl) / (numb + 1);
197 pvl = (*par)[oin - 1] - delta;
198 for (count1 = 1; count1 <= numb; count1++)
199 {
200 (*par)[oin - 1 - count1] = pvl;
201 pvl -= delta;
202 }
203 }
204
205 /* Make derivative indicators */
206
207 /* We used new0array which initializes all elements with zeroes
208 * and then this code is redundant.
209 *
210 * for (count1 = 0; count1 < oin; count1++)
211 * (*der)[count1] = 0;
212 */
213 /* Knots produced */
214
215 goto out;
216
217
218 /* Not enough memory. */
219
220 err101:
221 *jstat = -101;
222 s6err ("s1890", *jstat, kpos);
223 goto out;
224
225 /* Error in description of B-spline. */
226
227 err112:
228 *jstat = -112;
229 s6err ("s1890", *jstat, kpos);
230 goto out;
231
232 out:
233 return;
234 }
235