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: s1500.c,v 1.2 2001-03-19 15:58:49 afr Exp $
45 *
46 */
47
48
49 #define S1500
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1500(double base[],double norm[],double axisA[],double alpha,double ratio,int idim,int inumb,double carray[],int * jstat)55 s1500(double base[],double norm[],double axisA[],double alpha,
56 double ratio,int idim,int inumb,double carray[],int *jstat)
57 #else
58 void s1500(base,norm,axisA,alpha,ratio,idim,inumb,carray,jstat)
59 double base[];
60 double norm[];
61 double axisA[];
62 double alpha;
63 double ratio;
64 int idim;
65 int inumb;
66 double carray[];
67 int *jstat;
68 #endif
69 /*
70 *********************************************************************
71 *
72 * PURPOSE : To make a matrix of dimension (idim+1)x(idim+1)
73 * describing an elliptic cone as an implicit function.
74 *
75 *
76 * INPUT : base - The base point of the cone
77 * (the centre of the base ellipse)
78 * norm - Direction of cone axis
79 * axisA - One of the axes of the ellipse
80 * alpha - The opening angle of axisA
81 * ratio - The ratio of axisA to axisB
82 * idim - The dimension of the space the cylinder lies in
83 * inumb - The number of copies that are to be made of the
84 * matrix.
85 *
86 *
87 *
88 * OUTPUT : carray - The description of the elliptic cone. Outside
89 * this function the space for this array must be
90 * allocated. The need is (idim+1)*(idim+1)*inumb
91 * dimension 4x4 (xinarr)
92 * jstat - status messages
93 * > 0 : warning
94 * = 0 : ok
95 * < 0 : error
96 *
97 *
98 * METHOD :
99 *
100 * If the base point of the cone is denoted p=(p1,p2,p3),
101 * the axis vector of the cone is denoted n=(n1,n2,n3),
102 * the direction vector of one ellipse axis is denoted A=(a1,a2,a3),
103 * the second B = (n X A) / |n x A| is B=(b1,b2,b3),
104 * alpha is the opening angle of the cone at the axis A,
105 * then the matrix describing
106 * the cone is A:
107 *
108 *
109 * WHERE
110 * A11 = cos2*a1*a1+cos2*r*b1*b1-sin2*n1*n1
111 * A21 = A12 = cos2*a1*a2+cos2*r*b1*b2-sin2*n1*n2
112 * A31 = A13 = cos2*a1*a3+cos2*r*b1*b3-sin2*n1*n3
113 * A41 = A14 = -cos2*(A.p)*a1-cos2*r*(B.p)*b1+a*cos*sin*n1+sin2*(n.p)*n1
114 * A22 = cos2*a2*a2+cos2*r*b2*b2-sin2*n2*n2
115 * A32 = A23 = cos2*a2*a3+cos2*r*b2*b3-sin2*n2*n3
116 * A42 = A24 = -cos2*(A.p)*a2-cos2*r*(B.p)*b2+a*cos*sin*n2+sin2*(n.p)*n2
117 * A33 = cos2*a3*a3+cos2*r*b3*b3-sin2*n3*n3
118 * A43 = A34 = -cos2*(A.p)*a3-cos2*r*(B.p)*b3+a*cos*sin*n3+sin2*(n.p)*n3
119 * A44 = cos2*(A.p)**2+cos2*r*(B.p)**2-(cos*a+sin*(n.p))**2
120 *
121 * aND
122 *
123 * cos=cos(alpha), sin=sin(alpha), cos2=cos*cos, sin2=sin*sin
124 * a=|A|, r=(|A|/|B|)**2,
125 * A.p=a1*p1+a2*p2+a3*p3 (dot product), similarly for B.p, n.p,
126 * and B=n*A/|nxA| (cross product).
127 *
128 *
129 * REFERENCES :
130 *
131 *-
132 * CALLS :
133 *
134 * WRITTEN BY : Mike Floater, SI, Oslo, Norway, 16-October-1988
135 *
136 *********************************************************************
137 */
138 {
139 int kdimp1; /* Dimension of matrix kdimp1 = idim + 1 */
140 int kstop; /* Stop condition for for loop */
141 int ki,kj,kl; /* Running variables in loop */
142 int kpos=0; /* Position of error */
143 int kstat; /* Local status variable */
144 double a1,a2,a3; /* Coordinates of A vector */
145 double b1,b2,b3; /* Coordinates of B vector */
146 double n1,n2,n3; /* Coordinates of ndirec vector */
147 double temp; /* Temporary storage variable */
148 double cosT; /* cos(alpha) */
149 double cosT2; /* The square of cos(alpha) */
150 double sinT; /* sin(alpha) */
151 double sinT2; /* The square of sin(alpha) */
152 double cosTsinT; /* cos(alpha)*sin(alpha) */
153 double ndirec[3]; /* Normalized direction of cone axis */
154 double A[3]; /* Normalized elliptic axisA */
155 double B[3]; /* Normalized elliptic axisB */
156 double asize; /* Length of axisA */
157 double r; /* ratio*ratio */
158 double Adotp; /* scalar product of A and base */
159 double Bdotp; /* scalar product of B and base */
160 double ndotp; /* scalar product of n and base */
161
162
163 /* Test i legal input */
164 if (inumb <1 ) inumb = 1;
165 if (idim != 3 ) goto err104;
166
167 kdimp1 = idim + 1;
168 kstop = kdimp1*kdimp1;
169
170 for (ki=0;ki<kstop;ki++)
171 {
172 carray[ki] = DZERO;
173 }
174
175 /* Normalise direction vector of cone axis */
176
177 (void)s6norm(norm,idim,ndirec,&kstat);
178
179 /* Test if norm degenerate */
180 if(kstat == 0) goto err174;
181
182
183 /* Normalise elliptic axisA and find its length */
184
185 asize=s6norm(axisA,idim,A,&kstat);
186
187 /* Test if norm degenerate */
188 if(kstat == 0) goto err174;
189
190 /* Calculate the other elliptic axis from ndirec and A */
191 /* Note that norm and axisA are assumed to be orthogonal */
192
193 (void)s6crss(ndirec,A,B);
194
195 /* Calculate cosine and sine of angle alpha */
196
197 cosT=cos(alpha);
198 sinT=sin(alpha);
199 cosT2=cosT*cosT;
200 sinT2=sinT*sinT;
201 cosTsinT=cosT*sinT;
202
203 /* Find the square of the ratio of lengths of axisA and axisB */
204
205 r=ratio*ratio;
206
207
208 /* set up coordinates of the normalised axes */
209
210 a1 = A[0];
211 a2 = A[1];
212 a3 = A[2];
213 b1 = B[0];
214 b2 = B[1];
215 b3 = B[2];
216 n1 = ndirec[0];
217 n2 = ndirec[1];
218 n3 = ndirec[2];
219
220 /* set up dot (scalar) products of axes and base point */
221
222 Adotp = s6scpr(A,base,idim);
223 Bdotp = s6scpr(B,base,idim);
224 ndotp = s6scpr(ndirec,base,idim);
225
226
227 /* Make element (1,1) */
228
229 carray[0] = cosT2*(a1*a1+r*b1*b1)-sinT2*n1*n1;
230
231 /* Make element (1,1) */
232
233 carray[5] = cosT2*(a2*a2+r*b2*b2)-sinT2*n2*n2;
234
235 /* Make element (1,1) */
236
237 carray[10] = cosT2*(a3*a3+r*b3*b3)-sinT2*n3*n3;
238
239 /* Make element (1,2) and (2,1) */
240
241 temp = cosT2*(a1*a2+r*b1*b2)-sinT2*n1*n2;
242 carray[1] = temp;
243 carray[4] = temp;
244
245
246 /* Make element (1,3) and (3,1) */
247
248 temp = cosT2*(a1*a3+r*b1*b3)-sinT2*n1*n3;
249 carray[2] = temp;
250 carray[8] = temp;
251
252
253 /* Make element (2,3) and (3,2) */
254
255 temp = cosT2*(a2*a3+r*b2*b3)-sinT2*n2*n3;
256 carray[6] = temp;
257 carray[9] = temp;
258
259
260
261 /* Make element (1,4) and (4,1) */
262
263 temp = -cosT2*(Adotp*a1+r*Bdotp*b1)+(asize*cosTsinT+sinT2*ndotp)*n1;
264 carray[3] = temp;
265 carray[12] = temp;
266
267
268 /* Make element (2,4) and (4,2) */
269
270 temp = -cosT2*(Adotp*a2+r*Bdotp*b2)+(asize*cosTsinT+sinT2*ndotp)*n2;
271 carray[7] = temp;
272 carray[13] = temp;
273
274
275 /* Make element (3,4) and (4,3) */
276
277 temp = -cosT2*(Adotp*a3+r*Bdotp*b3)+(asize*cosTsinT+sinT2*ndotp)*n3;
278 carray[11] = temp;
279 carray[14] = temp;
280
281 /* Make element (4,4) */
282
283 temp = cosT*asize+sinT*ndotp;
284 carray[15] = cosT2*(Adotp*Adotp+r*Bdotp*Bdotp)-temp*temp;
285
286
287 /* Make extra copies of cone */
288
289 kj = kstop;
290 for (ki=1;ki<inumb;ki++)
291 {
292 for (kl=0;kl<kstop;kl++,kj++)
293 {
294 carray[kj] = carray[kl];
295 }
296 }
297
298 *jstat = 0;
299 goto out;
300
301 /* Dimension less than 1 */
302 err104: *jstat = -104;
303 s6err("s1500",*jstat,kpos);
304 goto out;
305
306 /* norm axis or elliptic axis is zero */
307 err174: *jstat = -174;
308 s6err("s1500",*jstat,kpos);
309 goto out;
310
311 out:
312 return;
313 }
314