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: s1377.c,v 1.3 2001-03-19 15:58:48 afr Exp $
45 *
46 */
47
48
49 #define S1377
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1377(SISLCurve * pcurv,double econic[],int ideg,int idim,SISLCurve ** rcurv,int * jstat)55 s1377 (SISLCurve * pcurv, double econic[], int ideg, int idim,
56 SISLCurve ** rcurv, int *jstat)
57 #else
58 void
59 s1377 (pcurv, econic, ideg, idim, rcurv, jstat)
60 SISLCurve *pcurv;
61 double econic[];
62 int ideg;
63 int idim;
64 SISLCurve **rcurv;
65 int *jstat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 * PURPOSE : To put a curve description into the descripiton of
71 * a torus surface described by the input array econic.
72 *
73 *
74 * INPUT : pcurv - Pointer to input curve
75 * econic - Description of torus
76 * ideg - Type of conic: torus: ideg=1001
77 * idim - Dimension of object space
78 *
79 *
80 * OUTPUT : rcurv - The resulting curve
81 * jstat - status messages
82 * > 0 : warning
83 * = 0 : ok
84 * < 0 : error
85 *
86 *
87 * METHOD : We first make the appropriate knot vector, then we calulate
88 * parametervalues for the interpolation, then the appropriate
89 * values of the curve put into the conic equation are found,
90 * and at last the curve is interpolated.
91 *
92 * REFERENCES :
93 *
94 * CALLS : s1376,s1890,s1221,s6scpr,s1891,s6err.
95 *
96 * WRITTEN BY : Tor Dokken, SI, 1988-11
97 * REVISED BY : Mike Floater, SI, 1991-04 for a rational curve.
98 *
99 *********************************************************************
100 */
101 {
102 int ikind; /* Type of curve pcurv is */
103 int kn; /* Number of vertices of pcurv */
104 int kk; /* Order in pcurv */
105 int kjkk; /* Order of interpolated basis */
106 int kjkn; /* Number of vertices in interpolated basis */
107 int kdim; /* Number of dimensions in pcurv */
108 int kstat; /* Local status variable */
109 int kpos = 0; /* Position indicator for errors */
110 int kzero = 0; /* Value 0 needed in call s1891 */
111 int kone = 1; /* Value 1 needed in call s1891 */
112 int cuopen; /* Open/Closed flag */
113 int ki, kj; /* Loop control variable */
114 int kleft = 0; /* Pointer into knot vector */
115 int *der = SISL_NULL; /* Derivate indicators. */
116 double *st = SISL_NULL; /* First knot vector is pcurv */
117 double *scentr = econic; /* Center of torus */
118 double *saxis = econic + 3; /* Axis of torus */
119 double tbigr = *(econic + 6); /* Big radius of torus */
120 double tsmalr = *(econic + 7);/* Small radius of torus */
121 double tbigr2 = tbigr * tbigr;/* Square of big radius */
122 double tdiffr2 = tbigr2 - tsmalr * tsmalr; /* Difference of square of radii*/
123 double *sval1 = SISL_NULL; /* Array of values of curve put into torus eq. */
124 double *sval2 = SISL_NULL;
125 double *sgt = SISL_NULL; /* Knot vector of curve put into torus surface */
126 double sy[3]; /* Difference between point and torus center */
127 double tzn; /* Projection of sy onto torus axis */
128 double tyy; /* Square of length of sy */
129 double tzz; /* Square of length of sz */
130 double ty; /* Component of sy */
131 double tz; /* Component of sz */
132 double sder[4]; /* Point on the curve */
133 double ww; /* the weight of sder squared if pcurv is rational */
134 double *par = SISL_NULL;
135 SISLCurve *tempcurv = SISL_NULL; /* only used for rational curves */
136
137 *jstat = 0;
138 if (idim != pcurv->idim) goto err104;
139 if (ideg != 1001) goto err200;
140
141 /* Make local pointers. */
142
143 kn = pcurv->in;
144 kk = pcurv->ik;
145 kdim = pcurv->idim;
146 st = pcurv->et;
147 ikind = pcurv->ikind;
148
149 if (ikind == 2 || ikind == 4)
150 {
151 tempcurv = newCurve (kn, kk, st, pcurv->rcoef, ikind - 1, kdim + 1, 0);
152 if (tempcurv == SISL_NULL)
153 goto err171;
154 tempcurv->cuopen = pcurv->cuopen;
155 }
156 else
157 {
158 tempcurv = pcurv;
159 }
160
161
162 /* Test input. */
163
164 if (kdim != 3)
165 goto err104;
166
167
168 /* Make description of knot array for interpolation. */
169
170 s1376 (st, kn, kk, &sgt, &kjkn, &kjkk, &kstat);
171 if (kstat < 0)
172 goto error;
173
174
175 /* Make parameter values and derivative indicators. */
176
177 s1890 (sgt, kjkk, kjkn, &par, &der, &kstat);
178 if (kstat < 0)
179 goto error;
180
181
182 /* Allocate array for values of curve put into torus equation. */
183
184 sval1 = newarray (kjkn, DOUBLE);
185 if (sval1 == SISL_NULL)
186 goto err101;
187
188
189 /* Calculate values to be interpolated. */
190
191 for (ki = 0; ki < kjkn; ki++)
192 {
193 /* Calculate values on 3-D curve. */
194
195 s1221 (tempcurv, 0, par[ki], &kleft, sder, &kstat);
196 if (kstat < 0)
197 goto error;
198
199 /*
200 * The calculation of a point on the torus surface can be done in the
201 * following way.
202 *
203 * y = p - scentr
204 * z = y - (y saxis) saxis
205 *
206 * The equation of the torus can be written
207 *
208 * 2 2
209 * (y - R z/sqrt(z z) ) - r = 0
210 *
211 *
212 * or by eliminating the square root:
213 *
214 * f =
215 *
216 * 2 2 2 2 2 2 2
217 * (yy) + 2 (yy)(R -r ) - 4R zz + (R -r ) = 0
218 *
219 *
220 *
221 * or in 4-D homogeneous coordinates:
222 *
223 *
224 * f =
225 *
226 * 2 2 2 2 2 2 4 2 2 2
227 * (yy) + 2 w (yy)(R -r ) - 4w R zz + w (R -r ) = 0
228 *
229 * where y = T - w*scentr, p=T/w
230 *
231 * We thus need to calculate yy and zz:
232 */
233
234 if (ikind == 2 || ikind == 4)
235 {
236 for (kj = 0; kj < 3; kj++)
237 sy[kj] = sder[kj] - sder[3] * scentr[kj];
238 ww = sder[3] * sder[3];
239 }
240 else
241 {
242 for (kj = 0; kj < 3; kj++)
243 sy[kj] = sder[kj] - scentr[kj];
244 ww = (double) 1.0;
245 }
246
247 tzn = s6scpr (sy, saxis, 3);
248
249 tyy = (double) 0.0;
250 tzz = (double) 0.0;
251
252
253 /* Make z and necessary derivatives of z */
254
255 for (kj = 0; kj < 3; kj++)
256 {
257 ty = sy[kj];
258 tz = ty - tzn * saxis[kj];
259 tyy += ty * ty;
260 tzz += tz * tz;
261 }
262
263 /*
264 * 2 2 2
265 * Now tyy = yy and tzz = zz, tbigr2 = R ,tdiffr2 = R - r
266 * --------------------------------------------------------
267 */
268
269 sval1[ki] = tyy * tyy + ((double) 2.0 * ww * tyy
270 + ww * ww * tdiffr2) * tdiffr2
271 - (double) 4.0 *ww * tbigr2 * tzz;
272 }
273
274 cuopen = TRUE;
275
276 s1891 (par, sval1, kone, kjkn, kone, der, cuopen, sgt, &sval2, &kjkn, kjkk,
277 kzero, kzero, &kstat);
278 if (kstat < 0)
279 goto error;
280
281 *rcurv = SISL_NULL;
282 *rcurv = newCurve (kjkn, kjkk, sgt, sval2, 1, 1, 1);
283 if (*rcurv == SISL_NULL)
284 goto err171;
285 (*rcurv)->cuopen = pcurv->cuopen;
286
287
288 /* Ok ! */
289
290 goto out;
291
292
293 /* Error in lower level function */
294
295 error:
296 *jstat = kstat;
297 s6err ("s1377", *jstat, kpos);
298 goto out;
299
300 /* Error in space allocation */
301
302 err101:
303 *jstat = -101;
304 s6err ("s1377", *jstat, kpos);
305 goto out;
306
307 /* Dimension not equal to 3 or conflicting dimensions */
308
309 err104:
310 *jstat = -104;
311 s6err ("s1377", *jstat, kpos);
312 goto out;
313
314 /* Could not create curve. */
315
316 err171:
317 *jstat = -171;
318 s6err ("s1377", *jstat, kpos);
319 goto out;
320
321 /* Wrong implicit type (ideg). */
322
323 err200:
324 *jstat = -200;
325 s6err ("s1377", *jstat, kpos);
326 goto out;
327
328 out:
329
330 /* Release allocated arrays */
331
332 if (sgt != SISL_NULL)
333 freearray (sgt);
334 if (par != SISL_NULL)
335 freearray(par);
336 if (der != SISL_NULL)
337 freearray(der);
338 if (sval1 != SISL_NULL)
339 freearray (sval1);
340 if (sval2 != SISL_NULL)
341 freearray (sval2);
342 if ((ikind == 2 || ikind == 4) && (tempcurv != SISL_NULL))
343 freeCurve (tempcurv);
344
345 return;
346 }
347