1 struct GRAVMAG_XY {
2 double x, y;
3 };
4
gaussian(double rad,double half_width)5 double gaussian(double rad, double half_width) {
6 double y, gauss_ct = -0.5 / (half_width * half_width);
7 y = exp(rad * rad * gauss_ct);
8 return (y);
9 }
10
helper_fun(struct GMTGRAVMAG3D_CTRL * Ctrl,struct GRAVMAG_XY * ellipse[2],int m,int j1,int j2,int j3,int i1,int i2,int i3,double z1,double z2,double z3)11 void helper_fun(struct GMTGRAVMAG3D_CTRL *Ctrl, struct GRAVMAG_XY *ellipse[2], int m, int j1, int j2, int j3, int i1, int i2, int i3, double z1, double z2, double z3) {
12 /* j1, j2 & j2 can only be 0 ot 1 */
13 Ctrl->raw_mesh[m].t1[0] = ellipse[j1][i1].x;
14 Ctrl->raw_mesh[m].t1[1] = -ellipse[j1][i1].y;
15 Ctrl->raw_mesh[m].t1[2] = -z1; /* -1 because Z is positive down in Okabe */
16
17 Ctrl->raw_mesh[m].t2[0] = ellipse[j2][i2].x;
18 Ctrl->raw_mesh[m].t2[1] = -ellipse[j2][i2].y;
19 Ctrl->raw_mesh[m].t2[2] = -z2;
20
21 Ctrl->raw_mesh[m].t3[0] = ellipse[j3][i3].x;
22 Ctrl->raw_mesh[m].t3[1] = -ellipse[j3][i3].y;
23 Ctrl->raw_mesh[m].t3[2] = -z3;
24 }
25
prism(struct GMT_CTRL * GMT,struct GMTGRAVMAG3D_CTRL * Ctrl,int nb)26 int prism(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int nb) {
27 double a, b, c, z_c, x0, y0, z_top, z_bot;
28 int i_tri = Ctrl->n_raw_triang;
29
30 a = Ctrl->M.params[PRISM][nb][0]; b = Ctrl->M.params[PRISM][nb][1]; c = Ctrl->M.params[PRISM][nb][2];
31 z_c = Ctrl->M.params[PRISM][nb][3]; x0 = Ctrl->M.params[PRISM][nb][4]; y0 = Ctrl->M.params[PRISM][nb][5];
32
33 Ctrl->raw_mesh = gmt_M_memory (GMT, Ctrl->raw_mesh, Ctrl->n_raw_triang + 12, struct GMTGRAVMAG3D_RAW);
34 z_top = -(z_c + c); z_bot = -z_c ; /* -1 because Z is positive down in Okabe */
35
36 /* vertex of top rectangle */
37 /* first triangle */
38 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_top;
39 Ctrl->raw_mesh[i_tri].t2[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
40 Ctrl->raw_mesh[i_tri].t3[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_top;
41 /* second triangle */
42 i_tri++;
43 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_top;
44 Ctrl->raw_mesh[i_tri].t2[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
45 Ctrl->raw_mesh[i_tri].t3[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_top;
46 /* vertex of east rectangle */
47 /* first triangle */
48 i_tri++;
49 Ctrl->raw_mesh[i_tri].t1[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
50 Ctrl->raw_mesh[i_tri].t2[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
51 Ctrl->raw_mesh[i_tri].t3[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_top;
52 /* second triangle */
53 i_tri++;
54 Ctrl->raw_mesh[i_tri].t1[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
55 Ctrl->raw_mesh[i_tri].t2[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
56 Ctrl->raw_mesh[i_tri].t3[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_bot;
57 /* vertex of north rectangle */
58 /* first triangle */
59 i_tri++;
60 Ctrl->raw_mesh[i_tri].t1[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
61 Ctrl->raw_mesh[i_tri].t2[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
62 Ctrl->raw_mesh[i_tri].t3[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_top;
63 /* second triangle */
64 i_tri++;
65 Ctrl->raw_mesh[i_tri].t1[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
66 Ctrl->raw_mesh[i_tri].t2[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
67 Ctrl->raw_mesh[i_tri].t3[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_bot;
68 /* vertex of west rectangle */
69 /* first triangle */
70 i_tri++;
71 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
72 Ctrl->raw_mesh[i_tri].t2[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
73 Ctrl->raw_mesh[i_tri].t3[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_top;
74 /* second triangle */
75 i_tri++;
76 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
77 Ctrl->raw_mesh[i_tri].t2[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
78 Ctrl->raw_mesh[i_tri].t3[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_bot;
79 /* vertex of south rectangle */
80 /* first triangle */
81 i_tri++;
82 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
83 Ctrl->raw_mesh[i_tri].t2[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
84 Ctrl->raw_mesh[i_tri].t3[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_top;
85 /* second triangle */
86 i_tri++;
87 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
88 Ctrl->raw_mesh[i_tri].t2[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_top;
89 Ctrl->raw_mesh[i_tri].t3[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_bot;
90 /* vertex of bottom rectangle */
91 /* first triangle */
92 i_tri++;
93 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
94 Ctrl->raw_mesh[i_tri].t2[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_bot;
95 Ctrl->raw_mesh[i_tri].t3[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_bot;
96 /* second triangle */
97 i_tri++;
98 Ctrl->raw_mesh[i_tri].t1[0] = -a/2 + x0; Ctrl->raw_mesh[i_tri].t1[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t1[2] = z_bot;
99 Ctrl->raw_mesh[i_tri].t2[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t2[1] = b/2 - y0; Ctrl->raw_mesh[i_tri].t2[2] = z_bot;
100 Ctrl->raw_mesh[i_tri].t3[0] = a/2 + x0; Ctrl->raw_mesh[i_tri].t3[1] = -b/2 - y0; Ctrl->raw_mesh[i_tri].t3[2] = z_bot;
101
102 Ctrl->n_raw_triang += 12;
103 return 0;
104 }
105
five_psoid(struct GMT_CTRL * GMT,struct GMTGRAVMAG3D_CTRL * Ctrl,int body_type,int nb,bool cone,bool piram,bool sino,bool hemi)106 int five_psoid(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int body_type, int nb, bool cone, bool piram, bool sino, bool hemi) {
107 /* Constructs either a sphere, ellipsoid, cone, pyramid, or a bell */
108 /* as a union of triangular facets. Returns number of triangles. */
109 int i, j, j1, k, l, m, m1, m2, n = 0, n_tri, i_tri, npts_circ, n_slices, n_sigmas = 2;
110 bool first = true;
111 double a, b, c, z_c, x0, y0, z_top, z_bot;
112 double dfi, d_sli, ai0, ai1, bi0, bi1, zi0, zi1;
113 double d_tet, half_width_x, half_width_y, dx, dy, dz, rad_x, rad_y;
114 struct GRAVMAG_XY *ellipse[2];
115
116 i_tri = Ctrl->n_raw_triang; /* Start over any previous raw triang collection */
117
118 /* Fish the these parameters from Ctrl struct */
119 switch (body_type) {
120 case BELL:
121 n_sigmas = Ctrl->M.params[body_type][nb][6];
122 npts_circ = Ctrl->M.params[body_type][nb][7];
123 n_slices = Ctrl->M.params[body_type][nb][8];
124 break;
125 case ELLIPSOID:
126 npts_circ = Ctrl->M.params[body_type][nb][6];
127 n_slices = Ctrl->M.params[body_type][nb][7];
128 break;
129 case CONE:
130 npts_circ = Ctrl->M.params[body_type][nb][4];
131 n_slices = Ctrl->n_slices;
132 break;
133 case SPHERE:
134 npts_circ = Ctrl->M.params[body_type][nb][4];
135 n_slices = Ctrl->M.params[body_type][nb][5];
136 break;
137 default:
138 npts_circ = Ctrl->npts_circ;
139 n_slices = Ctrl->n_slices;
140 break;
141 }
142
143 a = Ctrl->M.params[body_type][nb][0]; b = Ctrl->M.params[body_type][nb][1];
144 c = Ctrl->M.params[body_type][nb][2]; z_c = Ctrl->M.params[body_type][nb][3];
145 x0 = Ctrl->M.params[body_type][nb][4]; y0 = Ctrl->M.params[body_type][nb][5];
146 z_top = z_c + c; z_bot = z_c;
147
148 n_tri = (hemi) ? 2 * npts_circ * n_slices : 2 * (npts_circ * (n_slices*2 - 1));
149 ellipse[0] = (struct GRAVMAG_XY *) calloc((size_t) (npts_circ+1), sizeof(struct GRAVMAG_XY));
150 ellipse[1] = (struct GRAVMAG_XY *) calloc((size_t) (npts_circ+1), sizeof(struct GRAVMAG_XY));
151
152 Ctrl->n_raw_triang += n_tri;
153 Ctrl->raw_mesh = gmt_M_memory(GMT, Ctrl->raw_mesh, Ctrl->n_raw_triang, struct GMTGRAVMAG3D_RAW);
154
155 dfi = (TWO_PI / npts_circ); d_tet = (M_PI_2 / n_slices);
156 d_sli = c / n_slices;
157 half_width_x = 0.5 * a; half_width_y = 0.5 * b;
158 dx = n_sigmas * a / n_slices; dy = n_sigmas * b / n_slices;
159 dz = c / n_slices; /* repeated but ok */
160
161 for (j = 0; j < Ctrl->n_slices; j++) {
162 j1 = j + 1;
163 if (cone || piram) {
164 ai0 = j * a / n_slices; bi0 = j * b / n_slices;
165 zi0 = z_top - j * d_sli;
166 ai1 = j1 * a / n_slices; bi1 = j1 * b / n_slices;
167 zi1 = z_top - j1 * d_sli;
168 }
169 else if (sino) { /* Bell shaped volume */
170 rad_x = j * dx; rad_y = j * dy;
171 ai0 = rad_x; bi0 = rad_y;
172 zi0 = z_bot + c * gaussian(rad_x, half_width_x);
173 rad_x = j1 * dx; rad_y = j1 * dy;
174 ai1 = rad_x; bi1 = rad_y;
175 zi1 = z_bot + c * gaussian(rad_x, half_width_x);
176 }
177 else { /* Ellipsoide or Sphere */
178 ai0 = a*cos(M_PI_2-j*d_tet); bi0 = b*cos(M_PI_2-j*d_tet);
179 zi0 = z_top - c * (1. - sqrt(1. - (ai0/a)*(ai0/a)));
180 ai1 = a*cos(M_PI_2-j1*d_tet); bi1 = b*cos(M_PI_2-j1*d_tet);
181 zi1 = z_top - c * (1. - sqrt(1. - (ai1/a)*(ai1/a)));
182 }
183 for (i = 0; i < npts_circ; i++) { /* compute slice j */
184 ellipse[0][i].x = x0 + ai0 * cos (i*dfi);
185 ellipse[0][i].y = y0 + bi0 * sin (i*dfi);
186 ellipse[1][i].x = x0 + ai1 * cos (i*dfi);
187 ellipse[1][i].y = y0 + bi1 * sin (i*dfi);
188 }
189 ellipse[0][npts_circ].x = ellipse[0][0].x; /* close slice "contour" */
190 ellipse[0][npts_circ].y = ellipse[0][0].y;
191 ellipse[1][npts_circ].x = ellipse[1][0].x;
192 ellipse[1][npts_circ].y = ellipse[1][0].y;
193
194 /* Calculates vertex of triangles in slice j */
195 i = 0;
196 if (first) {
197 for (m = 0; m < npts_circ ; m++) {
198 helper_fun(Ctrl, ellipse, m + i_tri, 0, 1, 1, i, i+1, i, zi0, zi1, zi1);
199 i++;
200 }
201 }
202 else {
203 for (m = (j-1)*npts_circ + i_tri; m < j*npts_circ + i_tri; m++) {
204 /* First triangle */
205 m1 = 2 * m + npts_circ; m2 = 2 * m + 1 + npts_circ;
206 helper_fun(Ctrl, ellipse, m1, 0, 1, 1, i, i+1, i, zi0, zi1, zi1);
207 /* Second triangle */
208 helper_fun(Ctrl, ellipse, m2, 0, 0, 1, i, i+1, i+1, zi0, zi0, zi1);
209 i++;
210 }
211 }
212 first = false;
213 }
214 /* First half is ready. Now, either close it and return or construct the other half by simetry */
215 if (cone || piram || sino || hemi) { /* close the base and return */
216 if (sino && fabs ((zi1 - z_c) / z_c) > 0.01) {
217 /* bell's last slice is 1% far from base, so we add a vertical wall */
218 z_top = zi1; /* update z_top value */
219 for (n = 0; n < npts_circ; n++) {
220 /* First triangle */
221 j = npts_circ * (n_slices * 2 - 1) + 2 * (n + i_tri);
222 helper_fun(Ctrl, ellipse, j, 1, 1, 1, n, n+1, n, z_top, z_top, z_bot);
223 /* Second triangle */
224 j1 = npts_circ * (n_slices * 2 - 1) + 2 * (n + i_tri) + 1;
225 helper_fun(Ctrl, ellipse, j1, 1, 1, 1, n+1, n+1, n, z_top, z_top, z_bot);
226 }
227 }
228 else if (sino) /* slightly change base depth to force a closed volume */
229 z_c = zi1;
230
231 i = npts_circ * (n_slices*2 - 1) + 2 * n;
232 for (k = i, l = 0; k < i+npts_circ; k++, l++) {
233 helper_fun(Ctrl, ellipse, k + i_tri, 1, 1, 1, 0, l, l+1, z_c, z_c, z_c);
234 Ctrl->raw_mesh[k + i_tri].t1[0] = x0; /* These two fell of the general case of the helper_fun() */
235 Ctrl->raw_mesh[k + i_tri].t1[1] = y0;
236 }
237 free((void *)ellipse[0]);
238 free((void *)ellipse[1]);
239 return (k);
240 }
241 n_tri = npts_circ * (n_slices*2 - 1);
242 for (j = n_tri-1, i = n_tri; j >= 0; j--, i++) {
243 Ctrl->raw_mesh[i+i_tri].t1[0] = Ctrl->raw_mesh[j+i_tri].t1[0];
244 Ctrl->raw_mesh[i+i_tri].t1[1] = Ctrl->raw_mesh[j+i_tri].t1[1];
245 Ctrl->raw_mesh[i+i_tri].t1[2] = z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t1[2]);
246
247 Ctrl->raw_mesh[i+i_tri].t2[0] = Ctrl->raw_mesh[j+i_tri].t2[0];
248 Ctrl->raw_mesh[i+i_tri].t2[1] = Ctrl->raw_mesh[j+i_tri].t2[1];
249 Ctrl->raw_mesh[i+i_tri].t2[2] = z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t2[2]);
250
251 Ctrl->raw_mesh[i+i_tri].t3[0] = Ctrl->raw_mesh[j+i_tri].t3[0];
252 Ctrl->raw_mesh[i+i_tri].t3[1] = Ctrl->raw_mesh[j+i_tri].t3[1];
253 Ctrl->raw_mesh[i+i_tri].t3[2] = z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t3[2]);
254 }
255 free ((void *)ellipse[0]);
256 free ((void *)ellipse[1]);
257 n_tri = 2 * (npts_circ * (n_slices*2 - 1));
258 return (n_tri);
259 }
260
cilindro(struct GMT_CTRL * GMT,struct GMTGRAVMAG3D_CTRL * Ctrl,int nb)261 int cilindro (struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int nb) {
262 double z_top, z_bot, dfi, rad_c, height_c, z_c, x0, y0;
263 int i, j, j1, n_tri, i_tri, npts_circ;
264 struct GRAVMAG_XY *circ;
265
266 i_tri = Ctrl->n_raw_triang; /* Start over any previous raw triang collection */
267
268 rad_c = Ctrl->M.params[CYLINDER][nb][0]; height_c = Ctrl->M.params[CYLINDER][nb][1];
269 z_c = Ctrl->M.params[CYLINDER][nb][2];
270 x0 = Ctrl->M.params[CYLINDER][nb][3]; y0 = Ctrl->M.params[CYLINDER][nb][4];
271 npts_circ = Ctrl->M.params[CYLINDER][nb][5];
272 z_top = -(z_c + height_c); z_bot = -z_c ;
273
274 n_tri = Ctrl->npts_circ * 4;
275 circ = (struct GRAVMAG_XY *) calloc((size_t) (Ctrl->npts_circ+1), sizeof(struct GRAVMAG_XY));
276
277 Ctrl->n_raw_triang += n_tri;
278 Ctrl->raw_mesh = gmt_M_memory(GMT, Ctrl->raw_mesh, Ctrl->n_raw_triang, struct GMTGRAVMAG3D_RAW);
279
280 dfi = (TWO_PI / npts_circ);
281
282 for (i = 0; i < npts_circ; i++) { /* compute circle */
283 circ[i].x = x0 + rad_c * cos (i*dfi);
284 circ[i].y = y0 + rad_c * sin (i*dfi);
285 }
286 circ[npts_circ].x = circ[0].x; circ[npts_circ].y = circ[0].y;
287
288 for (i = 0; i < Ctrl->npts_circ; i++) { /* Calculates vertex of top circle */
289 Ctrl->raw_mesh[i+i_tri].t1[0] = x0;
290 Ctrl->raw_mesh[i+i_tri].t1[1] = -y0;
291 Ctrl->raw_mesh[i+i_tri].t1[2] = z_top;
292
293 Ctrl->raw_mesh[i+i_tri].t2[0] = circ[i+1].x;
294 Ctrl->raw_mesh[i+i_tri].t2[1] = -circ[i+1].y;
295 Ctrl->raw_mesh[i+i_tri].t2[2] = z_top;
296
297 Ctrl->raw_mesh[i+i_tri].t3[0] = circ[i].x;
298 Ctrl->raw_mesh[i+i_tri].t3[1] = -circ[i].y;
299 Ctrl->raw_mesh[i+i_tri].t3[2] = z_top;
300 }
301 for (i = 0; i < npts_circ; i++) { /* Computes vertex of side rectangle, where each one is decomposed into two triangles */
302 /* First triangle */
303 j = 2 * (i + i_tri) + npts_circ;
304 Ctrl->raw_mesh[j].t1[0] = circ[i].x; Ctrl->raw_mesh[j].t1[1] = -circ[i].y; Ctrl->raw_mesh[j].t1[2] = z_top;
305 Ctrl->raw_mesh[j].t2[0] = circ[i+1].x; Ctrl->raw_mesh[j].t2[1] = -circ[i+1].y; Ctrl->raw_mesh[j].t2[2] = z_top;
306 Ctrl->raw_mesh[j].t3[0] = circ[i].x; Ctrl->raw_mesh[j].t3[1] = -circ[i].y; Ctrl->raw_mesh[j].t3[2] = z_bot;
307 /* Second triangle */
308 j1 = 2 * (i + i_tri) + npts_circ + 1;
309 Ctrl->raw_mesh[j1].t1[0] = circ[i+1].x; Ctrl->raw_mesh[j1].t1[1] = -circ[i+1].y; Ctrl->raw_mesh[j1].t1[2] = z_top;
310 Ctrl->raw_mesh[j1].t2[0] = circ[i+1].x; Ctrl->raw_mesh[j1].t2[1] = -circ[i+1].y; Ctrl->raw_mesh[j1].t2[2] = z_bot;
311 Ctrl->raw_mesh[j1].t3[0] = circ[i].x; Ctrl->raw_mesh[j1].t3[1] = -circ[i].y; Ctrl->raw_mesh[j1].t3[2] = z_bot;
312 }
313 for (i = 0; i < npts_circ; i++) { /* Now the vertex of bottom circle */
314 j = (i + i_tri) + 3 * npts_circ;
315 Ctrl->raw_mesh[j].t1[0] = x0; Ctrl->raw_mesh[j].t1[1] = -y0; Ctrl->raw_mesh[j].t1[2] = z_bot;
316 Ctrl->raw_mesh[j].t2[0] = circ[i].x; Ctrl->raw_mesh[j].t2[1] = -circ[i].y; Ctrl->raw_mesh[j].t2[2] = z_bot;
317 Ctrl->raw_mesh[j].t3[0] = circ[i+1].x; Ctrl->raw_mesh[j].t3[1] = -circ[i+1].y; Ctrl->raw_mesh[j].t3[2] = z_bot;
318 }
319 free ((void *)circ);
320 return n_tri;
321 }
322