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: s1323.c,v 1.2 2001-03-19 15:58:44 afr Exp $
45 *
46 */
47
48
49 #define S1323
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1323(double etop[],double eaxis[],double econe[],int idim,int inumb,double carray[],int * jstat)55 s1323(double etop[],double eaxis[],double econe[],int idim,
56 int inumb,double carray[],int *jstat)
57 #else
58 void s1323(etop,eaxis,econe,idim,inumb,carray,jstat)
59 double etop[];
60 double eaxis[];
61 double econe[];
62 int idim;
63 int inumb;
64 double carray[];
65 int *jstat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 * PURPOSE : To make a matrix of dimension (idim+1)x(idim+1)
71 * describing a cone as an implicit function.
72 *
73 *
74 * INPUT : etop - The top point of the cone
75 * edirec - Direction of cylinder axis
76 * econe - A point on the cone surface different from the
77 * top point
78 * idim - The dimension of the space the cylinder lies
79 * inumb - The number of copies that are to be made of the
80 * matrix.
81 *
82 *
83 *
84 * OUTPUT : carray - The description of the cone. Outside
85 * this function the space for this array must be
86 * allocated. The need is (idim+1)*(idim+1)*inumb
87 * dimension 4x4 (xinarr)
88 * jstat - status messages
89 * > 0 : warning
90 * = 0 : ok
91 * < 0 : error
92 *
93 *
94 * METHOD :
95 *
96 * If the top point of the cone is denoted (X0,Y0,Z0), the direction
97 * vector of the cone axis is denoted (WX,WY,WZ) and COS(T) is
98 * the cosine of the opining angle of the cone, the matrix describing
99 * the cone is:
100 *
101 * I- -I
102 * I WX*WX -WX*WY -WX*WZ I
103 * I 1 - ---------- ----------- ----------- A I
104 * I (COS(T)**2) (COS(T)**2) (COS(T)**2) I
105 * I I
106 * I -WX*WY WY*WY -WY*WZ I
107 * I ----------- 1 - ----------- ----------- B I
108 * I (COS(T)**2) (COS(T)**2) (COS(T)**2) I
109 * I I
110 * I -WX*WZ -WY*WZ WZ*WZ I
111 * I ----------- ----------- 1 - ----------- C I
112 * I (COS(T)**2) (COS(T)**2) (COS(T)**2) I
113 * I I
114 * I A B * D I
115 * I- -I
116 *
117 * WHERE
118 * A = (X0*WX*WX+WX*(Y0*WY+Z0*WZ))/(COS(T)**2)-X0
119 * B = (Y0*WY*WY+WY*(Z0*WZ+X0*WX))/(COS(T)**2)-Y0
120 * C = (Z0*WZ*WZ+WZ*(X0*WX+Y0*WY))/(COS(T)**2)-Z0
121 * D = X0*X0+Y0*Y0+Z0*Z0-(X0*X0*WX*WY+Y0*Y0*WY*WY+Z0*Z0*WZ*WZ
122 * +2*X0*Y0*WX*WY+2*Y0*Z0*WY*WZ+2*Z0*X0*WZ*WX)/(COS(T)**2)
123 *
124 * The matrix is described in the following way: (X0,Y0,Z0) point
125 * on cylinder axis, (WX,WY,WZ) direction of cylinder axis and
126 * R radius of cylinder:
127 *
128 * I- -I
129 * I 1-WX*WX -WX*WY -WX*WZ A I
130 * I I
131 * I -WX*WY 1-WY*WY -WY*WZ B I
132 * I I
133 * I -WX*WZ -WY*WZ 1-WZ*WZ * I
134 * I I
135 * I A B C D I
136 * I- -I
137 *
138 * where
139 *
140 * A = X0*(WX*WX-1)+WX*(Y0*WY+Z0*WZ)
141 * B = Y0*(WY*WY-1)+WY*(Z0*WZ+X0*WX)
142 * C = Z0*(WZ*WZ-1)+WZ*(X0*WX+Y0*WY)
143 * D = X0*X0+Y0*Y0+Z0*Z0-X0*X0*WX*WX-Y0*Y0*WY*WY-Z0*Z0*WZ*WZ
144 * -2*X0*Y0*WX*WY-2*Y0*Z0*WY*WZ-2*Z0*X0*WZ*WX-R*R
145 *
146 * REFERENCES :
147 *
148 *-
149 * CALLS :
150 *
151 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway, 28-June-1988
152 *
153 *********************************************************************
154 */
155 {
156 int kdimp1; /* Dimension of matrix kdimp1 = idim + 1 */
157 int kdimp2; /* idim + 2 */
158 int kstop; /* Stop condition for for loop */
159 int ki,kj,kl; /* Running variables in loop */
160 int kpos=0; /* Position of error */
161 int kstat; /* Local status variable */
162 double twx,twy,twz; /* Local version of normalized direction vector */
163 double tx0,ty0,tz0; /* Local version of point on axis */
164 double temp; /* Temporary storage variable */
165 double tcost2; /* The square of the cosine of the opening angle */
166 double sdirec[3]; /* Normalized direction of cone axis */
167 double sdcone[3]; /* Normalized vector from top point to cone point*/
168
169
170 /* Test i legal input */
171 if (inumb <1 ) inumb = 1;
172 if (idim != 3 ) goto err104;
173
174 kdimp1 = idim + 1;
175 kdimp2 = idim + 2;
176 kstop = kdimp1*kdimp1;
177
178 for (ki=0;ki<kstop;ki++)
179 {
180 carray[ki] = DZERO;
181 }
182
183 /* Normalize direction vector of axis */
184
185 s6diff(etop,eaxis,idim,sdirec);
186 (void)s6norm(sdirec,idim,sdirec,&kstat);
187
188 /* Normalize vector from top point to point on conic surface */
189
190 s6diff(etop,econe,idim,sdcone);
191 (void)s6norm(sdcone,idim,sdcone,&kstat);
192
193 /* Make cosinus of angle between the two normalized vectors */
194
195 temp = s6scpr(sdirec,sdcone,idim);
196 tcost2 = temp*temp;
197
198 /* Test if cone degenerate */
199 if (DEQUAL(tcost2,DZERO)) goto err174;
200
201 /* Make diagonal elements */
202
203 for (ki=0,kl=0 ; kl<idim ; kl++,ki+=kdimp2)
204 {
205 temp = sdirec[kl];
206 carray[ki] = (double)1.0 - temp*temp/tcost2;
207 }
208
209 twx = sdirec[0];
210 twy = sdirec[1];
211 twz = sdirec[2];
212 tx0 = etop[0];
213 ty0 = etop[1];
214 tz0 = etop[2];
215
216
217 /* Make element (1,4) and (4,1) */
218
219 temp = (tx0*twx*twx + twx*(ty0*twy+tz0*twz))/tcost2 - tx0;
220
221 carray[3] = temp;
222 carray[12] = temp;
223
224
225 /* Make element (2,4) and (4,2) */
226
227 temp = (ty0*twy*twy + twy*(tz0*twz+tx0*twx))/tcost2 - ty0;
228
229 carray[7] = temp;
230 carray[13] = temp;
231
232
233 /* Make element (3,4) and (4,3) */
234
235 temp = (tz0*twz*twz + twz*(tx0*twx+ty0*twy))/tcost2 - tz0;
236
237 carray[11] = temp;
238 carray[14] = temp;
239
240 /* Make element (4,4) */
241
242 temp = tx0*tx0 + ty0*ty0 + tz0*tz0
243 - ( tx0*tx0*twx*twx + ty0*ty0*twy*twy + tz0*tz0*twz*twz
244 + (double)2.0*tx0*ty0*twx*twy + (double)2.0*ty0*tz0*twy*twz
245 + (double)2.0*tz0*tx0*twz*twx )/tcost2;
246 carray[15] = temp;
247
248
249 /* Make element (1,2) and (2,1) */
250
251 temp = -twx*twy/tcost2;
252 carray[1] = temp;
253 carray[4] = temp;
254
255
256 /* Make element (1,3) and (3,1) */
257
258 temp = -twx*twz/tcost2;
259 carray[2] = temp;
260 carray[8] = temp;
261
262
263 /* Make element (2,3) and (3,2) */
264
265 temp = -twy*twz/tcost2;
266 carray[6] = temp;
267 carray[9]= temp;
268
269
270 /* Make extra copies of cylinder */
271
272 kj = kstop;
273 for (ki=1;ki<inumb;ki++)
274 {
275 for (kl=0;kl<kstop;kl++,kj++)
276 {
277 carray[kj] = carray[kl];
278 }
279 }
280
281 *jstat = 0;
282 goto out;
283
284 /* Dimension less than 1 */
285 err104: *jstat = -104;
286 s6err("s1323",*jstat,kpos);
287 goto out;
288
289 /* Degenerate cond */
290 err174: *jstat = -174;
291 s6err("s1323",*jstat,kpos);
292 goto out;
293 out:
294 return;
295 }
296