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: s1331.c,v 1.2 2001-03-19 15:58:45 afr Exp $
45 *
46 */
47
48
49 #define S1331
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1331(double ep[],double eimpli[],int ideg,int ider,double gder[],double gnorm[],int * jstat)55 s1331(double ep[],double eimpli[],int ideg,int ider,
56 double gder[],double gnorm[],int *jstat)
57 #else
58 void s1331(ep,eimpli,ideg,ider,gder,gnorm,jstat)
59 double ep[];
60 double eimpli[];
61 int ideg;
62 int ider;
63 double gder[];
64 double gnorm[];
65 int *jstat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 * PURPOSE : To calculate the position and derivatives up to second
71 * order of a point (and corresponding derivatives) put into
72 * the equation of an implicit surface.
73 *
74 *
75 * INPUT : ep - 0-2 order derivatives of B-spline surface.
76 * For ideg=1,2 and 1001 the sequence is position,
77 * first derivative in first parameter direction,
78 * first derivative in second parameter direction,
79 * (2,0) derivative, (1,1) derivative, (0,2) derivative
80 * and normal. (21 numbers)
81 * For ideg=1003,1004,1005 the second derivatives are followed
82 * by the third derivatives and the normal (33 numbers)
83 * Compatible with output of s1421
84 * eimpli - Description of the implicit surface
85 * ideg - The degree of the implicit surface
86 * ideg=1: Plane
87 * ideg=2; Quadric surface
88 * ideg=1001: Torus surface
89 * ideg=1003: Silhouette line parallel projection
90 * ideg=1004: Silhouette line perspective projection
91 * ideg=1005: Silhouette line circular projection
92 * ider - Derivatives to be calculated
93 * ider=-1: Normal
94 * ider=0: Value + normal
95 * ider=1: Value + first derivatives + normal
96 * ider=2: Value + 1.st + 2.nd + normal
97 * Note if ideg=1003,1004,1005 then ider is set to max(ider,1).
98 *
99 *
100 * OUTPUT :
101 * jstat - status messages
102 * = 0 : ok, curvature radius
103 * < 0 : error
104 *
105 * gder - Specified derivatives in order: value,
106 * 1.st derivatives and 2.nd derivatives
107 * gnorm - Normal of the implicit surface.
108 * When ideg=1003,1004,1005 i.e. for silhouette curves
109 * gnorm = f P - f P which is the direction of
110 * v u u v
111 * the silhouette curve along the surface P.
112 *
113 * METHOD
114 *
115 * USE: This function is only working for 3-D input
116 *
117 * REFERENCES :
118 *-
119 * CALLS : s6scpr,fabs,s1307,s6err
120 *
121 *
122 * WRITTEN BY : Tor Dokken, SI, Oslo , Norway, 10 October 1988
123 * REVISED BY : Mike Floater, SI, Oslo , Norway, 31 January 1991
124 * Added perspective silhouette (ideg=1004) and
125 * circular silhouette (ideg=1005).
126 * Introduced a proper useable definition of gnorm for silhouettes.
127 * It's a little different from the usual cases. See s1331
128 * in order to see how it is used.
129 *
130 *********************************************************************
131 */
132 {
133 int kstat = 0; /* Local status variable */
134 int ki,kj,kl; /* Control variables in loop */
135 int kpos = 0; /* Position of error */
136 int ksize; /* Number of doubles for storage of derivateves
137 and normal vector */
138 int ksizem3; /* ksize - 3 */
139 double *spu; /* Pointer to dP/ds */
140 double *spv; /* Pointer to dP/dt */
141 double *spuu; /* Pointer to ddP/(dsds) */
142 double *spuv; /* Pointer to ddP/(dsdt) */
143 double *spvv; /* Pointer to ddP/(dtdt) */
144
145
146 /* If ideg=1,2 or 1001 then only derivatives up to second order
147 are calculated, then 18 doubles for derivatives and 3 for the
148 normal vector are to be used for calculation of points in the
149 spline surface. For ideg=1003,1004,1005 we have a silhouette curve and
150 derivatives up to the third are to be calculated,
151 thus 30 +3 a total of 33 doubles may be calculated.
152 Note, gnorm in these cases is replaced by f P - f P
153 v u u v
154 so we need ideg to be at least 1. */
155
156 if (ideg == 1003 || ideg == 1004 || ideg == 1005)
157 {
158 ksize = 33;
159 ider=max(ider,1);
160 }
161 else
162 {
163 ksize = 21;
164 }
165
166 ksizem3 = ksize -3;
167
168 /* The calculation of the derivatives of one parameter direction with
169 * respect to the other is dependent on the degree of the implicit equation
170 */
171
172 /* Set local pointers */
173
174 spu = ep + 3;
175 spv = ep + 6;
176 spuu = ep + 9;
177 spuv = ep + 12;
178 spvv = ep + 15;
179
180 if (ideg == 1)
181 {
182 /* First degree implicit geometry.
183 *
184 * Let A = (eimpli[0],eimpli[1],eimpli[2],eimpli[3]) and
185 * Q= (P(s,t),1),
186 * then putting Q into the implicit equation gives
187 *
188 * f(u,v) = A Q
189 *
190 * df dQ
191 * -- = A -- = N P , where N is the normal vector of the plane.
192 * du du u
193 *
194 * df dQ
195 * -- = A -- = N P
196 * dv dv v
197 *
198 * 2 2
199 * d f d Q
200 * -- = A --- = N P
201 * 2 2 uu
202 * du du
203 *
204 * 2 2
205 * d f d Q
206 * -- = A --- = N P
207 * dudv dudv uv
208 *
209 * 2 2
210 * d f d Q
211 * -- = A --- = N P
212 * 2 2 vv
213 * dv dv
214 *
215 */
216 if (ider>=0)
217 {
218 gder[0] = s6scpr(ep,eimpli,3) + eimpli[3];
219 if (ider>=1)
220 {
221 gder[1] = s6scpr(spu,eimpli,3);
222 gder[2] = s6scpr(spv,eimpli,3);
223 if (ider>=2)
224 {
225 gder[3] = s6scpr(spuu,eimpli,3);
226 gder[4] = s6scpr(spuv,eimpli,3);
227 gder[5] = s6scpr(spvv,eimpli,3);
228 }
229 }
230 }
231 memcopy(gnorm,eimpli,3,DOUBLE);
232 }
233 else if (ideg == 2)
234 {
235 double sq[4],sduq[4],sdvq[4]; /*vectors*/
236 double tsum1,tsum2,tsum3; /*Accumulation of matrix product */
237
238 /* Second degree implicit geometry.
239 * Denote the 4x4 matrix representing the implicit surface a A.*
240 *
241 * We can now calculate the u and v derivatives of the point put into
242 * the left hand side of the implicit surface equation:
243 *
244 * T T
245 * f(u,v) = Q A Q = sq Q
246 *
247 * d d T dQ T dQ
248 * - f = -- (Q A Q ) = 2 -- A Q = 2 sq --
249 * du du du du
250 *
251 * d d T dQ T dQ
252 * - f = -- (Q A Q ) = 2 -- A Q = 2 sq --
253 * dv dv dv dv
254 *
255 * 2 2 2 T 2
256 * d d T d Q T dQ dQ dQ dQ
257 * -- f = -- (Q A Q ) = 2 --- A Q + 2 -- A -- = 2(sq--- + sduq--)
258 * 2 2 2 du du 2 du
259 * du du du du
260 *
261 * 2 2 2 T 2
262 * d d T d Q T dQ dQ d Q dQ
263 * -- f = -- (Q A Q ) = 2 ----AQ + 2 --A-- = 2(sq---- + sduq--)
264 * dudv dudv dudv du dv dudv dv
265 *
266 * 2 2 2 T 2
267 * d d T d Q T dQ dQ dQ dQ
268 * -- f = -- (Q A Q ) = 2 --- A Q + 2 -- A -- = 2(sq-- + sdvq--)
269 * 2 2 2 dv dv 2 dv
270 * dv dv dv dv
271 *
272 * dQ dQ
273 * First calculate sq = QA, sduq = --A and sdvq = --A
274 * du dv
275 *
276 */
277 for (ki=0;ki<4;ki++)
278 {
279 tsum1 = eimpli[ki+12];
280 tsum2 = DZERO;
281 tsum3 = DZERO;
282 for (kj=0,kl=ki;kj<3;kj++,kl+=4)
283 {
284 tsum1 += eimpli[kl]*ep[kj];
285 if (ider>=1)
286 {
287 tsum2 += eimpli[kl]*spu[kj];
288 tsum3 += eimpli[kl]*spv[kj];
289 }
290 }
291 sq[ki] = tsum1;
292 sduq[ki] = tsum2;
293 sdvq[ki] = tsum3;
294 }
295
296 /* Make value and partial derivatives */
297
298 if (ider>=0)
299 {
300 gder[0] = s6scpr(sq,ep,3) + sq[3];
301 if(ider>=1)
302 {
303 gder[1] = (double)2.0*s6scpr(sq,spu,3);
304 gder[2] = (double)2.0*s6scpr(sq,spv,3);
305 if (ider>1)
306 {
307 gder[3] = (double)2.0*(s6scpr(sq,spuu,3) +
308 s6scpr(sduq,spu,3));
309 gder[4] = (double)2.0*(s6scpr(sq,spuv,3) +
310 s6scpr(sduq,spv,3));
311 gder[5] = (double)2.0*(s6scpr(sq,spvv,3) +
312 s6scpr(sdvq,spv,3));
313 }
314 }
315 }
316 memcopy(gnorm,sq,3,DOUBLE);
317 }
318
319 else if (ideg == 1001)
320 {
321 double sy[3],sz[3],szu[3],szv[3],szuu[3],szuv[3],szvv[3];/* Derivatives */
322 double tzn,tzun,tzvn,tzuun,tzuvn,tzvvn; /* Scalar products */
323 double tzduz,tzdvz,tzduduz,tzdudvz,tzdvdvz; /* Scalar products */
324 double tduzduz,tdvzdvz,tduzdvz; /* Scalar products */
325 double tlenz,tlenz3,tlenz5; /* Vector lengths */
326 double sp[3],sdup[3],sdvp[3]; /* temporary vectrs*/
327 double sdudup[3],sdudvp[3],sdvdvp[3]; /* temporary vectrs*/
328 double *scentr,*saxis,tbigr,tsmalr; /* Torus descript */
329
330 /* Torus surface. */
331
332 scentr = eimpli;
333 saxis = eimpli + 3;
334 tbigr = *(eimpli+6);
335 tsmalr = *(eimpli+7);
336
337
338 /* Intersection of implicit function of 1 and implicit torus surface
339 * and the the surface. Make the following temporary variables.
340 *
341 * y = p - scentr
342 * z = y - (y saxis) saxis
343 *
344 * Put the point and accompanying derivatives into the implicit
345 * representation of the torus
346 *
347 * 2 2
348 * f(u,v) = (y - R z/sqrt(z z) ) - r
349 *
350 *
351 *
352 * R z (z z )
353 * df R z u R z u
354 * -- = 2(y - ---------)(y - --------- + ----------)
355 * du sqrt(z z) u sqrt(z z) 3
356 * sqrt(z z)
357 *
358 * = 2 sp sdup
359 *
360 *
361 * R z (z z )
362 * df R z v R z v
363 * -- = 2(y - ---------)(y - --------- + ----------)
364 * dv sqrt(z z) v sqrt(z z) 3
365 * sqrt(z z)
366 *
367 * = 2 sp sdvp
368 *
369 *
370 *
371 *
372 * 2 R z (z z )
373 * d f u R z u 2
374 * -- = 2(y - --------- + -----------) +
375 * 2 u sqrt(z z) 3
376 * du sqrt(z z)
377 *
378 *
379 * R z z (z z )
380 * R z uu 2R u u
381 * 2(y - ---------)(y - --------- + ----------- +
382 * sqrt(z z) uu sqrt(z z) 3
383 * sqrt(z z)
384 *
385 *
386 * 2
387 * R z(z z +zz ) z(zz )
388 * u u uu R u
389 * -------------- - 3 --------- )
390 * 3 5
391 * sqrt(z z) sqrt(z z)
392 *
393 * = 2 sdup sdup + sp sdudup
394 *
395 *
396 * 2 R z (z z ) R z (z z )
397 * d f u R z u v R z v
398 * ---- = 2(y - --------- + -----------)(y - --------- + -----------) +
399 * u sqrt(z z) 3 v sqrt(z z) 3
400 * dudv sqrt(z z) sqrt(z z)
401 *
402 *
403 * R z z (z z ) z (z z )
404 * R z uv R u v R v u
405 * 2(y - ---------)(y - --------- + ------------ + ------------ +
406 * sqrt(z z) uv sqrt(z z) 3 3
407 * sqrt(z z) sqrt(z z)
408 *
409 *
410 *
411 * R z(z z +zz ) z(zz )(zz )
412 * v u uv R u v
413 * -------------- - 3 -------------)
414 * 3 5
415 * sqrt(z z) sqrt(z z)
416 *
417 *
418 * = 2 sdup sdvp + sp sdudvp
419 *
420 *
421 * 2 R z (z z )
422 * d f v R z v 2
423 * -- = 2(y - --------- + -----------) +
424 * 2 v sqrt(z z) 3
425 * dv sqrt(z z)
426 *
427 *
428 * R z z (z z )
429 * R z vv 2R v v
430 * 2(y - ---------)(y - --------- + ------------ +
431 * sqrt(z z) vv sqrt(z z) 3
432 * sqrt(z z)
433 *
434 *
435 * 2
436 * R z(z z +zz ) z(zz )
437 * v v vv R v
438 * -------------- - 3 --------- )
439 * 3 5
440 * sqrt(z z) sqrt(z z)
441 *
442 *
443 * = 2 sdvp sdvp + sp sdvdvp
444 *
445 *
446 *
447 * y = (p - scentr) = p
448 * u u u
449 *
450 * y = (p - scentr) = p
451 * v v v
452 *
453 * y = (p - scentr) = p
454 * uu uu uu
455 *
456 * y = (p - scentr) = p
457 * uv uv uv
458 *
459 * y = (p - scentr) = p
460 * vv vv vv
461 *
462 * z = (y - (y N)N) = p - (p N)N N = saxis (Torus axis)
463 * u u u u
464 *
465 * z = (y - (y N)N) = p - (p N)N
466 * v v v v
467 *
468 * z = (y - (y N)N) = p - (p N)N
469 * uu uu uu uu
470 *
471 * z = (y - (y N)N) = p - (p N)N
472 * uv uv uv uv
473 *
474 * z = (y - (y N)N) = p - (p N)N
475 * vv vv vv vv
476 *
477 */
478 for (ki=0;ki<3;ki++)
479 sy[ki] = ep[ki] - scentr[ki];
480
481 tzn = s6scpr(sy ,saxis,3);
482 tzun = s6scpr(spu ,saxis,3);
483 tzvn = s6scpr(spv ,saxis,3);
484 if (ider>1)
485 {
486 tzuun = s6scpr(spuu,saxis,3);
487 tzuvn = s6scpr(spuv,saxis,3);
488 tzvvn = s6scpr(spvv,saxis,3);
489 }
490
491 /* Make z and necessary derivatives of z */
492 for (ki=0;ki<3;ki++)
493 {
494 sz[ki] = sy[ki] - tzn*saxis[ki];
495 szu[ki] = spu[ki] - tzun*saxis[ki];
496 szv[ki] = spv[ki] - tzvn*saxis[ki];
497 if (ider>1)
498 {
499 szuu[ki] = spuu[ki] - tzuun*saxis[ki];
500 szuv[ki] = spuv[ki] - tzuvn*saxis[ki];
501 szvv[ki] = spvv[ki] - tzvvn*saxis[ki];
502 }
503 }
504
505 /* Make a number of necessary scalar products */
506
507 tzduz = s6scpr(sz,szu,3);
508 tzdvz = s6scpr(sz,szv,3);
509 tduzduz = s6scpr(szu,szu,3);
510 tduzdvz = s6scpr(szu,szv,3);
511 tdvzdvz = s6scpr(szv,szv,3);
512 if (ider>1)
513 {
514 tzduduz = s6scpr(sz,szuu,3);
515 tzdudvz = s6scpr(sz,szuv,3);
516 tzdvdvz = s6scpr(sz,szvv,3);
517 }
518
519 /* Find lengt of sz */
520
521 tlenz = s6length(sz,3,&kstat);
522
523 if (kstat<0) goto error;
524 if (DEQUAL(tlenz,DZERO)) tlenz =(double)1.0;
525 tlenz3 = tlenz*tlenz*tlenz;
526 tlenz5 = tlenz3*tlenz*tlenz;
527
528 /* Make a number of necessary vectors */
529
530 for (ki=0;ki<3;ki++)
531 {
532 sp[ki] = sy[ki] - tbigr*sz[ki]/tlenz;
533 if (ider>=1)
534 {
535 sdup[ki] = spu[ki] - tbigr*(szu[ki]/tlenz - sz[ki]*tzduz/tlenz3);
536 sdvp[ki] = spv[ki] - tbigr*(szv[ki]/tlenz - sz[ki]*tzdvz/tlenz3);
537 if (ider>=2)
538 {
539 sdudup[ki] = spuu[ki] -
540 tbigr*(szuu[ki]/tlenz -
541 ((double)2.0*szu[ki]*tzduz+
542 sz[ki]*(tduzduz+tzduduz))/tlenz3 +
543 (double)3.0*sz[ki]*tzduz*tzduz/tlenz5);
544
545 sdudvp[ki] = spuv[ki] -
546 tbigr*(szuv[ki]/tlenz -
547 (szu[ki]*tzdvz+szv[ki]*tzduz+
548 sz[ki]*(tduzdvz+tzdudvz))/tlenz3 +
549 (double)3.0*sz[ki]*tzduz*tzdvz/tlenz5);
550
551 sdvdvp[ki] = spvv[ki] -
552 tbigr*(szvv[ki]/tlenz -
553 ((double)2.0*szv[ki]*tzdvz+
554 sz[ki]*(tdvzdvz+tzdvdvz))/tlenz3 +
555 (double)3.0*sz[ki]*tzdvz*tzdvz/tlenz5);
556 }
557 }
558 }
559
560 /* Make the derivatives */
561 if (ider>=0)
562 {
563 gder[0] = s6scpr(sp,sp,3) - tsmalr*tsmalr;
564 if (ider>=1)
565 {
566 gder[1] = (double)2.0*s6scpr(sp,sdup,3);
567 gder[2] = (double)2.0*s6scpr(sp,sdvp,3);
568 if (ider>=2)
569 {
570 gder[3] = (double)2.0*(s6scpr(sdup,sdup,3)+
571 s6scpr(sp,sdudup,3));
572 gder[4] = (double)2.0*(s6scpr(sdup,sdvp,3)+
573 s6scpr(sp,sdudvp,3));
574 gder[5] = (double)2.0*(s6scpr(sdvp,sdvp,3)+
575 s6scpr(sp,sdvdvp,3));
576 }
577 }
578 }
579 memcopy(gnorm,sp,3,DOUBLE);
580 }
581
582 else if (ideg == 1003)
583
584 {
585
586 /* Silhouette curve, the first three elements of eimpli describes
587 the viewing direction.
588
589 The silhouette line is descibed by the implicit equation, when
590 Q(u,v) is the description of the surface:
591
592 f(u,v) = (Q x Q ) view = 0
593 u v
594
595 f = (Q x Q ) view + (Q x Q ) view
596 u uu v u uv
597
598
599 f = (Q x Q ) view + (Q x Q ) view
600 v uv v u vv
601
602
603
604 f = (Q x Q ) view + 2 (Q x Q ) view + (Q x Q ) view
605 uu uuu v uu uv u uuv
606
607
608 f = (Q x Q )view + (Q x Q )view + (Q x Q )view + (Q xQ)view
609 uv uuv v uu vv uv uv u uvv
610
611
612 f = (Q x Q ) view + 2 (Q x Q ) view + (Q x Q ) view
613 vv uvv v uv vv u vvv
614
615
616 Note that Q x Q has zero length.
617 uv uv
618 */
619 double *spuuu = ep + 18;
620 double *spuuv = ep + 21;
621 double *spuvv = ep + 24;
622 double *spvvv = ep + 27;
623 double sdum1[3],sdum2[3],sdum3[3];
624
625 if (ider>=0)
626 {
627 s6crss(spu,spv,sdum1);
628 gder[0] = s6scpr(sdum1,eimpli,3);
629 if (ider>=1)
630 {
631 s6crss(spuu,spv,sdum1);
632 s6crss(spu,spuv,sdum2);
633 gder[1] = s6scpr(sdum1,eimpli,3) + s6scpr(sdum2,eimpli,3);
634 s6crss(spuv,spv,sdum1);
635 s6crss(spu,spvv,sdum2);
636 gder[2] = s6scpr(sdum1,eimpli,3) + s6scpr(sdum2,eimpli,3);
637 if (ider>=2)
638 {
639 s6crss(spuuu,spv, sdum1);
640 s6crss(spuu, spuv, sdum2);
641 s6crss(spu, spuuv,sdum3);
642 gder[3] = s6scpr(sdum1,eimpli,3) +
643 (double)2.0*s6scpr(sdum2,eimpli,3)
644 + s6scpr(sdum3,eimpli,3);;
645 s6crss(spuuv,spv, sdum1);
646 s6crss(spuu, spvv, sdum2);
647 s6crss(spu, spuvv,sdum3);
648 gder[4] = s6scpr(sdum1,eimpli,3) + s6scpr(sdum2,eimpli,3)
649 + s6scpr(sdum3,eimpli,3);;
650 s6crss(spuvv,spv, sdum1);
651 s6crss(spuv, spvv, sdum2);
652 s6crss(spu, spvvv,sdum3);
653 gder[5] = s6scpr(sdum1,eimpli,3) +
654 (double)2.0*s6scpr(sdum2,eimpli,3)
655 + s6scpr(sdum3,eimpli,3);;
656 }
657 }
658 }
659 }
660
661 else if (ideg == 1004)
662
663 {
664
665 /* Perspective silhouette curve, the first three elements of eimpli
666 describe the eye point E.
667
668 The silhouette line is descibed by the implicit equation, when
669 P(u,v) is the description of the surface and
670 N(u,v) = P x P is the normal to the surface:
671 u v
672
673 f(u,v) = N(u,v) . (P(u,v) - E) = 0
674
675 Differentiating N gives:
676
677 N = P x P + P x P
678 u uu v u uv
679
680 N = P x P + P x P
681 v uv v u vv
682
683 N = P x P + 2 P x P + P x P
684 uu uuu v uu uv u uuv
685
686 N = P x P + P x P + P x P
687 uv uuu v uu vv u uuv
688
689 N = P x P + 2 P x P + P x P
690 vv uvv v uv vv u vvv
691
692
693
694 Differentiating f gives:
695
696 f = N . (P - E)
697 u u
698
699 f = N . (P - E)
700 v v
701
702 f = N . (P - E) + N . P
703 uu uu u u
704
705 f = N . (P - E) + N . P
706 uv uv u v
707
708 f = N . (P - E) + N . P
709 vv vv v v
710
711
712 */
713 double *spuuu = ep + 18;
714 double *spuuv = ep + 21;
715 double *spuvv = ep + 24;
716 double *spvvv = ep + 27;
717 double norm[3],normu[3],normv[3],normuu[3],normuv[3],normvv[3];
718 double cprod1[3],cprod2[3],cprod3[3];
719 double pediff[3];
720
721
722 if (ider>=0)
723 {
724 s6diff(ep,eimpli,3,pediff);
725
726 s6crss(spu,spv,norm);
727 gder[0] = s6scpr(norm,pediff,3);
728
729 if (ider>=1)
730 {
731 s6crss(spuu,spv,cprod1);
732 s6crss(spu,spuv,cprod2);
733 for (ki=0;ki<3;ki++)
734 {
735 normu[ki] = cprod1[ki] + cprod2[ki];
736 }
737 s6crss(spuv,spv,cprod1);
738 s6crss(spu,spvv,cprod2);
739 for (ki=0;ki<3;ki++)
740 {
741 normv[ki] = cprod1[ki] + cprod2[ki];
742 }
743
744 gder[1] = s6scpr(normu,pediff,3);
745 gder[2] = s6scpr(normv,pediff,3);
746
747 if (ider>=2)
748 {
749 s6crss(spuuu,spv,cprod1);
750 s6crss(spuu,spuv,cprod2);
751 s6crss(spu,spuuv,cprod3);
752 for (ki=0;ki<3;ki++)
753 {
754 normuu[ki] = cprod1[ki] + 2.0 * cprod2[ki] + cprod3[ki];
755 }
756
757 /* Note that spuv x spuv = 0 in the formula for normuv */
758 s6crss(spuuv,spv,cprod1);
759 s6crss(spuu,spvv,cprod2);
760 s6crss(spu,spuvv,cprod3);
761 for (ki=0;ki<3;ki++)
762 {
763 normuv[ki] = cprod1[ki] + cprod2[ki] + cprod3[ki];
764 }
765
766 s6crss(spuvv,spv,cprod1);
767 s6crss(spuv,spvv,cprod2);
768 s6crss(spu,spvvv,cprod3);
769 for (ki=0;ki<3;ki++)
770 {
771 normvv[ki] = cprod1[ki] + 2.0 * cprod2[ki] + cprod3[ki];
772 }
773
774 gder[3] = s6scpr(normuu,pediff,3) + s6scpr(normu,spu,3);
775 gder[4] = s6scpr(normuv,pediff,3) + s6scpr(normu,spv,3);
776 gder[5] = s6scpr(normvv,pediff,3) + s6scpr(normv,spv,3);
777
778 }
779 }
780 }
781 }
782
783 else if (ideg == 1005)
784
785 {
786
787 /* Perspective silhouette curve, the first three elements of eimpli
788 describe Q, the next three describe B.
789
790 The silhouette line is descibed by the implicit equation, when
791 P(u,v) is the description of the surface and
792 N(u,v) = P x P is the normal to the surface:
793 u v
794
795 f(u,v) = N(u,v) x (P(u,v) - Q) . B = 0
796
797 Differentiating N gives:
798
799 N = P x P + P x P
800 u uu v u uv
801
802 N = P x P + P x P
803 v uv v u vv
804
805 N = P x P + 2 P x P + P x P
806 uu uuu v uu uv u uuv
807
808 N = P x P + P x P + P x P
809 uv uuu v uu vv u uuv
810
811 N = P x P + 2 P x P + P x P
812 vv uvv v uv vv u vvv
813
814
815
816 Differentiating f gives:
817
818 f = { N x (P - Q) + N x P } . B
819 u u u
820
821 f = { N x (P - Q) + N x P } . B
822 v v v
823
824 f = { N x (P - Q) + 2 N x P + N x P } . B
825 uu uu u u uu
826
827 f = { N x (P - Q) + N x P + N x P + N x P } . B
828 uv uv u v v u uv
829
830 f = { N x (P - Q) + 2 N x P + N x P } . B
831 vv vv v v vv
832
833
834 */
835 double *spuuu = ep + 18;
836 double *spuuv = ep + 21;
837 double *spuvv = ep + 24;
838 double *spvvv = ep + 27;
839 double *bvec = eimpli + 3;
840 double norm[3],normu[3],normv[3],normuu[3],normuv[3],normvv[3];
841 double cprod1[3],cprod2[3],cprod3[3],cprod4[3];
842 double pqdiff[3],sum[3];
843
844
845 if (ider>=0)
846 {
847 s6diff(ep,eimpli,3,pqdiff);
848
849 s6crss(spu,spv,norm);
850 s6crss(norm,pqdiff,cprod1);
851 gder[0] = s6scpr(cprod1,bvec,3);
852
853 if (ider>=1)
854 {
855 s6crss(spuu,spv,cprod1);
856 s6crss(spu,spuv,cprod2);
857 for (ki=0;ki<3;ki++)
858 {
859 normu[ki] = cprod1[ki] + cprod2[ki];
860 }
861 s6crss(spuv,spv,cprod1);
862 s6crss(spu,spvv,cprod2);
863 for (ki=0;ki<3;ki++)
864 {
865 normv[ki] = cprod1[ki] + cprod2[ki];
866 }
867
868 s6crss(normu,pqdiff,cprod1);
869 s6crss(norm,spu,cprod2);
870 for (ki=0;ki<3;ki++)
871 {
872 sum[ki] = cprod1[ki] + cprod2[ki];
873 }
874 gder[1] = s6scpr(sum,bvec,3);
875
876 s6crss(normv,pqdiff,cprod1);
877 s6crss(norm,spv,cprod2);
878 for (ki=0;ki<3;ki++)
879 {
880 sum[ki] = cprod1[ki] + cprod2[ki];
881 }
882 gder[2] = s6scpr(sum,bvec,3);
883
884 if (ider>=2)
885 {
886 s6crss(spuuu,spv,cprod1);
887 s6crss(spuu,spuv,cprod2);
888 s6crss(spu,spuuv,cprod3);
889 for (ki=0;ki<3;ki++)
890 {
891 normuu[ki] = cprod1[ki] + 2.0 * cprod2[ki] + cprod3[ki];
892 }
893
894 /* Note that spuv x spuv = 0 in the formula for normuv */
895 s6crss(spuuv,spv,cprod1);
896 s6crss(spuu,spvv,cprod2);
897 s6crss(spu,spuvv,cprod3);
898 for (ki=0;ki<3;ki++)
899 {
900 normuv[ki] = cprod1[ki] + cprod2[ki] + cprod3[ki];
901 }
902
903 s6crss(spuvv,spv,cprod1);
904 s6crss(spuv,spvv,cprod2);
905 s6crss(spu,spvvv,cprod3);
906 for (ki=0;ki<3;ki++)
907 {
908 normvv[ki] = cprod1[ki] + 2.0 * cprod2[ki] + cprod3[ki];
909 }
910
911 s6crss(normuu,pqdiff,cprod1);
912 s6crss(normu,spu,cprod2);
913 s6crss(norm,spuu,cprod3);
914 for (ki=0;ki<3;ki++)
915 {
916 sum[ki] = cprod1[ki] + 2.0 * cprod2[ki] + cprod3[ki];
917 }
918 gder[3] = s6scpr(sum,bvec,3);
919
920 s6crss(normuv,pqdiff,cprod1);
921 s6crss(normu,spv,cprod2);
922 s6crss(normv,spu,cprod3);
923 s6crss(norm,spuv,cprod4);
924 for (ki=0;ki<3;ki++)
925 {
926 sum[ki] = cprod1[ki] + cprod2[ki] + cprod3[ki] + cprod4[ki];
927 }
928 gder[4] = s6scpr(sum,bvec,3);
929
930 s6crss(normvv,pqdiff,cprod1);
931 s6crss(normv,spv,cprod2);
932 s6crss(norm,spvv,cprod3);
933 for (ki=0;ki<3;ki++)
934 {
935 sum[ki] = cprod1[ki] + 2.0 * cprod2[ki] + cprod3[ki];
936 }
937 gder[5] = s6scpr(sum,bvec,3);
938
939 }
940 }
941 }
942 }
943
944 if(ideg == 1003 || ideg == 1004 || ideg == 1005)
945 {
946 /* The normal vector has no meaning in this case
947 so we calculate the direction D(u,v) of the solution curve
948 on the surface P(u,v)
949 of f(u,v) = N(u,v) . d = 0 (ideg = 1003)
950 or f(u,v) = N(u,v) . (P(u,v) - E) = 0 (ideg = 1004)
951 or f(u,v) = N(u,v) x (P(u,v) - Q) . B = 0 (ideg = 1005)
952 (through the present solution point)
953 directly. This direction D is in fact f P - f P
954 v u u v
955 where P(u,v) is the parameterised
956 surface, N is the normal P x P , and d is the
957 u v
958 view direction. The direction D is used by s1331. */
959
960 for(ki=0; ki<3; ki++)
961 {
962 gnorm[ki] = gder[2] * spu[ki] - gder[1] * spv[ki];
963 }
964 (void)s6norm(gnorm,3,gnorm,&kstat);
965 }
966
967 /* Everyting is ok */
968
969 *jstat = 0;
970 goto out;
971
972
973 /* Error in lower level function */
974 error:
975 *jstat = kstat;
976 s6err("s1331",*jstat,kpos);
977 goto out;
978
979 out:
980 return;
981 }
982