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: s1023.c,v 1.3 2001-03-19 15:58:41 afr Exp $
45 *
46 */
47
48
49 #define S1023
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1023(double center[],double axis[],double equator[],int latitude,int longitude,SISLSurf ** sphere,int * stat)55 s1023(double center[], double axis[], double equator[], int latitude,
56 int longitude, SISLSurf **sphere, int *stat)
57 #else
58 void s1023(center, axis, equator, latitude, longitude, sphere, stat)
59 double center[];
60 double axis[];
61 double equator[];
62 int latitude;
63 int longitude;
64 SISLSurf **sphere;
65 int *stat;
66 #endif
67 /*
68 *********************************************************************
69 *
70 * PURPOSE : To describe octants of a sphere as a NURBS. This can also
71 * be used to describe the complete sphere.
72 *
73 *
74 * INPUT : center - Center point of the sphere
75 * axis - Axis of the sphere (towards the north pole)
76 * equator - Vector from center to start point on the equator
77 * latitude - Flag indicating number of octants in north/south
78 * direction:
79 * = 1 : Octants in the northern hemisphere
80 * = 2 : Octants in both hemispheres
81 * longitude - Flag indicating number of octants along the
82 * equator:
83 * = 1 : Octants in 1. quadrant
84 * = 2 : Octants in 1. and 2. quadrant
85 * = 3 : Octants in 1., 2. and 3. quadrant
86 * = 4 : Octants in all quadrants
87 * This is counted counter-clockwise from equator
88 *
89 *
90 * OUTPUT :
91 * stat - status messages
92 * > 0 : warning
93 * = 0 : ok
94 * < 0 : error
95 * spher - Pointer to the sphere produced
96 *
97 * METHOD :
98 *
99 *
100 * REFERENCES :
101 *
102 *-
103 * CALLS :
104 *
105 * WRITTEN BY : Johannes Kaasa, SI, Oslo, Norway, Jan. 93
106 * Revised by : Paal Fugelli and Johannes Kaasa, SINTEF, Oslo, Norway, 08-94.
107 *
108 *********************************************************************
109 */
110 {
111 int kstat; /* Status variable. */
112 int kpos=0; /* Position of error. */
113 int ki, kj, kl; /* Indexes in for loops. */
114 int in1; /* Number of vertices along the longitude. */
115 int in2; /* Number of vertices along the latitude. */
116 int ik1 = 3; /* Order along the longitude. */
117 int ik2 = 3; /* Order along the latitude. */
118 double *et1 = SISL_NULL; /* Knot vector along the longitude. */
119 double *et2 = SISL_NULL; /* Knot vector along the latitude. */
120 double *rcoef = SISL_NULL; /* Coefficients of the sphere. */
121 int kind = 2; /* Rational Bspline surface. */
122 double weight; /* Rational weight. */
123 double radius; /* Radius of the sphere. */
124 double norm; /* Length of vectors. */
125 double x_axis[3]; /* axis normalized with radius length. */
126 double z_axis[3]; /* Radius vector normal to x_axis and equator. */
127 double w1, w2; /* Rational weights in both directions. */
128 double x_comp; /* Component in local x direction. */
129 double y_comp; /* Component in local y direction. */
130 double z_comp; /* Component in local z direction. */
131
132 /* Do necessary initiation and allocation. */
133
134 *sphere = SISL_NULL;
135 weight = (double)1.0/sqrt(2.0);
136 in1 = 1 + 2*latitude;
137 in2 = 1 + 2*longitude;
138
139 radius = s6length(equator, 3, &kstat);
140 if (kstat < 0) goto error;
141 norm = s6length(axis, 3, &kstat);
142 if (kstat < 0) goto error;
143 for (ki = 0; ki < 3; ki++)
144 x_axis[ki] = radius*axis[ki]/norm;
145 s6crss(x_axis, equator, z_axis);
146 norm = s6length(z_axis, 3, &kstat);
147 if (kstat < 0) goto error;
148 for (ki = 0; ki < 3; ki++)
149 z_axis[ki] = radius*z_axis[ki]/norm;
150
151 if((et1 = newarray(in1 + ik1, DOUBLE)) == SISL_NULL) goto err101;
152 if((et2 = newarray(in2 + ik2, DOUBLE)) == SISL_NULL) goto err101;
153 if((rcoef = newarray(4*in1*in2, DOUBLE)) == SISL_NULL) goto err101;
154
155 /* Initiate the knot vectors. */
156
157 for (ki = 0; ki < ik1; ki++)
158 et1[ki] = (double)0.;
159 for (ki = 0; ki < latitude; ki++)
160 {
161 et1[ik1 + 2*ki] = (ki + 1)*PIHALF;
162 et1[ik1 + 2*ki + 1] = (ki + 1)*PIHALF;
163 }
164 et1[in1 + ik1 - 1] = latitude*PIHALF;
165
166 for (ki = 0; ki < ik2; ki++)
167 et2[ki] = (double)0.;
168 for (ki = 0; ki < longitude; ki++)
169 {
170 et2[ik2 + 2*ki] = (ki + 1)*PIHALF;
171 et2[ik2 + 2*ki + 1] = (ki + 1)*PIHALF;
172 }
173 et2[in2 + ik2 - 1] = longitude*PIHALF;
174
175 /* Initiate the coefficient vector. */
176
177 for (ki = 0; ki < in2; ki++)
178 {
179 if (ki == 1 || ki == 3 || ki == 5 || ki == 7)
180 w2 = weight;
181 else
182 w2 = (double)1.;
183
184 if (ki == 0 || ki == 1 || ki == 7 || ki == 8)
185 y_comp = (double)1.;
186 else if (ki == 3 || ki == 4 || ki == 5)
187 y_comp = - (double)1.;
188 else
189 y_comp = (double)0.;
190
191 if (ki == 1 || ki == 2 || ki == 3)
192 z_comp = (double)1.;
193 else if (ki == 5 || ki == 6 || ki == 7)
194 z_comp = - (double)1.;
195 else
196 z_comp = (double)0.;
197
198 for (kj = 0; kj < in1; kj++)
199 {
200 if (kj == 1 || kj == 3)
201 w1 = weight;
202 else
203 w1 = (double)1.;
204
205 if (kj == 0 || kj == 1)
206 x_comp = (double)1.;
207 else if (kj == 3 || kj == 4)
208 x_comp = - (double)1.;
209 else
210 x_comp = (double)0.;
211
212 w1 *= w2;
213
214 if (kj == 0 || kj == 4)
215 {
216 for (kl = 0; kl < 3; kl++)
217 rcoef[4*(ki*in1 + kj) + kl] = w1*(center[kl] + x_comp*x_axis[kl]);
218 }
219 else
220 {
221 for (kl = 0; kl < 3; kl++)
222 rcoef[4*(ki*in1 + kj) + kl] = w1*(center[kl] + x_comp*x_axis[kl]
223 + y_comp*equator[kl] + z_comp*z_axis[kl]);
224 }
225 rcoef[4*(ki*in1 + kj) + 3] = w1;
226 }
227 }
228 (*sphere) = newSurf(in1, in2, ik1, ik2, et1, et2, rcoef, kind, 3, 1);
229 if ((*sphere) == SISL_NULL) goto err101;
230
231 if (et1 != SISL_NULL) freearray(et1);
232 if (et2 != SISL_NULL) freearray(et2);
233 if (rcoef != SISL_NULL) freearray(rcoef);
234
235 *stat = 0;
236 goto out;
237
238 /* Error in curve allocation. */
239
240 err101:
241 *stat = -101;
242 s6err("s1023",*stat,kpos);
243 goto out;
244
245 /* Error in lower level routine. */
246
247 error:
248 *stat = kstat;
249 s6err("s1023", *stat, kpos);
250 goto out;
251
252 out:
253 return;
254 }
255