1 // $Id: drawmap.cxx 1064 2010-10-22 19:38:20Z martin $
2 //
3 // drawmap.cpp - routine for DRAWxtl V5.5 - the GUI version
4 //
5 // Coded using the FLTK 1.1.6 widget set
6 //
7 // Larry W. Finger, Martin Kroeker and Brian Toby
8 //
9 // This module includes Fourier support routines
10 //
11 // routines contained within this file:
12 //
13 // generate_map - Generate Marching Cubes triangles from a map
14 // InterpolateMap - interpolate map value
15 // LookupMap - return the value of the mmap at one of the calculated points
16 // Maximize_rho - position the cursor at the local maximum (Q & D routine)
17
18 #include "drawxtl.h"
19 #include <stdio.h>
20 #include <string.h>
21 #include <math.h>
22 #include <stdlib.h>
23 #include <time.h>
24 #include "draw_ext.h"
25 #include "drawmap.h"
26 #include "DRAWxtlViewUI.h"
27 #include "MC.h"
28
29 #include "DRAWxtl_proto.h"
30
31 /* *************************************************************************************** */
32 /* Generate Marching Cubes triangles from a map with contour minValue*/
33 /* *************************************************************************************** */
34
35 void
generate_map(float minValue,int Solid,char * Color,char * BackColor)36 generate_map (float minValue, int Solid, char *Color, char *BackColor)
37 {
38 int numOfTriangles;
39
40 float step[3];
41
42 float Red, Green, Blue;
43
44 int nxMin, nxMax, nyMin, nyMax, nzMin, nzMax; /* max and min steps along a, b & c */
45
46 int snx, sny, snz; /* number of steps along a, b & c */
47
48 unsigned int sny_snz; /* sny x snz */
49
50 int i, j, k, ijk, ijkmax, ii, jj;
51
52 unsigned int ni, nj, nk, si, sj, sk, sni, snj, sind;
53
54 mp4Vector *PointsList;
55
56 float fvert[3], cvert[3];
57
58 TRIANGLE *pTriangles;
59
60 char string[40];
61
62 char Color2[40];
63
64 float d1, d2, d3;
65
66 int iproj = 0;
67
68 int iMax;
69
70 int jMax;
71
72 int kMax;
73
74 if (doVrml) {
75 if (!drvui->fpoutv) {
76 Error_Box ("Invalid vrml file path in drawmap.");
77 return;
78 }
79 }
80 if (!FourierPt)
81 return;
82
83 drvui->mainWindow->cursor (FL_CURSOR_WAIT);
84 strcpy (string, Color);
85 strcpy (Color2, Color);
86 Transform_VRML_Color (string);
87 Transform_POV_Color (Color2);
88 (void) sscanf (string, "%f %f %f", &Red, &Green, &Blue);
89
90 step[0] = 1.0f / (float) mapstep_a;
91 step[1] = 1.0f / (float) mapstep_b;
92 step[2] = 1.0f / (float) mapstep_c;
93 nxMin = (int) (max (drvui->frames[drvui->frame_no].map_lim[0], drvui->frames[drvui->frame_no].cryst_lim[0]) * mapstep_a);
94 nxMax = (int) (min (drvui->frames[drvui->frame_no].map_lim[3], drvui->frames[drvui->frame_no].cryst_lim[3]) * mapstep_a);
95 nyMin = (int) (max (drvui->frames[drvui->frame_no].map_lim[1], drvui->frames[drvui->frame_no].cryst_lim[1]) * mapstep_b);
96 nyMax = (int) (min (drvui->frames[drvui->frame_no].map_lim[4], drvui->frames[drvui->frame_no].cryst_lim[4]) * mapstep_b);
97 nzMin = (int) (max (drvui->frames[drvui->frame_no].map_lim[2], drvui->frames[drvui->frame_no].cryst_lim[2]) * mapstep_c);
98 nzMax = (int) (min (drvui->frames[drvui->frame_no].map_lim[5], drvui->frames[drvui->frame_no].cryst_lim[5]) * mapstep_c);
99 snx = abs (nxMax - nxMin + 1);
100 sny = abs (nyMax - nyMin + 1);
101 snz = abs (nzMax - nzMin + 1);
102 iMax = nxMax + 1;
103 jMax = nyMax + 1;
104 kMax = nzMax + 1;
105 if (snx == 1) {
106 iproj = 1;
107 snx++;
108 iMax++;
109 }
110 if (sny == 1) {
111 iproj = 2;
112 sny++;
113 jMax++;
114 }
115 if (snz == 1) {
116 iproj = 3;
117 snz++;
118 kMax++;
119 }
120 sny_snz = sny * snz;
121
122 ijkmax= (int) (drvui->frames[drvui->frame_no].cryst_lim[3] * (mapstep_a-1)) * mapstep_b * mapstep_c
123 + (int) (drvui->frames[drvui->frame_no].cryst_lim[4] * (mapstep_b-1)) * mapstep_c
124 + (int) (drvui->frames[drvui->frame_no].cryst_lim[5] * (mapstep_c-1));
125
126 /* create a new set of Fourier points */
127
128 PointsList =
129 (mp4Vector *) zalloc (sizeof (mp4Vector) * (snx + 1) * (sny + 1) * (snz + 1));
130 if (PointsList == NULL) {
131 Error_Box ("Error allocating space for PointsList vectors.");
132 return;
133 }
134 for (i = nxMin, si = 0; i <= iMax; i++, si++) {
135 ni = i;
136 ni = (5 * mapstep_a + ni) % mapstep_a; /* this will 'wrap' around any value (negative or positive) */
137 sni = si * sny_snz;
138 fvert[0] = i * step[0];
139 if (iproj == 1) {
140 ni = nxMin;
141 fvert[0] = ni * step[0] + 0.001f;
142 }
143 for (j = nyMin, sj = 0; j <= jMax; j++, sj++) {
144 nj = j;
145 nj = (5 * mapstep_b + nj) % mapstep_b;
146 snj = sj * snz;
147 fvert[1] = j * step[1];
148 if (iproj == 2) {
149 nj = nyMin;
150 fvert[1] = nj * step[1] + 0.001f;
151 }
152 for (k = nzMin, sk = 0; k <= kMax; k++, sk++) {
153 nk = k;
154 nk = (5 * mapstep_c + nk) % mapstep_c;
155 sind = sni + snj + sk;
156 fvert[2] = k * step[2];
157 if (iproj == 3) {
158 nk = nzMin;
159 fvert[2] = nk * step[2] + 0.001f;
160 }
161
162 /* convert Fractional to Orthonormal Coords */
163
164 for (ii = 0; ii <= 2; ++ii) { /* convert vertex coordinates to Cartesian */
165 cvert[ii] = 0.0f;
166 for (jj = 0; jj <= 2; ++jj)
167 cvert[ii] +=
168 (float) drvui->b_mat[ii][jj] * (fvert[jj] - origin[jj]);
169 }
170 PointsList[sind].x = cvert[0];
171 PointsList[sind].y = cvert[1];
172 PointsList[sind].z = cvert[2];
173 ijk = ni * mapstep_b * mapstep_c + nj * mapstep_c + nk;
174 if (ijk > ijkmax)
175 ijk-=ijkmax;
176 PointsList[sind].val = FourierPt[ijk];
177 }
178 }
179 }
180
181 /* create contours from grid of points with matching cubes */
182
183 pTriangles =
184 MC (snx - 1, sny - 1, snz - 1, step[0], step[1], step[2], minValue, PointsList,
185 numOfTriangles);
186 free (PointsList);
187
188 if (numOfTriangles <= 0)
189 return;
190
191 glPushMatrix ();
192 if (!drvui->frames[drvui->frame_no].slice) {
193 if (!Solid) {
194 glPolygonMode (GL_FRONT_AND_BACK, GL_LINE);
195 glDisable (GL_LIGHTING);
196 } else
197 glMaterialf (GL_FRONT, GL_SHININESS, 0.0); /* disable shinyness */
198 glColor3f (Red, Green, Blue);
199
200 glBegin (GL_TRIANGLES);
201
202 for (i = 0; i < numOfTriangles; i++) {
203 d1 = (pTriangles[i].p[0].x - pTriangles[i].p[1].x) *
204 (pTriangles[i].p[0].x - pTriangles[i].p[1].x) +
205 (pTriangles[i].p[0].y - pTriangles[i].p[1].y) *
206 (pTriangles[i].p[0].y - pTriangles[i].p[1].y) +
207 (pTriangles[i].p[0].z - pTriangles[i].p[1].z) *
208 (pTriangles[i].p[0].z - pTriangles[i].p[1].z);
209 d2 = (pTriangles[i].p[2].x - pTriangles[i].p[1].x) *
210 (pTriangles[i].p[2].x - pTriangles[i].p[1].x) +
211 (pTriangles[i].p[2].y - pTriangles[i].p[1].y) *
212 (pTriangles[i].p[2].y - pTriangles[i].p[1].y) +
213 (pTriangles[i].p[2].z - pTriangles[i].p[1].z) *
214 (pTriangles[i].p[2].z - pTriangles[i].p[1].z);
215 d3 = (pTriangles[i].p[2].x - pTriangles[i].p[0].x) *
216 (pTriangles[i].p[2].x - pTriangles[i].p[0].x) +
217 (pTriangles[i].p[2].y - pTriangles[i].p[0].y) *
218 (pTriangles[i].p[2].y - pTriangles[i].p[0].y) +
219 (pTriangles[i].p[2].z - pTriangles[i].p[0].z) *
220 (pTriangles[i].p[2].z - pTriangles[i].p[0].z);
221
222 // skip excursions completely across the cell when a contour reaches the edge
223 if ((d1 > 5.*snx*snx* step[0] ) || (d2 > 5.*sny*sny* step[1]) ||
224 (d3 > 5.*snz*snz* step[2]))
225 continue;
226
227 if (minValue > 0) {
228 if (doVrml) {
229 if (!Vrml2) {
230 fprintf (drvui->fpoutv, " Separator {\n");
231 fprintf (drvui->fpoutv, " Material { diffuseColor %s }\n", string);
232 fprintf (drvui->fpoutv, "Coordinate3 { point [\n");
233 } else {
234 fprintf (drvui->fpoutv, " Shape {\n");
235 if (Solid) {
236 fprintf (drvui->fpoutv, " appearance Appearance {\n");
237 fprintf (drvui->fpoutv,
238 " material Material {diffuseColor %s}\n", string);
239 fprintf (drvui->fpoutv,
240 " }\n geometry IndexedFaceSet { coord Coordinate{ point [\n");
241 } else
242 fprintf (drvui->fpoutv,
243 " geometry IndexedLineSet { coord Coordinate{ point [\n");
244 }
245 }
246 if (doPOV) {
247 if (!Solid) {
248 if (d1 > 0.00001f) {
249 fprintf (drvui->fpoutp, " cylinder { <%f,%f,%f>,<%f,%f,%f>,%f\n",
250 pTriangles[i].p[0].x, pTriangles[i].p[0].y,
251 pTriangles[i].p[0].z, pTriangles[i].p[1].x,
252 pTriangles[i].p[1].y, pTriangles[i].p[1].z,
253 0.0001 * Scale);
254 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n }\n",
255 Color2);
256 }
257 if (d2 > 0.00001f) {
258 fprintf (drvui->fpoutp, " cylinder { <%f,%f,%f>,<%f,%f,%f>,%f\n",
259 pTriangles[i].p[1].x, pTriangles[i].p[1].y,
260 pTriangles[i].p[1].z, pTriangles[i].p[2].x,
261 pTriangles[i].p[2].y, pTriangles[i].p[2].z,
262 0.0001 * Scale);
263 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n }\n",
264 Color2);
265 }
266 if (d3 > 0.00001f) {
267 fprintf (drvui->fpoutp, " cylinder { <%f,%f,%f>,<%f,%f,%f>,%f\n",
268 pTriangles[i].p[2].x, pTriangles[i].p[2].y,
269 pTriangles[i].p[2].z, pTriangles[i].p[0].x,
270 pTriangles[i].p[0].y, pTriangles[i].p[0].z,
271 0.0001 * Scale);
272 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n }\n",
273 Color2);
274 }
275 }
276 if (Solid)
277 fprintf (drvui->fpoutp, "smooth_triangle {\n");
278 }
279 for (j = 0; j < 3; j++) {
280 if (Solid)
281 glNormal3f (pTriangles[i].norm[j].x, pTriangles[i].norm[j].y,
282 pTriangles[i].norm[j].z);
283 glVertex3f (pTriangles[i].p[j].x, pTriangles[i].p[j].y,
284 pTriangles[i].p[j].z);
285 if (Solid && doPOV) {
286 fprintf (drvui->fpoutp, "<%8.5f, %8.5f, %8.5f>,", pTriangles[i].p[j].x,
287 pTriangles[i].p[j].y, pTriangles[i].p[j].z);
288 fprintf (drvui->fpoutp, "<%8.5f, %8.5f, %8.5f>", pTriangles[i].norm[j].x,
289 pTriangles[i].norm[j].y, pTriangles[i].norm[j].z);
290 }
291 if (doVrml)
292 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f", pTriangles[i].p[j].x,
293 pTriangles[i].p[j].y, pTriangles[i].p[j].z);
294
295 if (j < 2) {
296 if (Solid && doPOV)
297 fprintf (drvui->fpoutp, ",");
298 if (doVrml)
299 fprintf (drvui->fpoutv, ",\n");
300 } else {
301 if (Solid && doPOV) {
302 fprintf (drvui->fpoutp, "\n texture{pigment{color %s }}",
303 Color2);
304 if (BackColor)
305 fprintf (drvui->fpoutp, "\n interior_texture{pigment{color %s }}\n }\n",
306 BackColor);
307 else
308 fprintf(drvui->fpoutp, "\n }\n");
309 }
310 if (doVrml)
311 fprintf (drvui->fpoutv, "]\n");
312 }
313 } // for (j...
314 if (doAsy) {
315 if (!Solid)
316 fprintf (drvui->fpouta, "draw(pic, (%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--cycle,rgb(%.2f,%.2f,%.2f)+linewidth(%.2f));\n",
317 pTriangles[i].p[0].x, pTriangles[i].p[0].y,pTriangles[i].p[0].z,
318 pTriangles[i].p[1].x, pTriangles[i].p[1].y,pTriangles[i].p[1].z,
319 pTriangles[i].p[2].x, pTriangles[i].p[2].y,pTriangles[i].p[2].z,
320 Red,Green,Blue,0.001*Scale);
321 else
322 fprintf (drvui->fpouta, "draw(pic, surface( (%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--cycle),rgb(%.2f,%.2f,%.2f));\n",
323 pTriangles[i].p[0].x, pTriangles[i].p[0].y,pTriangles[i].p[0].z,
324 pTriangles[i].p[1].x, pTriangles[i].p[1].y,pTriangles[i].p[1].z,
325 pTriangles[i].p[2].x, pTriangles[i].p[2].y,pTriangles[i].p[2].z,
326 Red,Green,Blue);
327 }
328 if (doVrml) {
329 if (Vrml2) {
330 if (Solid) {
331 fprintf (drvui->fpoutv,
332 " }\n coordIndex [ 0,1,2,-1]\n solid FALSE\n convex"
333 " TRUE\n creaseAngle 1.57075\n }\n }\n");
334 } else {
335 fprintf (drvui->fpoutv,
336 " }\n coordIndex [ 0,1,2,-1]\n color Color { color [%s]"
337 "}\n colorIndex [0]\n colorPerVertex FALSE\n }\n }\n",
338 string);
339 }
340 } else {
341 if (Solid) {
342 fprintf (drvui->fpoutv,
343 "}\n IndexedFaceSet { coordIndex [0,1,2,-1] }\n}\n");
344 } else {
345 fprintf (drvui->fpoutv,
346 "}\n IndexedLineSet { coordIndex [0,1,2,-1] }\n}\n");
347 }
348 }
349 }
350 } else { // if (minValue > 0)
351 if (doVrml) {
352 if (!Vrml2) {
353 fprintf (drvui->fpoutv, " Separator {\n");
354 fprintf (drvui->fpoutv, " Material { diffuseColor %s }\n", string);
355 fprintf (drvui->fpoutv, "Coordinate3 { point [\n");
356 } else {
357 fprintf (drvui->fpoutv, " Shape {\n");
358 if (Solid) {
359 fprintf (drvui->fpoutv, " appearance Appearance {\n");
360 fprintf (drvui->fpoutv,
361 " material Material {diffuseColor %s}\n", string);
362 fprintf (drvui->fpoutv,
363 " }\n geometry IndexedFaceSet { coord Coordinate{ point [\n");
364 } else
365 fprintf (drvui->fpoutv,
366 " geometry IndexedLineSet { coord Coordinate{ point [\n");
367 }
368 }
369 if (doPOV) {
370 if (!Solid) {
371 if (d1 > 0.00001f) { // skip zero length cylinder
372 fprintf (drvui->fpoutp, " cylinder { <%f,%f,%f>,<%f,%f,%f>,%f\n",
373 pTriangles[i].p[0].x, pTriangles[i].p[0].y,
374 pTriangles[i].p[0].z, pTriangles[i].p[1].x,
375 pTriangles[i].p[1].y, pTriangles[i].p[1].z,
376 0.0001 * Scale);
377 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n }\n",
378 Color2);
379 }
380 if (d2 > 0.00001f) { // skip zero length cylinder
381 fprintf (drvui->fpoutp, " cylinder { <%f,%f,%f>,<%f,%f,%f>,%f\n",
382 pTriangles[i].p[1].x, pTriangles[i].p[1].y,
383 pTriangles[i].p[1].z, pTriangles[i].p[2].x,
384 pTriangles[i].p[2].y, pTriangles[i].p[2].z,
385 0.0001 * Scale);
386 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n }\n",
387 Color2);
388 }
389 if (d3 > 0.00001f) { // skip zero length cylinder
390 fprintf (drvui->fpoutp, " cylinder { <%f,%f,%f>,<%f,%f,%f>,%f\n",
391 pTriangles[i].p[2].x, pTriangles[i].p[2].y,
392 pTriangles[i].p[2].z, pTriangles[i].p[0].x,
393 pTriangles[i].p[0].y, pTriangles[i].p[0].z,
394 0.0001 * Scale);
395 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n }\n",
396 Color2);
397 }
398 }
399 if (Solid)
400 fprintf (drvui->fpoutp, "smooth_triangle {\n");
401 }
402 for (j = 2; j >= 0; j--) {
403 if (Solid)
404 glNormal3f (pTriangles[i].norm[j].x, pTriangles[i].norm[j].y,
405 pTriangles[i].norm[j].z);
406 glVertex3f (pTriangles[i].p[j].x, pTriangles[i].p[j].y,
407 pTriangles[i].p[j].z);
408
409 if (Solid && doPOV) {
410 fprintf (drvui->fpoutp, "<%8.5f, %8.5f, %8.5f>,", pTriangles[i].p[j].x,
411 pTriangles[i].p[j].y, pTriangles[i].p[j].z);
412 fprintf (drvui->fpoutp, "<%8.5f, %8.5f, %8.5f>", pTriangles[i].norm[j].x,
413 pTriangles[i].norm[j].y, pTriangles[i].norm[j].z);
414 }
415 if (doVrml)
416 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f", pTriangles[i].p[j].x,
417 pTriangles[i].p[j].y, pTriangles[i].p[j].z);
418 if (j > 0) {
419 if (Solid && doPOV)
420 fprintf (drvui->fpoutp, ",");
421 if (doVrml)
422 fprintf (drvui->fpoutv, ",\n");
423 } else {
424 if (Solid && doPOV) {
425 fprintf (drvui->fpoutp, "\n texture{pigment{color %s }}",
426 Color2);
427 if (BackColor)
428 fprintf (drvui->fpoutp, "\n interior_texture{pigment{color %s }}\n }\n",
429 BackColor);
430 else
431 fprintf(drvui->fpoutp, "\n }\n");
432 }
433 if (doVrml)
434 fprintf (drvui->fpoutv, "]\n");
435 }
436 }
437 if (doAsy) {
438 if (!Solid)
439 fprintf (drvui->fpouta, "draw(pic, (%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--cycle,rgb(%.2f,%.2f,%.2f)+linewidth(%.2f));\n",
440 pTriangles[i].p[0].x, pTriangles[i].p[0].y,pTriangles[i].p[0].z,
441 pTriangles[i].p[1].x, pTriangles[i].p[1].y,pTriangles[i].p[1].z,
442 pTriangles[i].p[2].x, pTriangles[i].p[2].y,pTriangles[i].p[2].z,
443 Red,Green,Blue,0.001*Scale);
444 else
445 fprintf (drvui->fpouta, "draw(pic, surface( (%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f)--cycle),rgb(%.2f,%.2f,%.2f));\n",
446 pTriangles[i].p[0].x, pTriangles[i].p[0].y,pTriangles[i].p[0].z,
447 pTriangles[i].p[1].x, pTriangles[i].p[1].y,pTriangles[i].p[1].z,
448 pTriangles[i].p[2].x, pTriangles[i].p[2].y,pTriangles[i].p[2].z,
449 Red,Green,Blue);
450 }
451 if (doVrml) {
452 if (Vrml2) {
453 if (Solid) {
454 fprintf (drvui->fpoutv,
455 "}\n coordIndex [ 1,2,3,-1]\n solid FALSE\n convex TRUE\n"
456 " creaseAngle 1.57075\n }\n }\n");
457 } else {
458 fprintf (drvui->fpoutv,
459 " }\n coordIndex [ 0,1,2,-1]\n color Color { color [%s]"
460 "}\n colorIndex [0]\n colorPerVertex FALSE\n }\n }\n",
461 string);
462 }
463 } else {
464 if (Solid) {
465 fprintf (drvui->fpoutv,
466 "}\n IndexedFaceSet { coordIndex [1,2,3,-1] }\n}\n");
467 } else {
468 fprintf (drvui->fpoutv,
469 "}\n IndexedLineSet { coordIndex [1,2,3,-1] }\n}\n");
470 }
471 }
472 }
473 } // if (minValue > 0)
474 }
475 glEnd ();
476 if (!Solid) {
477 glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);
478 glEnable (GL_LIGHTING);
479 }
480
481 } else { // 2D slices
482 mpVector p0,p1;
483
484 glPolygonMode (GL_FRONT_AND_BACK, GL_LINE);
485 glDisable (GL_LIGHTING);
486 glMaterialf (GL_FRONT, GL_SHININESS, 0.0); /* disable shinyness */
487 glColor3f (Red, Green, Blue);
488 glBegin (GL_LINES);
489 for (i = 0; i < numOfTriangles; i++) {
490 d1 = (pTriangles[i].p[0].x - pTriangles[i].p[1].x) *
491 (pTriangles[i].p[0].x - pTriangles[i].p[1].x) +
492 (pTriangles[i].p[0].y - pTriangles[i].p[1].y) *
493 (pTriangles[i].p[0].y - pTriangles[i].p[1].y) +
494 (pTriangles[i].p[0].z - pTriangles[i].p[1].z) *
495 (pTriangles[i].p[0].z - pTriangles[i].p[1].z);
496 d2 = (pTriangles[i].p[2].x - pTriangles[i].p[1].x) *
497 (pTriangles[i].p[2].x - pTriangles[i].p[1].x) +
498 (pTriangles[i].p[2].y - pTriangles[i].p[1].y) *
499 (pTriangles[i].p[2].y - pTriangles[i].p[1].y) +
500 (pTriangles[i].p[2].z - pTriangles[i].p[1].z) *
501 (pTriangles[i].p[2].z - pTriangles[i].p[1].z);
502 d3 = (pTriangles[i].p[2].x - pTriangles[i].p[0].x) *
503 (pTriangles[i].p[2].x - pTriangles[i].p[0].x) +
504 (pTriangles[i].p[2].y - pTriangles[i].p[0].y) *
505 (pTriangles[i].p[2].y - pTriangles[i].p[0].y) +
506 (pTriangles[i].p[2].z - pTriangles[i].p[0].z) *
507 (pTriangles[i].p[2].z - pTriangles[i].p[0].z);
508
509 // skip excursions completely across the cell when a contour reaches the edge
510 if ((d1 > 5.*snx*snx* step[0] ) || (d2 > 5.*sny*sny* step[1]) ||
511 (d3 > 5.*snz*snz* step[2]))
512 continue;
513
514 if (ContourFacet(pTriangles[i],&p0,&p1) == 2) {
515 glVertex3f(p0.x,p0.y,p0.z);
516 glVertex3f(p1.x,p1.y,p1.z);
517
518 if (doVrml) {
519 if (!Vrml2) {
520 fprintf (drvui->fpoutv, " Separator {\n");
521 fprintf (drvui->fpoutv, " Material { diffuseColor %s }\n", string);
522 fprintf (drvui->fpoutv, "Coordinate3 { point [\n");
523 } else {
524 fprintf (drvui->fpoutv, " Shape {\n");
525 fprintf (drvui->fpoutv,
526 " geometry IndexedLineSet { coord Coordinate{ point [\n");
527 }
528 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f,\n", p0.x,p0.y,p0.z);
529 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f]\n", p1.x,p1.y,p1.z);
530 if (Vrml2) {
531 fprintf (drvui->fpoutv,
532 " }\n coordIndex [ 0,1,-1]\n color Color { color [%s]"
533 "}\n colorIndex [0]\n colorPerVertex FALSE\n }\n }\n",
534 string);
535 } else {
536 fprintf (drvui->fpoutv,
537 "}\n IndexedLineSet { coordIndex [0,1,-1] }\n}\n");
538 }
539 }
540 if (doPOV) {
541 if (fabs(p0.x-p1.x) > 0.00001f || fabs(p0.y-p1.y) > 0.00001f
542 || fabs(p0.z-p1.z) >0.00001f) {
543 fprintf (drvui->fpoutp, " cylinder { <%f,%f,%f>,<%f,%f,%f>,%f\n",
544 p0.x, p0.y,p0.z, p1.x,p1.y,p1.z,
545 0.0001 * Scale);
546 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n }\n",
547 Color2);
548 }
549 }
550 if (doAsy) {
551 if (fabs(p0.x-p1.x) > 0.00001f || fabs(p0.y-p1.y) > 0.00001f
552 || fabs(p0.z-p1.z) >0.00001f) {
553 fprintf (drvui->fpouta, " draw(pic, (<%8.5f,%8.5f,%8.5f)--(%8.5f,%8.5f,%8.5f),\n",
554 p0.x, p0.y,p0.z, p1.x,p1.y,p1.z);
555 fprintf (drvui->fpouta,"rgb(%4.2f,%4.2f,%4.2f) );\n",Red,Green,Blue);
556 }
557 }
558 }
559 }
560 glEnd();
561 glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);
562 glEnable (GL_LIGHTING);
563 }
564 glPopMatrix ();
565 free (pTriangles); /* free the contour array */
566 drvui->mainWindow->cursor (FL_CURSOR_DEFAULT);
567 }
568
569 /* interpolate rho on a 3-D grid */
570
571 float
InterpolateMap(float x,float y,float z)572 InterpolateMap (float x, float y, float z)
573 {
574
575 /* for simplicity, use weighting appropriate for a cubic grid (should be better than no interpolation) */
576
577 int ix, iy, iz;
578
579 int ix1, iy1, iz1;
580
581 float Xmin, Ymin, Zmin, Xmax, Ymax, Zmax, AvgRho;
582
583 /* Compute the grid points on either size of x, then map those points onto the 3D array indices */
584
585 ix = ((int) (x * mapstep_a));
586 Xmin = ((float) ix) / mapstep_a;
587 Xmax = ((float) (ix + 1)) / mapstep_a;
588 ix = ((int) (x * mapstep_a)) % mapstep_a;
589 ix1 = (ix + 1) % mapstep_a;
590
591 /* ditto for y & z */
592
593 iy = ((int) (y * mapstep_b));
594 Ymin = ((float) iy) / mapstep_b;
595 Ymax = ((float) (iy + 1)) / mapstep_b;
596 iy = ((int) (y * mapstep_b)) % mapstep_b;
597 iy1 = (iy + 1) % mapstep_b;
598 iz = ((int) (z * mapstep_c));
599 Zmin = ((float) iz) / mapstep_c;
600 Zmax = ((float) (iz + 1)) / mapstep_c;
601 iz = ((int) (z * mapstep_c)) % mapstep_c;
602 iz1 = (iz + 1) % mapstep_c;
603
604 AvgRho = FourierPt[LookupMap (ix, iy, iz)] * (Xmax - x) * (Ymax - y) * (Zmax - z);
605 AvgRho += FourierPt[LookupMap (ix1, iy, iz)] * (x - Xmin) * (Ymax - y) * (Zmax - z);
606 AvgRho += FourierPt[LookupMap (ix1, iy1, iz)] * (x - Xmin) * (y - Ymin) * (Zmax - z);
607 AvgRho += FourierPt[LookupMap (ix1, iy, iz1)] * (x - Xmin) * (Ymax - y) * (z - Zmin);
608 AvgRho += FourierPt[LookupMap (ix, iy1, iz)] * (Xmax - x) * (y - Ymin) * (Zmax - z);
609 AvgRho += FourierPt[LookupMap (ix, iy1, iz1)] * (Xmax - x) * (y - Ymin) * (z - Zmin);
610 AvgRho += FourierPt[LookupMap (ix, iy, iz1)] * (Xmax - x) * (Ymax - y) * (z - Zmin);
611 AvgRho += FourierPt[LookupMap (ix1, iy1, iz1)] * (x - Xmin) * (y - Ymin) * (z - Zmin);
612
613 AvgRho /= (Xmax - Xmin) * (Ymax - Ymin) * (Zmax - Zmin);
614
615 return AvgRho;
616 }
617
618 /* lookup a point in the Fourier map */
619
620 int
LookupMap(int ix,int iy,int iz)621 LookupMap (int ix, int iy, int iz)
622 {
623
624 return abs (ix * mapstep_b * mapstep_c + iy * mapstep_c +
625 ((mapstep_c + iz % mapstep_c) % mapstep_c));
626 }
627
628 int
Maximize_rho(int Sense)629 Maximize_rho (int Sense)
630 {
631 // find min/max in rho - Sense is -1 if min, +1 if max
632 float rhomax = -1.0e15f;
633
634 int i, j, k, l;
635
636 float w0[3][3][3];
637
638 int maxpt;
639
640 float param[7];
641
642 int map[3];
643
644 if (!ReadFourMap)
645 return (-1);
646 // extract 3x3x3 array of points around the current position
647 for (i = 0; i < 7; i++)
648 param[i] = 0.0f;
649 map[0] = (int) ((float) mapstep_a * cur_cen[0] + 0.5);
650 map[1] = (int) ((float) mapstep_b * cur_cen[1] + 0.5);
651 map[2] = (int) ((float) mapstep_c * cur_cen[2] + 0.5);
652 while (1) {
653 rhomax = -1.0e15f;
654 maxpt = -1;
655 for (i = 0; i < 3; i++) {
656 for (j = 0; j < 3; j++) {
657 for (k = 0; k < 3; k++) {
658 l = (map[0] + i - 1) * mapstep_b * mapstep_c
659 + (map[1] + j - 1) * mapstep_c + map[2] + k - 1;
660 w0[i][j][k] = (float) Sense *FourierPt[l];
661
662 if (w0[i][j][k] > rhomax) {
663 rhomax = w0[i][j][k];
664 maxpt = (i * 100) + (j * 10) + k;
665 }
666 }
667 }
668 }
669 if (maxpt == 111)
670 break; // maximum in middle
671 k = maxpt % 10; // shift maximum to middle
672 j = maxpt / 10 % 10;
673 i = maxpt / 100;
674 map[2] += k - 1;
675 map[1] += j - 1;
676 map[0] += i - 1;
677 i = 0;
678 if (map[0] < 1) {
679 map[0] = 1;
680 i = 1;
681 }
682 if (map[1] < 1) {
683 map[1] = 1;
684 i = 1;
685 }
686 if (map[2] < 1) {
687 map[2] = 1;
688 i = 1;
689 }
690 if (map[0] > mapstep_a - 1) {
691 map[0] = mapstep_a - 1;
692 i = 1;
693 }
694 if (map[1] > mapstep_b - 1) {
695 map[1] = mapstep_b - 1;
696 i = 1;
697 }
698 if (map[2] > mapstep_c - 1) {
699 map[2] = mapstep_c - 1;
700 i = 1;
701 }
702 if (i)
703 break;
704 cur_cen[2] = map[2] / (float) mapstep_c;
705 cur_cen[1] = map[1] / (float) mapstep_b;
706 cur_cen[0] = map[0] / (float) mapstep_a;
707 }
708 param[3] = -w0[1][1][1] + (w0[0][1][1] + w0[2][1][1]) * 0.5f;
709 param[4] = -w0[1][1][1] + (w0[1][0][1] + w0[1][2][1]) * 0.5f;
710 param[5] = -w0[1][1][1] + (w0[1][1][0] + w0[1][1][2]) * 0.5f;
711 if (param[3] != 0.0f)
712 param[0] = (w0[0][1][1] - w0[2][1][1]) * 0.25f / param[3];
713 if (param[4] != 0.0f)
714 param[1] = (w0[1][0][1] - w0[1][2][1]) * 0.25f / param[4];
715 if (param[5] != 0.0f)
716 param[2] = (w0[1][1][0] - w0[1][1][2]) * 0.25f / param[5];
717 param[6] = w0[1][1][1] - param[3] * param[0] * param[0]
718 - param[4] * param[1] * param[1] - param[5] * param[2] * param[2];
719
720 cur_cen[0] += param[0] / mapstep_a;
721 cur_cen[1] += param[1] / mapstep_b;
722 cur_cen[2] += param[2] / mapstep_c;
723 i = 0;
724 for (k = 0; k < nvert; k++) {
725 if (fabs (o_vert[3 * k] - cur_cen[0]) + fabs (o_vert[3 * k + 1] - cur_cen[1])
726 + fabs (o_vert[3 * k + 2] - cur_cen[2]) < 0.001)
727 i = k;
728 }
729 if (i)
730 return i; // position already in vertex list
731 drvui->orig_atom_no[nvert] = -1; // new position - add it
732 for (k = 0; k < 3; k++) {
733 o_vert[3 * nvert + k] = cur_cen[k];
734 s_vert[3 * nvert + k] = 0.0f;
735 for (l = 0; l < 3; l++) { /* calculate cartesian coordinates */
736 s_vert[3 * nvert + k] +=
737 (float) drvui->b_mat[k][l] * (cur_cen[l] - origin[l]);
738 }
739 }
740 vert_sym_no[nvert++] = 0;
741 if (!check_vert_alloc (nvert, 1))
742 return (-1);
743 return (nvert - 1);
744 }
745
746 // 3D variant of Paul Bourke's conrec taken from
747 // http://local.wasp.uwa.edu.au/~pbourke/papers/conrec/
748 //
749 //-------------------------------------------------------------------------
750 // Create a contour slice through a 3 vertex facet "p"
751 // Given the normal of the cutting plane "n" and a point on the plane "p0"
752 // Return
753 // 0 if the contour plane doesn't cut the facet
754 // 2 if it does cut the facet, the contour line segment is p1->p2
755 // -1 for an unexpected occurence
756 // If a vertex touches the contour plane nothing need to be drawn!?
757 // Note: the following has been written as a "stand along" piece of
758 // code that will work but is far from efficient....
759 //
ContourFacet(TRIANGLE tri,mpVector * p1,mpVector * p2)760 int ContourFacet (TRIANGLE tri, mpVector *p1, mpVector *p2)
761 {
762 #define SIGN(x) (x>0? 1: 0)
763 double A,B,C,D;
764 double side[3];
765
766 // Determine the equation of the plane as
767 // Ax + By + Cz + D = 0
768 A = drvui->frames[drvui->frame_no].planeeq[0];
769 B = drvui->frames[drvui->frame_no].planeeq[1];
770 C = drvui->frames[drvui->frame_no].planeeq[2];
771 D = drvui->frames[drvui->frame_no].planeeq[3];
772
773 // Evaluate the equation of the plane for each vertex
774 // If side < 0 then it is on the side to be retained
775 // else it is to be clipped
776
777 side[0] = A*tri.p[0].x + B*tri.p[0].y + C*tri.p[0].z + D;
778 side[1] = A*tri.p[1].x + B*tri.p[1].y + C*tri.p[1].z + D;
779 side[2] = A*tri.p[2].x + B*tri.p[2].y + C*tri.p[2].z + D;
780
781 // Are all the vertices on the same side
782 if (side[0] >= 0 && side[1] >= 0 && side[2] >= 0)
783 return(0);
784 if (side[0] <= 0 && side[1] <= 0 && side[2] <= 0)
785 return(0);
786
787 // Is p0 the only point on a side by itself
788 if ((SIGN(side[0]) != SIGN(side[1])) && (SIGN(side[0]) != SIGN(side[2]))) {
789 p1->x = (float)(tri.p[0].x - side[0] * (tri.p[2].x - tri.p[0].x) / (side[2] - side[0]));
790 p1->y = (float)(tri.p[0].y - side[0] * (tri.p[2].y - tri.p[0].y) / (side[2] - side[0]));
791 p1->z = (float)(tri.p[0].z - side[0] * (tri.p[2].z - tri.p[0].z) / (side[2] - side[0]));
792 p2->x = (float)(tri.p[0].x - side[0] * (tri.p[1].x - tri.p[0].x) / (side[1] - side[0]));
793 p2->y = (float)(tri.p[0].y - side[0] * (tri.p[1].y - tri.p[0].y) / (side[1] - side[0]));
794 p2->z = (float)(tri.p[0].z - side[0] * (tri.p[1].z - tri.p[0].z) / (side[1] - side[0]));
795 return(2);
796 }
797
798 // Is p1 the only point on a side by itself
799 if ((SIGN(side[1]) != SIGN(side[0])) && (SIGN(side[1]) != SIGN(side[2]))) {
800 p1->x = (float)(tri.p[1].x - side[1] * (tri.p[2].x - tri.p[1].x) / (side[2] - side[1]));
801 p1->y = (float)(tri.p[1].y - side[1] * (tri.p[2].y - tri.p[1].y) / (side[2] - side[1]));
802 p1->z = (float)(tri.p[1].z - side[1] * (tri.p[2].z - tri.p[1].z) / (side[2] - side[1]));
803 p2->x = (float)(tri.p[1].x - side[1] * (tri.p[0].x - tri.p[1].x) / (side[0] - side[1]));
804 p2->y = (float)(tri.p[1].y - side[1] * (tri.p[0].y - tri.p[1].y) / (side[0] - side[1]));
805 p2->z = (float)(tri.p[1].z - side[1] * (tri.p[0].z - tri.p[1].z) / (side[0] - side[1]));
806 return(2);
807 }
808
809 // Is p2 the only point on a side by itself
810 if ((SIGN(side[2]) != SIGN(side[0])) && (SIGN(side[2]) != SIGN(side[1]))) {
811 p1->x = (float)(tri.p[2].x - side[2] * (tri.p[0].x - tri.p[2].x) / (side[0] - side[2]));
812 p1->y = (float)(tri.p[2].y - side[2] * (tri.p[0].y - tri.p[2].y) / (side[0] - side[2]));
813 p1->z = (float)(tri.p[2].z - side[2] * (tri.p[0].z - tri.p[2].z) / (side[0] - side[2]));
814 p2->x = (float)(tri.p[2].x - side[2] * (tri.p[1].x - tri.p[2].x) / (side[1] - side[2]));
815 p2->y = (float)(tri.p[2].y - side[2] * (tri.p[1].y - tri.p[2].y) / (side[1] - side[2]));
816 p2->z = (float)(tri.p[2].z - side[2] * (tri.p[1].z - tri.p[2].z) / (side[1] - side[2]));
817 return(2);
818 }
819
820 // Shouldn't get here
821 return(-1);
822 }
823
824 /* */
825 /* BW and "Hot-Cold" Colormaps from Paul Bourke's */
826 /* http://local.wasp.uwa.edu.au/~pbourke/texture_colour/colourramp/ */
827 /* */
colorramp(float rho,float * r,float * g,float * b)828 void colorramp(float rho,float *r, float *g, float *b)
829 {
830 float dv=Map_Info.rhomx-Map_Info.rhomn;
831
832 if (drvui->frames[drvui->frame_no].slice == 3) {
833 *r=*g=*b=(rho-Map_Info.rhomn)/(Map_Info.rhomx-Map_Info.rhomn);
834 return;
835 }
836 *r = *g = *b = 1.0f;
837 if (rho < (Map_Info.rhomn +0.25f *dv) ) {
838 *r = 0.;
839 *g = 4.0f * (rho - Map_Info.rhomn)/dv;
840 } else if (rho < (Map_Info.rhomn +0.5f * dv)) {
841 *r = 0.;
842 *b = 1.0f + 4.0f * (Map_Info.rhomn +0.25f * dv - rho) / dv;
843 } else if (rho < (Map_Info.rhomn + 0.75f * dv)) {
844 *r = 4.0f * (rho -Map_Info.rhomn - 0.5f * dv) / dv;
845 *b = 0.0f;
846 } else {
847 *g = 1.0f + 4.0f * (Map_Info.rhomn + 0.75f * dv - rho) /dv;
848 *b = 0.0f;
849 }
850 }
851
852 void
Add_mapslice(int type)853 Add_mapslice (int type)
854 {
855 // add a 2d map section through the last three atoms
856
857 float p[3],q[3],pq[3];
858
859 if (!ReadFourMap)
860 return;
861 if (cur_atom[3]>0) {
862 p[0]=o_vert[3*cur_atom[1]] - cur_cen[0];
863 p[1]=o_vert[3*cur_atom[1]+1] - cur_cen[1];
864 p[2]=o_vert[3*cur_atom[1]+2] - cur_cen[2];
865 q[0]=o_vert[3*cur_atom[2]] - cur_cen[0];
866 q[1]=o_vert[3*cur_atom[2]+1] - cur_cen[1];
867 q[2]=o_vert[3*cur_atom[2]+2] - cur_cen[2];
868 } else {
869 p[0]=o_vert[3*cur_atom[0]] - cur_cen[0];
870 p[1]=o_vert[3*cur_atom[0]+1] - cur_cen[1];
871 p[2]=o_vert[3*cur_atom[0]+2] - cur_cen[2];
872 q[0]=o_vert[3*cur_atom[1]] - cur_cen[0];
873 q[1]=o_vert[3*cur_atom[1]+1] - cur_cen[1];
874 q[2]=o_vert[3*cur_atom[1]+2] - cur_cen[2];
875 }
876 vcross(p,q,pq);
877 vnormalize(pq);
878 drvui->frames[drvui->frame_no].mapslice[0]=cur_cen[0];
879 drvui->frames[drvui->frame_no].mapslice[1]=cur_cen[1];
880 drvui->frames[drvui->frame_no].mapslice[2]=cur_cen[2];
881 drvui->frames[drvui->frame_no].mapnorm[0]=pq[0];
882 drvui->frames[drvui->frame_no].mapnorm[1]=pq[1];
883 drvui->frames[drvui->frame_no].mapnorm[2]=pq[2];
884 drvui->frames[drvui->frame_no].slice=type;
885 drvui->Str_File_Changed = 1;
886 Update_Str (0);
887 Generate_Drawing (0);
888 }
889
Mul_Rv(float R[3][3],float inp[3],float oup[3])890 void Mul_Rv(float R[3][3], float inp[3], float oup[3])
891 {
892 /* Multiply R * inp => oup */
893 int i, j;
894
895 for (i = 0; i < 3; i++) {
896 oup[i] = 0.0f;
897 for (j = 0; j < 3; j++)
898 oup[i] += R[i][j] * inp[j];
899 }
900 }
901
Mul_Rinvv(float R[3][3],float inp[3],float oup[3])902 void Mul_Rinvv(float R[3][3], float inp[3], float oup[3])
903 {
904 /* Multiply RT * inp => oup */
905 int i, j;
906
907 for (i = 0; i < 3; i++) {
908 oup[i] = 0.0f;
909 for (j = 0; j < 3; j++)
910 oup[i] += R[j][i] * inp[j];
911 }
912 }
913
generate_slice(void)914 void generate_slice (void)
915 {
916 /* Calculate the density in the plane described by the planeeq information.
917 * It lies a distance planeeq[3] from the origin, and has a normal N
918 * given by planeeq[0-2].
919 * Three different coordinate systems will be used: (1) The crystallographic
920 * system x, (2) the corresponding Cartesian system X, and (3) a rotated
921 * system X' where the slice is perpendicular to the Z' axis.
922 * Routine Convert_Cryst_Cart() converts from crystal to Cartesian.
923 * Routine Convert_Cart_Cryst() converts from Cartesian to crystal.
924 * Matrix R converts from X to X' and its transpose converts from X' to X.
925 */
926 int i,j,k,l;
927 int dx,dy,dxy;
928 float pp[4][3];
929 float q[4][3];
930 float rho[4],rhor[4],rhog[4],rhob[4];
931 float *prhor = NULL, *prhog = NULL, *prhob = NULL;
932 int vertn = 0;
933 double phi, chi, cosp, sinp, cosc, sinc;
934 float min_x = 999999.0f;
935 float min_y = 999999.0f;
936 float max_x = -999999.0f;
937 float max_y = -999999.0f;
938 float xlim[2], ylim[2], zlim[2];
939 float R[3][3], X[3], XP[3], x[3];
940 float bmat_inv[3][3];
941 int offset[4][2] = {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}};
942 float D;
943
944 drvui->mainWindow->cursor (FL_CURSOR_WAIT);
945 for (i = 0; i < 3; i++)
946 for (j = 0; j < 3; j++)
947 bmat_inv[i][j] = (float)drvui->b_mat[i][j];
948 matinv(bmat_inv);
949 /* The plane normal should be of unit length, but make sure */
950 vnormalize(drvui->frames[drvui->frame_no].planeeq);
951 D = drvui->frames[drvui->frame_no].planeeq[3];
952 /* Build R */
953 phi = atan2(drvui->frames[drvui->frame_no].planeeq[1],
954 drvui->frames[drvui->frame_no].planeeq[0]);
955 cosp = (float)cos(phi);
956 sinp = (float)sin(phi);
957 chi = atan2((cosp * drvui->frames[drvui->frame_no].planeeq[0] +
958 sinp * drvui->frames[drvui->frame_no].planeeq[1]),
959 -drvui->frames[drvui->frame_no].planeeq[2]);
960 sinc = (float)sin(chi);
961 cosc = (float)cos(chi);
962 R[0][0] = (float)(cosc * cosp);
963 R[0][1] = (float)(cosc * sinp);
964 R[0][2] = (float)sinc;
965 R[1][0] = -(float)sinp;
966 R[1][1] = (float)cosp;
967 R[1][2] = 0.0f;
968 R[2][0] = -(float)(sinc * cosp);
969 R[2][1] = -(float)(sinc * sinp);
970 R[2][2] = (float)cosc;
971 /* find extreme limits in X' for the display box in x */
972 xlim[0] = drvui->frames[drvui->frame_no].cryst_lim[0];
973 xlim[1] = drvui->frames[drvui->frame_no].cryst_lim[3];
974 ylim[0] = drvui->frames[drvui->frame_no].cryst_lim[1];
975 ylim[1] = drvui->frames[drvui->frame_no].cryst_lim[4];
976 zlim[0] = drvui->frames[drvui->frame_no].cryst_lim[2];
977 zlim[1] = drvui->frames[drvui->frame_no].cryst_lim[5];
978 for (i = 0; i < 2; i++) {
979 x[0] = xlim[i];
980 for (j = 0; j < 2; j++) {
981 x[1] = ylim[j];
982 for (k = 0; k < 2; k++) {
983 x[2] = zlim[k];
984 Convert_Cryst_Cart(drvui->b_mat, x, X, origin);
985 Mul_Rv(R, X, XP);
986 if (XP[0] < min_x)
987 min_x = XP[0];
988 if (XP[0] > max_x)
989 max_x = XP[0];
990 if (XP[1] < min_y)
991 min_y = XP[1];
992 if (XP[1] > max_y)
993 max_y = XP[1];
994 }
995 }
996 }
997
998 dx = 50 * (int)(max_x - min_x);
999 dy = 50 * (int)(max_y - min_y);
1000 if (dx > dy)
1001 dy = dx;
1002 else
1003 dx = dy;
1004 fprintf(drvui->flout, "No. of points in 2D section in each direction %d\n", dx);
1005 min_x = 3.0f * min_x;
1006 max_x = 3.0f * max_x;
1007 min_y = 3.0f * min_y;
1008 max_y = 3.0f * max_y;
1009
1010 // dxy = dx * dy * 4; // theoretical maximum if no clipping occurs
1011
1012 dxy=0;
1013 for (i = 0; i <= dx; i++) {
1014 XP[0] = min_x + (float)(i) * (max_x - min_x) / (float)(dx);
1015 for (j = 0; j <= dy; j++) {
1016 XP[1] = min_y + (float)(j) * (max_y - min_y) / (float)(dy);
1017 XP[2] = D; /* Now have slice point in XP */
1018 for (k = 0; k < 4; k++) {
1019 float XPP[3];
1020 XPP[2] = XP[2];
1021 XPP[1] = XP[1] + offset[k][1] * (max_y - min_y) / (float)(dx);
1022 XPP[0] = XP[0] + offset[k][0] * (max_x - min_x) / (float)(dy);
1023 Mul_Rinvv(R, XPP, X); /* Convert XPP to X */
1024 Convert_Cart_Cryst(bmat_inv, X, x, origin); /* and X to x */
1025 if (x[0] < drvui->frames[drvui->frame_no].cryst_lim[0] ||
1026 x[0] < drvui->frames[drvui->frame_no].map_lim[0] ||
1027 x[0] > drvui->frames[drvui->frame_no].cryst_lim[3] ||
1028 x[0] > drvui->frames[drvui->frame_no].map_lim[3])
1029 goto eout;
1030 if (x[1] < drvui->frames[drvui->frame_no].cryst_lim[1] ||
1031 x[1] < drvui->frames[drvui->frame_no].map_lim[1] ||
1032 x[1] > drvui->frames[drvui->frame_no].cryst_lim[4] ||
1033 x[1] > drvui->frames[drvui->frame_no].map_lim[4])
1034 goto eout;
1035 if (x[2] < drvui->frames[drvui->frame_no].cryst_lim[2] ||
1036 x[2] < drvui->frames[drvui->frame_no].map_lim[2] ||
1037 x[2] > drvui->frames[drvui->frame_no].cryst_lim[5] ||
1038 x[2] > drvui->frames[drvui->frame_no].map_lim[5])
1039 goto eout;
1040 }
1041 dxy += 4;
1042 eout:;
1043 }
1044 }
1045 fprintf(drvui->flout, "Allocation size for 2D variables (bytes) is %d\n", (int)(3*dxy*sizeof(float)));
1046
1047 if (dxy == 0) {
1048 Error_Box("Error: No part of the mapslice would be visible!\nChoose different parameters.");
1049 return;
1050 }
1051
1052
1053 if (doPOV || doVrml) {
1054
1055 prhor = (float *)malloc(dxy*sizeof(float));
1056 if (!prhor) {
1057 Error_Box("Error: Unable to allocate memory for map slice!");
1058 return;
1059 }
1060 prhog = (float *)malloc(dxy*sizeof(float));
1061 if (!prhog) {
1062 Error_Box("Error: Unable to allocate memory for map slice!");
1063 free(prhor);
1064 return;
1065 }
1066 prhob = (float *)malloc(dxy*sizeof(float));
1067 if (!prhob) {
1068 Error_Box("Error: Unable to allocate memory for map slice!");
1069 free(prhor);
1070 free(prhog);
1071 return;
1072 }
1073
1074 if (doPOV)
1075 fprintf(drvui->fpoutp,"mesh2{\n vertex_vectors { %d,\n",dxy);
1076
1077 if (doVrml) {
1078 if (!Vrml2) {
1079 fprintf (drvui->fpoutv, " Separator {\n");
1080 fprintf (drvui->fpoutv, "Coordinate3 { point [\n");
1081 } else {
1082 fprintf (drvui->fpoutv, " Shape {\n");
1083 fprintf (drvui->fpoutv,
1084 " geometry IndexedFaceSet { coord Coordinate{ point [\n");
1085 }
1086 }
1087 }
1088 glBegin(GL_TRIANGLES);
1089 glNormal3f(drvui->frames[drvui->frame_no].mapnorm[0],
1090 drvui->frames[drvui->frame_no].mapnorm[1],
1091 drvui->frames[drvui->frame_no].mapnorm[2]);
1092
1093 for (i = 0; i <= dx; i++) {
1094 XP[0] = min_x + (float)(i) * (max_x - min_x) / (float)(dx);
1095 for (j = 0; j <= dy; j++) {
1096 XP[1] = min_y + (float)(j) * (max_y - min_y) / (float)(dy);
1097 XP[2] = D; /* Now have slice point in XP */
1098 for (k = 0; k < 4; k++) {
1099 float XPP[3];
1100 XPP[2] = XP[2];
1101 XPP[1] = XP[1] + offset[k][1] * (max_y - min_y) / (float)(dx);
1102 XPP[0] = XP[0] + offset[k][0] * (max_x - min_x) / (float)(dy);
1103 Mul_Rinvv(R, XPP, X); /* Convert XPP to X */
1104 Convert_Cart_Cryst(bmat_inv, X, x, origin); /* and X to x */
1105 if (x[0] < drvui->frames[drvui->frame_no].cryst_lim[0] ||
1106 x[0] < drvui->frames[drvui->frame_no].map_lim[0] ||
1107 x[0] > drvui->frames[drvui->frame_no].cryst_lim[3] ||
1108 x[0] > drvui->frames[drvui->frame_no].map_lim[3])
1109 goto endy;
1110 if (x[1] < drvui->frames[drvui->frame_no].cryst_lim[1] ||
1111 x[1] < drvui->frames[drvui->frame_no].map_lim[1] ||
1112 x[1] > drvui->frames[drvui->frame_no].cryst_lim[4] ||
1113 x[1] > drvui->frames[drvui->frame_no].map_lim[4])
1114 goto endy;
1115 if (x[2] < drvui->frames[drvui->frame_no].cryst_lim[2] ||
1116 x[2] < drvui->frames[drvui->frame_no].map_lim[2] ||
1117 x[2] > drvui->frames[drvui->frame_no].cryst_lim[5] ||
1118 x[2] > drvui->frames[drvui->frame_no].map_lim[5])
1119 goto endy;
1120 for (l = 0; l < 3; l++) {
1121 q[k][l] = x[l];
1122 pp[k][l] = X[l];
1123 }
1124 }
1125 rho[0] = InterpolateMap(q[0][0], q[0][1], q[0][2]);
1126 colorramp(rho[0], &rhor[0], &rhog[0], &rhob[0]);
1127 rho[1] = InterpolateMap(q[1][0], q[1][1], q[1][2]);
1128 colorramp(rho[1], &rhor[1], &rhog[1], &rhob[1]);
1129 rho[2] = InterpolateMap(q[2][0], q[2][1], q[2][2]);
1130 colorramp(rho[2], &rhor[2], &rhog[2], &rhob[2]);
1131 rho[3] = InterpolateMap(q[3][0], q[3][1], q[3][2]);
1132 colorramp(rho[3], &rhor[3], &rhog[3], &rhob[3]);
1133 if (doPOV || doVrml) {
1134 prhor[vertn] = rhor[0];
1135 prhog[vertn] = rhog[0];
1136 prhob[vertn++] = rhob[0];
1137 prhor[vertn] = rhor[1];
1138 prhog[vertn] = rhog[1];
1139 prhob[vertn++] = rhob[1];
1140 prhor[vertn] = rhor[2];
1141 prhog[vertn] = rhog[2];
1142 prhob[vertn++] = rhob[2];
1143 prhor[vertn] = rhor[3];
1144 prhog[vertn] = rhog[3];
1145 prhob[vertn++] = rhob[3];
1146 }
1147
1148 glColor3f(rhor[0], rhog[0], rhob[0]);
1149 glVertex3f(pp[0][0], pp[0][1], pp[0][2]);
1150 glColor3f(rhor[1], rhog[1], rhob[1]);
1151 glVertex3f(pp[1][0],pp[1][1],pp[1][2]);
1152 glColor3f(rhor[2], rhog[2], rhob[2]);
1153 glVertex3f(pp[2][0], pp[2][1], pp[2][2]);
1154
1155 glColor3f(rhor[0], rhog[0], rhob[0]);
1156 glVertex3f(pp[0][0], pp[0][1], pp[0][2]);
1157 glColor3f(rhor[2], rhog[2], rhob[2]);
1158 glVertex3f(pp[2][0], pp[2][1], pp[2][2]);
1159 glColor3f(rhor[3], rhog[3], rhob[3]);
1160 glVertex3f(pp[3][0], pp[3][1], pp[3][2]);
1161
1162 if (doPOV) {
1163 fprintf (drvui->fpoutp, "<%8.5f, %8.5f, %8.5f> <%8.5f,"
1164 "%8.5f, %8.5f> <%8.5f, %8.5f, %8.5f> <%8.5f,%8.5f,%8.5f>\n",
1165 pp[0][0],pp[0][1],pp[0][2],pp[1][0],pp[1][1],pp[1][2],
1166 pp[2][0],pp[2][1],pp[2][2],pp[3][0],pp[3][1],pp[3][2]);
1167 }
1168 if (doAsy) {
1169 fprintf (drvui->fpouta, "draw (pic, surface ( (%8.5f, %8.5f, %8.5f)--(%8.5f,"
1170 "%8.5f, %8.5f)--(%8.5f, %8.5f, %8.5f)--(%8.5f,%8.5f,%8.5f)--cycle),\n",
1171 pp[0][0],pp[0][1],pp[0][2],pp[1][0],pp[1][1],pp[1][2],
1172 pp[2][0],pp[2][1],pp[2][2],pp[3][0],pp[3][1],pp[3][2]);
1173 fprintf(drvui->fpouta, "new pen[] {rgb(%.2f,%.2f,%.2f),rgb(%.2f,%.2f,%.2f),",
1174 rhor[0],rhog[0],rhob[0],rhor[1],rhog[1],rhob[1]);
1175 fprintf(drvui->fpouta, "rgb(%.2f,%.2f,%.2f),rgb(%.2f,%.2f,%.2f)});\n",
1176 rhor[2],rhog[2],rhob[2],rhor[3],rhog[3],rhob[3]);
1177 }
1178 if (doVrml) {
1179 fprintf (drvui->fpoutv, "%8.5f %8.5f %8.5f, %8.5f %8.5f"
1180 "%8.5f, %8.5f %8.5f %8.5f, %8.5f %8.5f %8.5f",
1181 pp[0][0],pp[0][1],pp[0][2],pp[1][0],pp[1][1],pp[1][2],
1182 pp[2][0],pp[2][1],pp[2][2],pp[3][0],pp[3][1],pp[3][2]);
1183 if (vertn < dxy-1) {
1184 fprintf(drvui->fpoutv,",\n");
1185 } else {
1186 fprintf(drvui->fpoutv,"]}\n");
1187 }
1188 }
1189 endy:;
1190 }
1191 }
1192 glEnd();
1193
1194 if (doPOV) {
1195 fprintf(drvui->fpoutp," }\n texture_list { %d,\n ",dxy);
1196 for (k=0;k<dxy-1;k++)
1197 fprintf(drvui->fpoutp," texture {pigment { color red %.2f green %.2f blue %.2f }},\n",
1198 prhor[k],prhog[k],prhob[k]);
1199 fprintf(drvui->fpoutp," texture {pigment { color red %.2f green %.2f blue %.2f }}\n }\n",
1200 prhor[dxy-1],prhog[dxy-1],prhob[dxy-1]);
1201
1202 fprintf(drvui->fpoutp," face_indices { %d,\n ",dxy/2);
1203 int kk=0;
1204 int k=0;
1205 while (kk<dxy-6) { // bracketed triple is coordinate, other texture
1206 fprintf(drvui->fpoutp,"<%d,%d,%d>,%d,%d,%d,",kk,kk+1,kk+2,kk,kk+1,kk+2);
1207 fprintf(drvui->fpoutp,"<%d,%d,%d>,%d,%d,%d,",kk,kk+2,kk+3,kk,kk+1,kk+2);
1208 kk+=4;
1209 k++;
1210 if (k==2) {
1211 fprintf(drvui->fpoutp,"\n ");
1212 k=0;
1213 }
1214 }
1215 fprintf(drvui->fpoutp,"<%d,%d,%d>,%d,%d,%d,<%d,%d,%d>,%d,%d,%d\n }\n}\n",
1216 dxy-4,dxy-3,dxy-2,dxy-4,dxy-3,dxy-2,dxy-4,dxy-2,dxy-1,dxy-4,dxy-2,dxy-1);
1217 }
1218
1219 if (doVrml) {
1220 if (!Vrml2) {
1221 fprintf (drvui->fpoutv, " MaterialBinding { value PER_VERTEX }\n Material { diffuseColor[ \n");
1222 for (k=0;k<dxy-1;k++)
1223 fprintf(drvui->fpoutv," %.2f %.2f %.2f,\n",
1224 prhor[k],prhog[k],prhob[k]);
1225 fprintf(drvui->fpoutv," %.2f %.2f %.2f] }\n",
1226 prhor[dxy-1],prhog[dxy-1],prhob[dxy-1]);
1227 fprintf(drvui->fpoutv," IndexedFaceSet{");
1228 }
1229 fprintf(drvui->fpoutv," coordIndex[");
1230 int kk=0;
1231 while (kk<dxy-6) {
1232 for (k=0;k<3;k++){
1233 fprintf(drvui->fpoutv,"%d,%d,%d,-1,%d,%d,%d,-1,",kk,kk+1,kk+2,kk,kk+2,kk+3);
1234 kk+=4;
1235 }
1236 fprintf(drvui->fpoutv,"\n ");
1237 }
1238 fprintf(drvui->fpoutv,"%d,%d,%d,-1,%d,%d,%d,-1]\n",
1239 dxy-4,dxy-3,dxy-2,dxy-4,dxy-2,dxy-1);
1240 if (Vrml2) {
1241 fprintf (drvui->fpoutv, " color Color { color[\n");
1242 for (k=0;k<dxy-1;k++)
1243 fprintf(drvui->fpoutv," %.2f %.2f %.2f,\n",
1244 prhor[k],prhog[k],prhob[k]);
1245 fprintf(drvui->fpoutv," %.2f %.2f %.2f]}\n",
1246 prhor[dxy-1],prhog[dxy-1],prhob[dxy-1]);
1247 fprintf(drvui->fpoutv," colorPerVertex TRUE\n");
1248 fprintf(drvui->fpoutv," convex TRUE\n");
1249 fprintf(drvui->fpoutv," solid FALSE\n");
1250 }
1251 fprintf(drvui->fpoutv," }\n }\n");
1252 }
1253 free (prhor);
1254 free (prhog);
1255 free (prhob);
1256 drvui->mainWindow->cursor (FL_CURSOR_DEFAULT);
1257 }
1258
MapLegend(void)1259 void MapLegend (void)
1260 {
1261 char label[80]="";
1262 char *p;
1263 GLfloat fw=0.0;
1264 float glr, glg, glb, rho;
1265 float dv;
1266 int i;
1267
1268 dv = Map_Info.rhomx - Map_Info.rhomn;
1269 glMatrixMode (GL_PROJECTION);
1270 glPushMatrix ();
1271 glLoadIdentity ();
1272 glOrtho (0, 1, 0, 1, -1, 1);
1273 glMatrixMode (GL_MODELVIEW);
1274 glPushMatrix ();
1275 glLoadIdentity ();
1276 glDisable (GL_LIGHTING);
1277 glBegin (GL_LINES);
1278 for (i = 0; i < 256; i++) {
1279 rho = Map_Info.rhomn + 0.0039f * (float)i * dv;
1280 colorramp (rho, &glr, &glg, &glb);
1281 glColor3f (glr, glg, glb);
1282 glVertex3f (0.05f, 0.55f + .00156f * (float)i, 0.0f);
1283 glVertex3f (0.1f, 0.55f + .00156f * (float)i, 0.0f);
1284 }
1285 glEnd ();
1286 glColor3f (0., 0., 0.);
1287 glBegin (GL_LINES);
1288 glVertex3f (0.03f, 0.55f + 0.48f * (0.0f - Map_Info.rhomn) / dv, 0.0f);
1289 glVertex3f (0.12f, 0.55f + 0.48f * (0.0f - Map_Info.rhomn) / dv, 0.0f);
1290 glEnd();
1291
1292 for (i = 0; i < 6; i++) {
1293 fw = (0.2f * (float)i *dv) + Map_Info.rhomn;
1294 sprintf (label, "% 5.3f", fw);
1295 glPushMatrix ();
1296 glTranslatef (0.1f, 0.55f + 0.08f * (float)i, 0.0f);
1297 glLineWidth (2.0);
1298 glScalef (drvui->label_scale * 0.0002f, drvui->label_scale * 0.0002f, drvui->label_scale * 0.0002f);
1299 for (p = label; *p; p++)
1300 glutStrokeCharacter (GLUT_STROKE_ROMAN, *p);
1301 glScalef (1.0f, 1.0f, 1.0f);
1302 glLineWidth (1.0f);
1303 glPopMatrix ();
1304 }
1305 glEnable (GL_LIGHTING);
1306 glMatrixMode (GL_PROJECTION);
1307 glPopMatrix ();
1308 glMatrixMode (GL_MODELVIEW);
1309 glPopMatrix ();
1310 }
1311