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: s1616.c,v 1.2 2001-03-19 15:58:52 afr Exp $
45 *
46 */
47
48
49 #define S1616
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
s1616(double epoint[],int inbpnt,int idim,int eptyp[],double econic[],int * jstat)54 void s1616(double epoint[], int inbpnt, int idim, int eptyp[],
55 double econic[], int *jstat)
56 #else
57 void s1616(epoint, inbpnt, idim, eptyp, econic, jstat)
58 double epoint[];
59 int inbpnt;
60 int idim;
61 int eptyp[];
62 double econic[];
63 int *jstat;
64 #endif
65 /*
66 *************************************************************************
67 *
68 * PURPOSE: To find if the conic equation of the points in the array
69 * epoint. The first and last points must be of type 1 or 2. The
70 * intermediate points can be points or tangents, but no tangents
71 * can be doubly defined. The routine s1614 should be run before this
72 * routine to ensure that the description of the points/tangents in the
73 * eptyp array are ok.
74 *
75 * INPUT:
76 * Epoint - The points/tangents describing the conic.
77 * Inbpnt - No. of points/tangents in the epoint array.
78 * Idim - The dimension of the space in which the points lie.
79 * Eptyp - Type indicator for the points/tangents :
80 * 1 - Ordinary point.
81 * 2 - Knuckle point. (Is treated as an ordinary point.)
82 * 3 - Tangent to next point.
83 * 4 - Tangent to prior point.
84 *
85 * Output:
86 * Econic - The conic coefficients of the points in Epoint.
87 * Jstat - status messages:
88 * < 0 : Error.
89 * = 0 : Ok.
90 * > 0 : Warning.
91 *
92 * Method:
93 * If Inbpnt is 3, a circle is produced.
94 * If Inbpnt is 4, a conic with the coeficient of x*y equal to
95 * zero is produced.
96 * If Inbpnt is 5, a general conic is produced.
97 * The coeficients econic have the following meaning:
98 * econic(0)*x*x + 2*econic(1)*x*y + econic(2)*y*y +
99 * econic(3)*x + 2*econic(4)*y + econic(5) = 0.0.
100 *-
101 * Calls: s6lufacp, s6lusolp, s1618, s6err.
102 *
103 * Written by: A.M. Ytrehus, SI Oslo Oct.91.
104 * After FORTRAN, (P10634), written by: T. Dokken SI.
105 * Revised by: J. Kaasa, SI, Aug. 1992 (Transposed smatrix and fixed the
106 * handling of singular matrix equations.)
107 * Revised by: J. Kaasa, SI, Aug. 1992 (Changed status 9005 to 105).
108 *****************************************************************
109 */
110 {
111 int ki, kj, kk;
112 int kdim;
113 int kn;
114 double tmax = (double) HUGE;
115 int kperm = 0;
116 int ktyp;
117 int kdir;
118 double *smatrix = SISL_NULL;
119 double *save = SISL_NULL;
120 double *sright = SISL_NULL;
121 double tx, ty, tdx, tdy;
122 int *npiv = SISL_NULL;
123 double solu[6];
124 double tdiff;
125 double tdum;
126 int kpos = 0;
127 int kstat = 0;
128
129 *jstat = 0;
130
131
132 /* Allocate local arrays. */
133
134 smatrix = newarray (inbpnt * inbpnt, DOUBLE);
135 if (smatrix == SISL_NULL)
136 goto err101;
137
138 sright = newarray (inbpnt, DOUBLE);
139 if (sright == SISL_NULL)
140 goto err101;
141
142 save = newarray (inbpnt * inbpnt, DOUBLE);
143 if (save == SISL_NULL)
144 goto err101;
145
146 npiv = newarray (inbpnt, INT);
147 if (npiv == SISL_NULL)
148 goto error;
149
150 kn = inbpnt + 1;
151 kdim = inbpnt * inbpnt;
152
153
154 /* Build interpolation matrix to find the coeficients of the conic. */
155
156 if (inbpnt == 3)
157 {
158
159 /* A circle is to be produced. We put econic(0) = econic(2) = 1.0,
160 and econic(1) = 0.0. (This is not done now but later.) */
161
162 for (ki = 0; ki < inbpnt; ki++)
163 {
164 kk = idim * ki;
165 ktyp = eptyp[ki];
166
167 if (ktyp < 3)
168 {
169 tx = epoint[kk];
170 ty = epoint[kk + 1];
171
172 smatrix[inbpnt*ki] = ((double) 2.0) * tx;
173 smatrix[inbpnt*ki + 1] = ((double) 2.0) * ty;
174 smatrix[inbpnt*ki + 2] = (double) 1.0;
175 sright[ki] = -(tx * tx + ty * ty);
176 }
177 else if (ktyp > 2)
178 {
179
180 /* Derivative condition. */
181
182 kdir = 1;
183 if (ktyp == 4)
184 kdir = -1;
185
186 tx = epoint[kk + idim * kdir];
187 ty = epoint[kk + idim * kdir + 1];
188 tdx = epoint[kk];
189 tdy = epoint[kk + 1];
190
191 smatrix[inbpnt*ki] = ((double) 2.0) * tdx;
192 smatrix[inbpnt*ki + 1] = ((double) 2.0) * tdy;
193 smatrix[inbpnt*ki + 2] = (double) 0.0;
194 sright[ki] = -tdx * tx * ((double) 2.0) - tdy * ty * ((double) 2.0);
195 }
196 }
197 }
198 else if (inbpnt == 4)
199 {
200 /* A conic with xy term equal to zero is to be produced.
201 We put econic(1) = 0.0. (This is not done now but later.) */
202
203 for (ki = 0; ki < inbpnt; ki++)
204 {
205 kk = idim * ki;
206 ktyp = eptyp[ki];
207 if (ktyp < 3)
208 {
209 tx = epoint[kk];
210 ty = epoint[kk + 1];
211
212 smatrix[inbpnt*ki] = tx * tx;
213 smatrix[inbpnt*ki + 1] = ty * ty;
214 smatrix[inbpnt*ki + 2] = ((double) 2.0) * tx;
215 smatrix[inbpnt*ki + 3] = ((double) 2.0) * ty;
216 sright[ki] = -(double) 1.0;
217 }
218 else if (ktyp > 2)
219 {
220
221 /* Derivative condition. */
222
223 kdir = 1;
224 if (ktyp == 4)
225 kdir = -1;
226
227 tx = epoint[kk + idim * kdir];
228 ty = epoint[kk + idim * kdir + 1];
229 tdx = epoint[kk];
230 tdy = epoint[kk + 1];
231
232 smatrix[inbpnt*ki] = ((double) 2.0) * tdx * tx;
233 smatrix[inbpnt*ki + 1] = ((double) 2.0) * tdy * ty;
234 smatrix[inbpnt*ki + 2] = ((double) 2.0) * tdx;
235 smatrix[inbpnt*ki + 3] = ((double) 2.0) * tdy;
236 sright[ki] = (double) 0.0;
237 }
238 }
239 }
240 else if (inbpnt == 5)
241 {
242
243 /* A general conic is to be produced. */
244
245 for (ki = 0; ki < inbpnt; ki++)
246 {
247 kk = idim * ki;
248 ktyp = eptyp[ki];
249 if (ktyp < 3)
250 {
251 tx = epoint[kk];
252 ty = epoint[kk + 1];
253
254 smatrix[inbpnt*ki] = tx * tx;
255 smatrix[inbpnt*ki + 1] = ((double) 2.0) * tx * ty;
256 smatrix[inbpnt*ki + 2] = ty * ty;
257 smatrix[inbpnt*ki + 3] = ((double) 2.0) * tx;
258 smatrix[inbpnt*ki + 4] = ((double) 2.0) * ty;
259 sright[ki] = -(double) 1.0;
260 }
261 else if (ktyp > 2)
262 {
263
264 /* Derivative condition. */
265
266 kdir = 1;
267 if (ktyp == 4)
268 kdir = -1;
269
270 tx = epoint[kk + idim * kdir];
271 ty = epoint[kk + idim * kdir + 1];
272 tdx = epoint[kk];
273 tdy = epoint[kk + 1];
274
275 smatrix[inbpnt*ki] = ((double) 2.0) * tdx * tx;
276 smatrix[inbpnt*ki + 1] = ((double) 2.0) * tdy * tx + ((double) 2.0) * tdx * ty;
277 smatrix[inbpnt*ki + 2] = ((double) 2.0) * tdy * ty;
278 smatrix[inbpnt*ki + 3] = ((double) 2.0) * tdx;
279 smatrix[inbpnt*ki + 4] = ((double) 2.0) * tdy;
280 sright[ki] = (double) 0.0;
281 }
282 }
283 }
284
285 /* Solve the equation system. If the system is not solvable, interchange
286 the r.h.side with one of the colomns on the l.h. side. */
287
288 for (ki = 0; ki < kn; ki++)
289 {
290 /* Remember interpolation matrix. */
291
292 for (kj = 0; kj < kdim; kj++)
293 save[kj] = smatrix[kj];
294
295
296 /* Find solution.
297 s6lusolp put the soltion into the r.h. side array given as input,
298 thus store the right hand side in econic. */
299
300 for (kj = 0; kj < inbpnt; kj++)
301 econic[kj] = sright[kj];
302
303 s6lufacp (smatrix, npiv, inbpnt, &kstat);
304
305 if (kstat >= 0 && kstat != 1)
306 s6lusolp (smatrix, econic, npiv, inbpnt, &kstat);
307 kstat = 0;
308
309 /* If we are here, we have been able to find a solution, econic.
310 Test this. First, restore smatrix and sright. */
311
312 for (kj = 0; kj < kdim; kj++)
313 smatrix[kj] = save[kj];
314
315 s1618 (smatrix, sright, econic, inbpnt, &tdiff);
316
317 if (tdiff < tmax)
318 {
319 /* If we are here, we have found the best solution until now. */
320
321 tmax = tdiff;
322 kperm = ki;
323
324 for (kj = 0; kj < inbpnt; kj++)
325 solu[kj] = econic[kj];
326
327 if (inbpnt == 3) break;
328 }
329
330 if (ki < (kn - 1))
331 {
332 /* Change the right hand side with one column of the matrix. */
333
334 for (kj = 0; kj < inbpnt; kj++)
335 {
336 kk = inbpnt * kj;
337 tdum = -sright[kj];
338 sright[kj] = -smatrix[kk + ki];
339 smatrix[kk + ki] = tdum;
340 }
341 }
342 }
343
344 /* If no solution found, make straight line. */
345
346 if (tmax > 0.0001)
347 {
348 *jstat = 105;
349 econic[0] = 0.;
350 econic[1] = 0.;
351 econic[2] = 0.;
352 econic[3] = (epoint[1] - epoint[2*inbpnt - 1])/2.;
353 econic[4] = (epoint[2*inbpnt - 2] - epoint[0])/2.;
354 econic[5] = epoint[0]*epoint[2*inbpnt - 1] -
355 epoint[1]*epoint[2*inbpnt - 2];
356 goto out;
357 }
358
359 /* Get the remembered best solution. */
360
361 for (kj = 0; kj < inbpnt; kj++)
362 econic[kj] = solu[kj];
363
364 econic[inbpnt] = (double) 1.0;
365
366
367 /* If the order of the colomns were changed, substitute back again. */
368
369 if (kperm != 0)
370 {
371 for (ki = 1; ki <= kperm; ki++)
372 {
373 kj = kperm - ki;
374 tdum = econic[kj];
375 econic[kj] = econic[inbpnt];
376 econic[inbpnt] = tdum;
377 }
378 }
379
380 /* Expand to full conic description when Inbpnt is less than 5. */
381
382 if (inbpnt == 3)
383 {
384 /* Circle. Only x, y and constant coeff. calculated. */
385
386 for (ki = 0; ki < inbpnt; ki++)
387 econic[ki + inbpnt] = econic[ki];
388
389 econic[0] = (double) 1.0;
390 econic[1] = (double) 0.0;
391 econic[2] = (double) 1.0;
392 }
393 else if (inbpnt == 4)
394 {
395 /* Conic with xy constant zero. */
396
397 for (ki = 5; ki > 1; ki--)
398 econic[ki] = econic[ki - 1];
399
400 econic[1] = (double) 0.0;
401 }
402
403 goto out;
404
405 /* Error in allocation. */
406
407 err101:
408 *jstat = -101;
409 s6err ("s1616", *jstat, kpos);
410 goto out;
411
412 /* Error in lower level routine. */
413
414 error:
415 *jstat = kstat;
416 s6err ("s1616", *jstat, kpos);
417 goto out;
418
419 out:
420
421 /* Free arrays. */
422
423 if (smatrix != SISL_NULL)
424 freearray (smatrix);
425 if (save != SISL_NULL)
426 freearray (save);
427 if (sright != SISL_NULL)
428 freearray (sright);
429 if (npiv != SISL_NULL)
430 freearray (npiv);
431
432 return;
433 }
434