1 // $Id: Fourier.cxx 1102 2011-01-03 15:17:35Z martin $
2 //
3 // Fourier.cxx - Source module for DRAWxtl V5.5
4 // Coded using the FLTK 1.1.6 widget set
5 //
6 // Larry W. Finger, Martin Kroeker and Brian Toby
7 //
8 //
9 // This module contains the following routines:
10 //
11 // fillcube - read a 512 byte record of a dn6 Fourier map
12 // read_dn6 - read O Format (dn6) Fourier file
13 // read_fcf - read CIF-standard Fo/Fc file and compute density
14 // read_grd - read GSAS (grd) Fourier file
15 // read_stf - read JANA (stf) Fourier file
16 // read_m80 - read JANA (m80) Fourier file
17 // read_m81 - read JANA (m81) density file
18 // read_flp - read FullProf (flp) Fourier file
19 // read_w2k - read WIEN2k calculated charge density file
20 // read_exc - read EXCITING calculated charge density or ELF file
21 // read_vasp - read VASP calculated charge density file
22 // read_aim - read WIEN2k calculated Bader AIM surface file
23 // read_xsf - read charge density file in XCrysDen format
24
25 #include <stdio.h>
26 #include <string.h>
27 #include <math.h>
28 #include <stdlib.h>
29 #if defined(WIN32)
30 #define snprintf _snprintf
31 #endif
32 #include <ctype.h>
33 #include "drawxtl.h"
34 #include "draw_ext.h"
35 #include "DRAWxtlViewUI.h"
36
37 #include "DRAWxtl_proto.h"
38
39 /* ************************************************************** */
40 /* ************************************************************** */
41
42 int
fillcube(uchar * cube,FILE * mapin)43 fillcube (uchar * cube, FILE * mapin)
44 {
45 /* routine to read the next record of a dn6 map.
46 *
47 * The first record is a header. After that, each record contains
48 * an 8 x 8 x 8 piece of the map scaled so that each map point
49 * occupies a single byte.
50 */
51
52 int i, val;
53
54 for (i = 0; i < 512; i++) {
55 if ((val = fgetc (mapin)) == EOF)
56 return 0;
57 cube[i] = (uchar) val;
58 }
59 return 1;
60 }
61
62 /* ************************************************************** */
63 /* ************************************************************** */
64
65 void
read_dn6(char * infile,int Quick)66 read_dn6 (char *infile, int Quick)
67 {
68 /* routine to read a dn6 or 'O' format Fourier map */
69 FILE *mapin;
70
71 float rhomin = 1.0e15f, rhomax = -1.0e15f;
72
73 char string[200];
74
75 uchar cube[512];
76
77 float rhoscal;
78
79 int rhoadd;
80
81 int iscal1, iscal2;
82
83 int i, j, k;
84
85 int ori[3];
86
87 int nga, ngb, ngc;
88
89 int kk, jj, ii;
90
91 int ncube;
92
93 int mappos;
94
95 int nz, ny, nx;
96
97 if ((mapin = fopen (infile, "rb")) == NULL) {
98 sprintf (string, "Cannot open Fourier map (.dn6) file %s", infile);
99 Error_Box (string);
100 return;
101 }
102 if (!fillcube (cube, mapin)) {
103 Error_Box ("Error reading DN6 Fourier File.");
104 fclose (mapin);
105 return;
106 }
107 /* TODO confirm the description of ori */
108 ori[0] = (int) cube[0] * 256 + (int) cube[1]; /* get starting point of map */
109 ori[1] = (int) cube[2] * 256 + (int) cube[3];
110 ori[2] = (int) cube[4] * 256 + (int) cube[5];
111 nga = (int) cube[6] * 256 + (int) cube[7];
112 ngb = (int) cube[8] * 256 + (int) cube[9];
113 ngc = (int) cube[10] * 256 + (int) cube[11];
114 mapstep_a = (int) cube[12] * 256 + (int) cube[13];
115 mapstep_b = (int) cube[14] * 256 + (int) cube[15];
116 mapstep_c = (int) cube[16] * 256 + (int) cube[17];
117 iscal1 = (int) cube[34] * 256 + (int) cube[35];
118 map_a = (float) ((int) cube[18] * 256 + (int) cube[19]) / (float) iscal1;
119 map_b = (float) ((int) cube[20] * 256 + (int) cube[21]) / (float) iscal1;
120 map_c = (float) ((int) cube[22] * 256 + (int) cube[23]) / (float) iscal1;
121 map_alpha = (float) ((int) cube[24] * 256 + (int) cube[25]) / (float) iscal1;
122 map_beta = (float) ((int) cube[26] * 256 + (int) cube[27]) / (float) iscal1;
123 map_gamma = (float) ((int) cube[28] * 256 + (int) cube[29]) / (float) iscal1;
124 iscal2 = (int) cube[36] * 256 + (int) cube[37];
125 rhoscal = (float) ((int) cube[30] * 256 + (int) cube[31]) / (float) iscal2;
126 rhoadd = (int) cube[32] * 256 + (int) cube[33];
127 Map_Info.lat_con[0] = map_a;
128 Map_Info.lat_con[1] = map_b;
129 Map_Info.lat_con[2] = map_c;
130 Map_Info.lat_con[3] = map_alpha;
131 Map_Info.lat_con[4] = map_beta;
132 Map_Info.lat_con[5] = map_gamma;
133 Map_Info.map_int[0] = mapstep_a;
134 Map_Info.map_int[1] = mapstep_b;
135 Map_Info.map_int[2] = mapstep_c;
136 Map_Info.xlim[0] = (float) ori[0] / (float) mapstep_a;
137 Map_Info.ylim[0] = (float) ori[1] / (float) mapstep_b;
138 Map_Info.zlim[0] = (float) ori[2] / (float) mapstep_c;
139 Map_Info.xlim[1] = (float) (nga + ori[0] - 1) / (float) mapstep_a++;
140 Map_Info.ylim[1] = (float) (ngb + ori[1] - 1) / (float) mapstep_b++;
141 Map_Info.zlim[1] = (float) (ngc + ori[2] - 1) / (float) mapstep_c++;
142
143 Map_Info.info_valid = 1;
144
145 if (FourierPt)
146 free (FourierPt);
147
148 FourierPt = (float *) zalloc (mapstep_a * mapstep_b * mapstep_c * sizeof (float));
149 if (FourierPt == NULL) {
150 Error_Box ("ERROR -- Unable to allocate space for map.\n");
151 fclose (mapin);
152 return;
153 }
154
155 /* read the Rho values - map file has x fastest, and z slowest,
156 * In FourierPt, we store with z fastest */
157
158 for (k = 0; k < ngc; k = k + 8) {
159 for (j = 0; j < ngb; j = j + 8) {
160 for (i = 0; i < nga; i = i + 8) {
161 if (!fillcube (cube, mapin)) { /* read the next 8 x 8 x 8 block */
162 Error_Box ("Error reading DN6 Fourier File.");
163 fclose (mapin);
164 return;
165 }
166 for (kk = 0; kk < min (8, ngc - k); kk++) { /* unpack the cube */
167 nz = k + kk + ori[2];
168 for (jj = 0; jj < min (8, ngb - j); jj++) {
169 ny = (j + jj + ori[1]) * mapstep_c;
170 for (ii = 0; ii < min (8, nga - i); ii = ii + 2) {
171 nx = (i + ii + ori[0]) * mapstep_c * mapstep_b;
172 mappos = nx + ny + nz;
173 ncube = 64 * kk + 8 * jj + ii;
174 FourierPt[mappos] =
175 (float) (cube[ncube + 1] - rhoadd) / rhoscal;
176 if (FourierPt[mappos] > rhomax)
177 rhomax = FourierPt[mappos];
178 if (FourierPt[mappos] < rhomin)
179 rhomin = FourierPt[mappos];
180 if (ii + 1 < min (8, nga - i)) {
181 FourierPt[mappos + 1] =
182 (float) (cube[ncube] - rhoadd) / rhoscal;
183 if (FourierPt[mappos + 1] > rhomax)
184 rhomax = FourierPt[mappos + 1];
185 if (FourierPt[mappos + 1] < rhomin)
186 rhomin = FourierPt[mappos + 1];
187 }
188 }
189 }
190 }
191 }
192 }
193 }
194 if (!Quick) {
195 fprintf (drvui->flout, "Reading Fourier map file %s\n"
196 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
197 "Range of rho values: %f to %f\n", infile,
198 map_a, map_b, map_c, map_alpha, map_beta, map_gamma,
199 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax);
200 }
201 Map_Info.rhomn = rhomin;
202 Map_Info.rhomx = rhomax;
203 ReadFourMap = 1;
204 fclose (mapin);
205 }
206
207 /* ************************************************************** */
208 /* ************************************************************** */
209
210 void
read_fcf(char * infile,int Quick)211 read_fcf (char *infile, int Quick)
212 {
213 // read CIF-standard Fo/Fc file and compute density
214
215 FILE *mapin;
216
217 float *fo, *fc, *f2, *ftmp1, *ftmp2, *ftmp3;
218
219 short *ih, *ik, *il, *itmp1, *itmp2, *itmp3;
220
221 float *Cklx;
222
223 float *Dklx;
224
225 float *Elxy;
226
227 float *Flxy;
228
229 double x, y, z;
230
231 char string[200];
232
233 char maptitle[61];
234
235 int squared, phases, ab;
236
237 int nval = 200;
238
239 int nr;
240
241 int i = 0, j, k, l;
242
243 int ijk = 0;
244
245 float rhomin = 1.0e15f, rhomax = -1.0e15f;
246
247 float factor;
248
249 short ih0, ik0, il0;
250
251 short ih1, ik1, il1;
252
253 short im, in, io;
254
255 float fo0, fc0, fc1, f20, f21;
256
257 int kmin, kmax, nk;
258
259 int lmin, lmax, nl;
260
261 float phase;
262
263 int skip;
264
265 int progress;
266
267 int dimension = 3;
268
269 if (FourierPt != NULL) {
270 return;
271 }
272
273 if ((mapin = fopen (infile, "r")) == NULL) {
274 sprintf (string, "Cannot open Fo/Fc (CIF fcf) file %s", infile);
275 Error_Box (string);
276 return;
277 }
278
279 fo = (float *) zalloc (200 * sizeof (float));
280 if (fo == NULL) {
281 sprintf (string, "Unable to obtain storage for Fo data");
282 Error_Box (string);
283 fclose (mapin);
284 return;
285 }
286 fc = (float *) zalloc (200 * sizeof (float));
287 if (fc == NULL) {
288 sprintf (string, "Unable to obtain storage for Fc data");
289 Error_Box (string);
290 free (fo);
291 fclose (mapin);
292 return;
293 }
294 f2 = (float *) zalloc (200 * sizeof (float));
295 if (f2 == NULL) {
296 sprintf (string, "Unable to obtain storage for B or phase data");
297 Error_Box (string);
298 free (fo);
299 free (fc);
300 fclose (mapin);
301 return;
302 }
303 ih = (short *) zalloc (200 * sizeof (short));
304 if (ih == NULL) {
305 sprintf (string, "Unable to obtain storage for h data");
306 Error_Box (string);
307 free (fo);
308 free (fc);
309 free (f2);
310 fclose (mapin);
311 return;
312 }
313 ik = (short *) zalloc (200 * sizeof (short));
314 if (ik == NULL) {
315 sprintf (string, "Unable to obtain storage for k data");
316 Error_Box (string);
317 free (fo);
318 free (fc);
319 free (f2);
320 free (ih);
321 fclose (mapin);
322 return;
323 }
324 il = (short *) zalloc (200 * sizeof (short));
325 if (il == NULL) {
326 sprintf (string, "Unable to obtain storage for l data");
327 Error_Box (string);
328 free (fo);
329 free (fc);
330 free (f2);
331 free (ih);
332 free (ik);
333 fclose (mapin);
334 return;
335 }
336
337 do {
338 fgets (string, 80, mapin);
339 if (strstr (string, "title")) {
340 memset (maptitle, 0, 61);
341 sscanf (string, "%*s %60c", maptitle);
342 trim_string (maptitle, 60);
343 strcpy (Map_Info.title, maptitle);
344 Map_Info.info_valid = 1;
345 }
346 } while (strncmp (string, " _refln", 6) && !feof(mapin) ); // search for start of _refln group
347
348 if (strncmp (string, " _refln_index_h", 14)) {
349 sprintf (string, "Unsupported field sequence in fcf file - need h k l first");
350 Error_Box (string);
351 free (fo);
352 free (fc);
353 free (f2);
354 free (ih);
355 free (ik);
356 free (il);
357 fclose (mapin);
358 return;
359 }
360
361 squared = 0;
362 phases = 0;
363 ab = 0;
364 while (!strncmp (string, " _refln", 6)) {
365 fgets (string, 80, mapin);
366 if (strstr (string, "squared"))
367 squared = 1;
368 if (strstr (string, "_phase_"))
369 phases = 1;
370 if (strstr (string, "A_calc"))
371 ab = 1;
372 if (strstr (string, "index_m"))
373 dimension = 4;
374 if (strstr (string, "index_n"))
375 dimension = 5;
376 if (strstr (string, "index_o"))
377 dimension = 6;
378 }
379
380 if (dimension > 3) {
381 sprintf (string,
382 "%d-dimensional data encountered - only main reflections will be used",
383 dimension);
384 Error_Box (string);
385 }
386
387 if (phases == 0 && ab == 0) {
388 sprintf (string, "fcf file does not contain phase angle or A,B data");
389 Error_Box (string);
390 free (fo);
391 free (fc);
392 free (f2);
393 free (ih);
394 free (ik);
395 free (il);
396 fclose (mapin);
397 return;
398 }
399
400 nr = 0;
401 do {
402 skip = 0;
403 if (dimension > 3) {
404 switch (dimension) {
405 case 4:
406 im = 0;
407 i = sscanf (string, "%hd %hd %hd %hd %f %f %f", &ih0, &ik0, &il0, &im,
408 &fo0, &fc0, &f20);
409 if (im != 0)
410 skip = 1;
411 break;
412 case 5:
413 im = in = 0;
414 i = sscanf (string, "%hd %hd %hd %hd %hd %f %f %f", &ih0, &ik0, &il0, &im,
415 &in, &fo0, &fc0, &f20);
416 if (im != 0 || in != 0)
417 skip = 1;
418 break;
419 case 6:
420 im = in = io = 0;
421 i = sscanf (string, "%hd %hd %hd %hd %hd %hd %f %f %f",
422 &ih0, &ik0, &il0, &im, &in, &io, &fo0, &fc0, &f20);
423 if (im != 0 || in != 0 || io != 0)
424 skip = 1;
425 break;
426 }
427 } else {
428 i = sscanf (string, "%hd %hd %hd %f %*f %f %f", &ih0, &ik0, &il0, &fo0, &fc0,
429 &f20);
430 }
431 if (i < dimension + 3)
432 break;
433 if (skip == 1) {
434 fgets (string, 80, mapin);
435 continue;
436 }
437
438 if (squared == 1)
439 fo0 = (float) sqrt (fo0);
440 if (ab) { // change all data to 'phase' type
441 fc1 = (float) sqrt (fc0 * fc0 + f20 * f20);
442 f20 = (float) atan2 (f20 / fc1, fc0 / fc1);
443 fc0 = fc1;
444 } else {
445 f20 *= (float) (PI / 180.0); // phase in radians
446 }
447
448 // apply spacegroup symmetry (including xyz - catch duplicates in input list)
449 for (k = 0; k < drvui->ng; ++k) {
450 skip = 0;
451
452 ih1 =
453 drvui->ss[k][0][0] * ih0 + drvui->ss[k][1][0] * ik0 +
454 drvui->ss[k][2][0] * il0;
455 ik1 =
456 drvui->ss[k][0][1] * ih0 + drvui->ss[k][1][1] * ik0 +
457 drvui->ss[k][2][1] * il0;
458 il1 =
459 drvui->ss[k][0][2] * ih0 + drvui->ss[k][1][2] * ik0 +
460 drvui->ss[k][2][2] * il0;
461 phase =
462 2.0f * (float) PI *(drvui->ts[k][0] * ih0 + drvui->ts[k][1] * ik0 +
463 drvui->ts[k][2] * il0);
464 f21 = f20 + phase;
465 f21 = (float) atan2 (sin (f21), cos (f21)); // get f21 in range -PI to PI
466
467 for (j = 0; j < nr; j++) {
468 if (ih[j] == ih1 && ik[j] == ik1 && il[j] == il1) {
469 skip = 1;
470 break;
471 }
472 }
473 if (skip == 0) {
474 ih[nr] = ih1;
475 ik[nr] = ik1;
476 il[nr] = il1;
477 fo[nr] = fo0;
478 fc[nr] = fc0;
479 f2[nr] = f21;
480 nr++;
481
482 }
483 ih1 *= -1; // add Friedel opposite
484 ik1 *= -1;
485 il1 *= -1;
486 f21 *= -1.0;
487 skip = 0;
488
489 for (j = 0; j < nr; j++) {
490 if (ih[j] == ih1 && ik[j] == ik1 && il[j] == il1) {
491 skip = 1;
492 break;
493 }
494 }
495 if (skip == 0) {
496 ih[nr] = ih1;
497 ik[nr] = ik1;
498 il[nr] = il1;
499 fo[nr] = fo0;
500 fc[nr] = fc0;
501 f2[nr] = f21;
502 nr++;
503 }
504
505 if (nr > nval - 10) {
506 nval += 100;
507 ftmp1 = (float *) realloc (fo, nval * sizeof (float));
508 ftmp2 = (float *) realloc (fc, nval * sizeof (float));
509 ftmp3 = (float *) realloc (f2, nval * sizeof (float));
510 itmp1 = (short *) realloc (ih, nval * sizeof (short));
511 itmp2 = (short *) realloc (ik, nval * sizeof (short));
512 itmp3 = (short *) realloc (il, nval * sizeof (short));
513 if (!ftmp1 || !ftmp2 || !ftmp3 || !itmp1 || !itmp2 || !itmp3) {
514 sprintf (string,
515 "Unable to expand storage for h, k, l, fo, fc or phase data");
516 free (fo);
517 free (fc);
518 free (f2);
519 free (ih);
520 free (ik);
521 free (il);
522 if (ftmp1)
523 free (ftmp1);
524 if (ftmp2)
525 free (ftmp2);
526 if (ftmp3)
527 free (ftmp3);
528 if (itmp1)
529 free (itmp1);
530 if (itmp2)
531 free (itmp2);
532 if (itmp3)
533 free (itmp3);
534 Error_Box (string);
535 fclose (mapin);
536 return;
537 } else {
538 fo = ftmp1;
539 fc = ftmp2;
540 f2 = ftmp3;
541 ih = itmp1;
542 ik = itmp2;
543 il = itmp3;
544 }
545 } // if we need to realloc
546 } // loop over all symmetry operators
547 fgets (string, 80, mapin);
548 } while (i > 0 && !feof (mapin));
549
550 fclose (mapin);
551
552
553 factor = 1.0f / (drvui->lat_con[0] * drvui->lat_con[1] * drvui->lat_con[2]
554 * (float) sqrt (1.0 - cos (drvui->lat_con[3] * PI / 180.0) *
555 cos (drvui->lat_con[3] * PI / 180.0) *
556 cos (drvui->lat_con[4] * PI / 180.0) *
557 cos (drvui->lat_con[4] * PI / 180.0) *
558 cos (drvui->lat_con[5] * PI / 180.0) *
559 cos (drvui->lat_con[5] * PI / 180.0)));
560 kmin = kmax = ik[0];
561 lmin = lmax = il[0];
562
563 for (l = 0; l < nr; l++) {
564 if (ik[l] < kmin)
565 kmin = ik[l]; // find range of k and l
566 if (ik[l] > kmax)
567 kmax = ik[l];
568 if (il[l] < lmin)
569 lmin = il[l];
570 if (il[l] > lmax)
571 lmax = il[l];
572 switch (Map_Info.map_type) { // get map coefficients/V
573 default:
574 case 0: // Fo map
575 break;
576 case 1: // Fc map
577 fo[l] = fc[l];
578 break;
579 case 2: // Fo - Fc map
580 fo[l] = fo[l] - fc[l];
581 break;
582 case 3: // 2Fo - Fc map
583 fo[l] = 2.0f * fo[l] - fc[l];
584 break;
585 case 4: // Fo2 (Patterson)
586 fo[l] = fo[l] * fo[l];
587 f2[l] = 0.0f;
588 break;
589 }
590 fo[l] *= factor; // update coefficients times 1/V
591 }
592
593 nk = kmax - kmin + 1;
594 nl = lmax - lmin + 1;
595 Cklx = (float *) zalloc (nk * nl * sizeof (float));
596 Dklx = (float *) zalloc (nk * nl * sizeof (float));
597 Elxy = (float *) zalloc (nl * sizeof (float));
598 Flxy = (float *) zalloc (nl * sizeof (float));
599
600 if (!Cklx || !Dklx || !Elxy || !Flxy) {
601 Error_Box ("Unable to allocate Beevers-Lipson arrays.");
602 if (Cklx) free (Cklx);
603 if (Dklx) free (Dklx);
604 if (Elxy) free (Elxy);
605 if (Flxy) free (Flxy);
606 return;
607 }
608 mapstep_a = (int) ((float)Map_Info.res * drvui->lat_con[0] + 0.5f);
609 mapstep_b = (int) ((float)Map_Info.res * drvui->lat_con[1] + 0.5f);
610 mapstep_c = (int) ((float)Map_Info.res * drvui->lat_con[2] + 0.5f);
611
612 FourierPt = (float *) zalloc (mapstep_a * mapstep_b * mapstep_c * sizeof (float));
613 if (FourierPt == NULL) {
614 Error_Box ("ERROR -- Unable to allocate space for map.");
615 return;
616 }
617 // start the Beevers-Lipson expansion
618
619 Progress_Window (-1, "Computing Density", (float) mapstep_a * mapstep_b);
620
621 progress = 0;
622 for (i = 0; i < mapstep_a; i++) {
623 x = (float) i / (float) mapstep_a;
624 memset (Cklx, 0, nk * nl * sizeof (float));
625 memset (Dklx, 0, nk * nl * sizeof (float));
626 for (l = 0; l < nr; l++) {
627 ijk = (il[l] - lmin) * nk + ik[l] - kmin;
628 Cklx[ijk] += fo[l] * (float) cos (2.0 * PI * ih[l] * x - f2[l]);
629 Dklx[ijk] -= fo[l] * (float) sin (2.0 * PI * ih[l] * x - f2[l]);
630 }
631 for (j = 0; j < mapstep_b; j++) {
632 Progress_Window (0, NULL, (float) ++progress);
633 y = (float) j / (float) mapstep_b;
634 memset (Elxy, 0, nl * sizeof (float));
635 memset (Flxy, 0, nl * sizeof (float));
636 for (l = 0; l < nl; l++) {
637 for (k = 0; k < nk; k++) {
638 ijk = l * nk + k;
639 Elxy[l] += Cklx[ijk] * (float) cos (2.0 * PI * (k + kmin) * y)
640 + Dklx[ijk] * (float) sin (2.0 * PI * (k + kmin) * y);
641 Flxy[l] += Dklx[ijk] * (float) cos (2.0 * PI * (k + kmin) * y)
642 - Cklx[ijk] * (float) sin (2.0 * PI * (k + kmin) * y);
643 }
644 }
645 for (k = 0; k < mapstep_c; k++) {
646 z = (float) k / (float) mapstep_c;
647 ijk = i * mapstep_b * mapstep_c + j * mapstep_c + k;
648
649 for (l = 0; l < nl; l++) {
650 FourierPt[ijk] += Elxy[l] * (float) cos (2.0 * PI * (l + lmin) * z)
651 + Flxy[l] * (float) sin (2.0 * PI * (l + lmin) * z);
652 }
653 if (FourierPt[ijk] < rhomin)
654 rhomin = FourierPt[ijk];
655 if (FourierPt[ijk] > rhomax)
656 rhomax = FourierPt[ijk];
657 }
658 }
659 }
660 if (!Quick)
661 fprintf (drvui->fcns, "Reading Fourier map file %s\n "
662 "map steps: a=%d b=%d c=%d\n "
663 "Range of rho values: %f to %f, gridsize %d\n", infile,
664 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax, ijk);
665
666
667 Progress_Window (-2, NULL, 0.0f);
668 ReadFourMap = 1;
669 Map_Info.rhomn = rhomin;
670 Map_Info.rhomx = rhomax;
671 Map_Info.map_int[0] = mapstep_a;
672 Map_Info.map_int[1] = mapstep_b;
673 Map_Info.map_int[2] = mapstep_c;
674 Map_Info.lat_con[0] = drvui->lat_con[0];
675 Map_Info.lat_con[1] = drvui->lat_con[1];
676 Map_Info.lat_con[2] = drvui->lat_con[2];
677 Map_Info.lat_con[3] = drvui->lat_con[3];
678 Map_Info.lat_con[4] = drvui->lat_con[4];
679 Map_Info.lat_con[5] = drvui->lat_con[5];
680 Map_Info.xlim[0] = 0.0f;
681 Map_Info.xlim[1] = 1.0f;
682 Map_Info.ylim[0] = 0.0f;
683 Map_Info.ylim[1] = 1.0f;
684 Map_Info.zlim[0] = 0.0f;
685 Map_Info.zlim[1] = 1.0f;
686
687 free (Cklx);
688 free (Dklx);
689 free (Elxy);
690 free (Flxy);
691 free (fo);
692 free (fc);
693 free (f2);
694 free (ih);
695 free (ik);
696 free (il);
697 }
698
699 /* ************************************************************** */
700 /* ************************************************************** */
701
702 void
read_m81(char * infile,int Quick)703 read_m81 (char *infile, int Quick)
704 {
705 // read JANA binary .m81 Fo or jpdf file
706 FILE *mapin;
707
708 char string[200];
709
710 int iheaders[6];
711
712 float fheaders[6];
713
714 int i, j, k, l, ijk;
715
716 int nskip, nskip5, nskip6;
717
718 int nrec, nmap, maptype;
719
720 int astart, aend, bstart, bend, cstart, cend;
721
722 int step1, step2, step3, step4, step5, step6;
723
724 int new_mapstep_a, new_mapstep_b, new_mapstep_c;
725
726 int ireslt = 0;
727
728 int fullcell;
729
730 float min1, max1, min2, max2, min3, max3;
731
732 float min4, max4, min5, max5, min6, max6;
733
734 float xmin, xmax, ymin, ymax, zmin, zmax;
735
736 float x4min = 0.f, x4max = 0.f, x5min = 0.f, x5max = 0.f, x6min = 0.f, x6max = 0.f;
737
738 float xstep4, xstep5, xstep6;
739
740 float rhomin, rhomax;
741
742 float *ftmp;
743
744 int axis[6], axes, modaxes;
745 const char *theaxes[] = { " ", "x", "y", "z", "x4", "x5", "x6" };
746 const char *maptypes[] = { "Patterson", "checking Patterson", "difference Patterson",
747 "Fourier", "checking Fourier", "difference Fourier", "shape function"
748 };
749
750 if (FourierPt)
751 free (FourierPt);
752
753 if ((mapin = fopen (infile, "rb")) == NULL) {
754 sprintf (string, "Cannot open JANA m81 file %s", infile);
755 Error_Box (string);
756 return;
757 }
758 i = fread (iheaders, 4, 6, mapin); // nx ny nz nx4 nx5 nx6
759 step1 = iheaders[0];
760 step2 = iheaders[1];
761 step3 = iheaders[2];
762 step4 = iheaders[3];
763 step5 = iheaders[4];
764 step6 = iheaders[5];
765 i = fread (iheaders, 4, 2, mapin); // nx*ny nmaps
766 nrec = iheaders[0];
767 nmap = iheaders[1];
768 i = fread (fheaders, 4, 6, mapin); //xmin..xmax...zmin..zmax
769 min1 = fheaders[0];
770 max1 = fheaders[1];
771 min2 = fheaders[2];
772 max2 = fheaders[3];
773 min3 = fheaders[4];
774 max3 = fheaders[5];
775 i = fread (fheaders, 4, 6, mapin); // qmin..qmax...smin..smax (ignored)
776 min4 = fheaders[0];
777 max4 = fheaders[1];
778 min5 = fheaders[2];
779 max5 = fheaders[3];
780 min6 = fheaders[4];
781 max6 = fheaders[5];
782 i = fread (fheaders, 4, 6, mapin); // step sizes (ignored)
783 xstep4 = fheaders[3];
784 xstep5 = fheaders[4];
785 xstep6 = fheaders[5];
786 if (max4 == 0.0f && xstep4 == 1.0f)
787 xstep4 = 0.0f;
788 if (max5 == 0.0f && xstep5 == 1.0f)
789 xstep5 = 0.0f;
790 if (max6 == 0.0f && xstep6 == 1.0f)
791 xstep6 = 0.0f;
792 i = fread (iheaders, 4, 6, mapin); // axis sequence, e.g. 2 1 3
793 axis[0] = iheaders[0];
794 axis[1] = iheaders[1];
795 axis[2] = iheaders[2];
796 axis[3] = iheaders[3];
797 axis[4] = iheaders[4];
798 axis[5] = iheaders[5];
799 axes = 100 * axis[0] + 10 * axis[1] + axis[2];
800 modaxes = 100 * axis[3] + 10 * axis[4] + axis[5];
801 i = fread (iheaders, 4, 2, mapin); // composite number, maptype
802 maptype = iheaders[0] - 1;
803
804 ftmp = (float *) zalloc (nrec * sizeof (float));
805 if (ftmp == NULL) {
806 Error_Box ("ERROR -- Unable to allocate space for map.\n");
807 fclose (mapin);
808 return;
809 }
810
811 switch (axes) {
812 case 123:
813 mapstep_a = step1;
814 mapstep_b = step2;
815 mapstep_c = step3;
816 xmin = min1;
817 xmax = max1;
818 ymin = min2;
819 ymax = max2;
820 zmin = min3;
821 zmax = max3;
822 break;
823 case 132:
824 mapstep_a = step1;
825 mapstep_b = step3;
826 mapstep_c = step2;
827 xmin = min1;
828 xmax = max1;
829 ymin = min3;
830 ymax = max3;
831 zmin = min2;
832 zmax = max2;
833 break;
834 case 312:
835 mapstep_a = step2;
836 mapstep_b = step3;
837 mapstep_c = step1;
838 xmin = min2;
839 xmax = max2;
840 ymin = min3;
841 ymax = max3;
842 zmin = min1;
843 zmax = max1;
844 break;
845 case 321:
846 mapstep_a = step3;
847 mapstep_b = step2;
848 mapstep_c = step1;
849 xmin = min3;
850 xmax = max3;
851 ymin = min2;
852 ymax = max2;
853 zmin = min1;
854 zmax = max1;
855 break;
856 case 213:
857 mapstep_a = step2;
858 mapstep_b = step1;
859 mapstep_c = step3;
860 xmin = min2;
861 xmax = max2;
862 ymin = min1;
863 ymax = max1;
864 zmin = min3;
865 zmax = max3;
866 break;
867 case 231:
868 mapstep_a = step3;
869 mapstep_b = step1;
870 mapstep_c = step2;
871 xmin = min3;
872 xmax = max3;
873 ymin = min1;
874 ymax = max1;
875 zmin = min2;
876 zmax = max2;
877 break;
878 default:
879 sprintf (string, "Primary axes must be some permutation of x y z,\n"
880 "while this file has %s %s %s %s %s %s.\n", theaxes[axis[0]],
881 theaxes[axis[1]], theaxes[axis[2]], theaxes[axis[3]], theaxes[axis[4]],
882 theaxes[axis[5]]);
883 Error_Box (string);
884 free (ftmp);
885 fclose (mapin);
886 return;
887 break;
888 }
889 switch (modaxes) {
890 case 0:
891 break;
892 case 400:
893 x4step = xstep4;
894 x4min = min4;
895 x4max = max4;
896 x5min = 0.;
897 x5max = 0.;
898 x6min = 0.;
899 x6max = 0.;
900 break;
901 case 450:
902 x4step = xstep4;
903 x4min = min4;
904 x4max = max4;
905 x5step = xstep5;
906 x5min = min5;
907 x5max = max5;
908 x6min = 0.;
909 x6max = 0.;
910 break;
911 case 540:
912 x4step = xstep5;
913 x4min = min5;
914 x4max = max5;
915 x5step = xstep4;
916 x5min = min4;
917 x5max = max4;
918 x6min = 0.;
919 x6max = 0.;
920 break;
921 case 456:
922 x4step = xstep4;
923 x4min = min4;
924 x4max = max4;
925 x5step = xstep5;
926 x5min = min5;
927 x5max = max5;
928 x6step = xstep6;
929 x6min = min6;
930 x6max = max6;
931 break;
932 case 465:
933 x4step = xstep4;
934 x4min = min4;
935 x4max = max4;
936 x5step = xstep6;
937 x5min = min6;
938 x5max = max6;
939 x6step = xstep5;
940 x6min = min5;
941 x6max = max5;
942 break;
943 case 546:
944 x4step = xstep5;
945 x4min = min5;
946 x4max = max5;
947 x5step = xstep4;
948 x5min = min4;
949 x5max = max4;
950 x6step = xstep6;
951 x6min = min6;
952 x6max = max6;
953 break;
954 case 564:
955 x4step = xstep6;
956 x4min = min6;
957 x4max = max6;
958 x5step = xstep4;
959 x5min = min4;
960 x5max = max4;
961 x6step = xstep5;
962 x6min = min5;
963 x6max = max5;
964 break;
965 case 654:
966 x4step = xstep6;
967 x4min = min6;
968 x4max = max6;
969 x5step = xstep5;
970 x5min = min5;
971 x5max = max5;
972 x6step = xstep4;
973 x6min = min4;
974 x6max = max4;
975 break;
976 case 645:
977 x4step = xstep5;
978 x4min = min5;
979 x4max = max5;
980 x5step = xstep6;
981 x5min = min6;
982 x5max = max6;
983 x6step = xstep4;
984 x6min = min4;
985 x6max = max4;
986 break;
987 default:
988 Error_Box ("Unhandled sequence of modulated axes.\n");
989 free (ftmp);
990 fclose (mapin);
991 return;
992 break;
993 }
994
995 // skip to end of first record - record length is max(nrec,34)*4
996
997 nskip = 34;
998 if (nrec > nskip)
999 nskip = nrec;
1000 nskip -= 34; // header values already read
1001 i = fread (ftmp, 4, nskip, mapin);
1002
1003 Map_Info.xlim[0] = xmin;
1004 Map_Info.xlim[1] = xmax;
1005 Map_Info.ylim[0] = ymin;
1006 Map_Info.ylim[1] = ymax;
1007 Map_Info.zlim[0] = zmin;
1008 Map_Info.zlim[1] = zmax;
1009
1010 fullcell = 1;
1011 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
1012 // not full cell along x
1013 new_mapstep_a =
1014 (int) ((1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0])) * (float) mapstep_a +
1015 0.5f);
1016 astart = (int) (Map_Info.xlim[0] * new_mapstep_a + 0.5f);
1017 aend = (int) (Map_Info.xlim[1] * new_mapstep_a + 0.5f);
1018 fullcell = 0;
1019 if (astart < 0 && aend - astart != mapstep_a)
1020 aend++;
1021 if (astart > 0 && aend - astart != mapstep_a)
1022 aend--;
1023 if (aend - astart != mapstep_a)
1024 Error_Box ("Error recalculating map a limits for full cell");
1025 } else {
1026 new_mapstep_a = mapstep_a;
1027 astart = 0;
1028 aend = mapstep_a;
1029 }
1030 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
1031 // not full cell along y
1032 new_mapstep_b =
1033 (int) ((1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0])) * (float) mapstep_b +
1034 0.5f);
1035 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b + 0.5f);
1036 bend = (int) (Map_Info.ylim[1] * new_mapstep_b + 0.5f);
1037 fullcell = 0;
1038 if (bstart < 0 && bend - bstart != mapstep_b)
1039 bend++;
1040 if (bstart > 0 && bend - bstart != mapstep_b)
1041 bend--;
1042 if (bend - bstart != mapstep_b)
1043 Error_Box ("Error recalculating map b limits for full cell");
1044 } else {
1045 new_mapstep_b = mapstep_b;
1046 bstart = 0;
1047 bend = mapstep_b;
1048 }
1049 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
1050 // not full cell along z
1051 new_mapstep_c =
1052 (int) ((1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0])) * (float) mapstep_c +
1053 0.5f);
1054 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c + 0.5f);
1055 cend = (int) (Map_Info.zlim[1] * new_mapstep_c + 0.5f);
1056 fullcell = 0;
1057 if (cstart < 0 && cend - cstart != mapstep_c)
1058 cend++;
1059 if (cstart > 0 && cend - cstart != mapstep_c)
1060 cend--;
1061 if (cend - cstart != mapstep_c)
1062 Error_Box ("Error recalculating map c limits for full cell");
1063 } else {
1064 new_mapstep_c = mapstep_c;
1065 cstart = 0;
1066 cend = mapstep_c;
1067 }
1068 // if (FourierPt) free(FourierPt);
1069 FourierPt =
1070 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c * sizeof (float));
1071 if (FourierPt == NULL) {
1072 ireslt = 0;
1073 Error_Box ("ERROR -- Unable to allocate space for map.");
1074 free (ftmp);
1075 fclose (mapin);
1076 return;
1077 }
1078
1079 switch (modaxes) {
1080 case 0:
1081 break;
1082 case 400:
1083 if (x4Val > 0.0f && x4Val <= x4max) {
1084 nskip = (int) ((x4Val - x4min) / x4step + .5);
1085 for (i = 0; i < nskip; i++) {
1086 for (j = 0; j < step3; j++) {
1087 ireslt = fread (ftmp, 4, nrec, mapin);
1088 if (ireslt <= 0) {
1089 // fprintf (stderr, "error parsing superspace map\n");
1090 Error_Box ("Error seeking section in superspace map.");
1091 free (ftmp);
1092 fclose (mapin);
1093 return;
1094 }
1095 }
1096 }
1097 }
1098 break;
1099 case 450:
1100 if (x5Val > 0.0f && x5Val <= x5max) {
1101 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1102 for (i = 0; i < nskip5; i++) {
1103 for (j = 0; j < step4; j++) {
1104 for (k = 0; k < step3; k++) {
1105 ireslt = fread (ftmp, 4, nrec, mapin);
1106 if (ireslt <= 0) {
1107 // fprintf (stderr, "error parsing superspace map\n");
1108 Error_Box ("Error seeking section in superspace map.");
1109 free (ftmp);
1110 fclose (mapin);
1111 return;
1112 }
1113 }
1114 }
1115 }
1116 }
1117 if (x4Val > 0.0f && x4Val <= x4max) {
1118 nskip = (int) ((x4Val - x4min) / x4step + .5);
1119 for (i = 0; i < nskip; i++) {
1120 for (j = 0; j < step3; j++) {
1121 ireslt = fread (ftmp, 4, nrec, mapin);
1122 if (ireslt <= 0) {
1123 // fprintf (stderr, "error parsing superspace map\n");
1124 Error_Box ("Error seeking section in superspace map.");
1125 free (ftmp);
1126 fclose (mapin);
1127 return;
1128 }
1129 }
1130 }
1131 }
1132 break;
1133 case 540:
1134 if (x4Val > 0.0f && x4Val <= x4max) {
1135 nskip = (int) ((x4Val - x4min) / x4step + .5);
1136 for (i = 0; i < nskip; i++) {
1137 for (j = 0; j < step5; j++) {
1138 for (k = 0; k < step3; k++) {
1139 ireslt = fread (ftmp, 4, nrec, mapin);
1140 if (ireslt <= 0) {
1141 // fprintf (stderr, "error parsing superspace map\n");
1142 Error_Box ("Error seeking section in superspace map.");
1143 free (ftmp);
1144 fclose (mapin);
1145 return;
1146 }
1147 }
1148 }
1149 }
1150 }
1151 if (x5Val > 0.0f && x5Val <= x5max) {
1152 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1153 for (i = 0; i < nskip5; i++) {
1154 for (j = 0; j < step3; j++) {
1155 ireslt = fread (ftmp, 4, nrec, mapin);
1156 if (ireslt <= 0) {
1157 // fprintf (stderr, "error parsing superspace map\n");
1158 Error_Box ("Error seeking section in superspace map.");
1159 free (ftmp);
1160 fclose (mapin);
1161 return;
1162 }
1163 }
1164 }
1165 }
1166 break;
1167 case 456:
1168 if (x6Val > 0.0f && x6Val <= x6max) {
1169 nskip6 = (int) ((x6Val - x6min) / x6step + .5);
1170 for (i = 0; i < nskip6; i++) {
1171 for (j = 0; j < step5; j++) {
1172 for (k = 0; k < step4; k++) {
1173 for (l = 0; l < step3; l++) {
1174 ireslt = fread (ftmp, 4, nrec, mapin);
1175 if (ireslt <= 0) {
1176 // fprintf (stderr, "error parsing superspace map\n");
1177 Error_Box ("Error seeking section in superspace map.");
1178 free (ftmp);
1179 fclose (mapin);
1180 return;
1181 }
1182 }
1183 }
1184 }
1185 }
1186 }
1187 if (x5Val > 0.0f && x5Val <= x5max) {
1188 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1189 for (i = 0; i < nskip5; i++) {
1190 for (j = 0; j < step4; j++) {
1191 for (k = 0; k < step3; k++) {
1192 ireslt = fread (ftmp, 4, nrec, mapin);
1193 if (ireslt <= 0) {
1194 // fprintf (stderr, "error parsing superspace map\n");
1195 Error_Box ("Error seeking section in superspace map.");
1196 free (ftmp);
1197 fclose (mapin);
1198 return;
1199 }
1200 }
1201 }
1202 }
1203 }
1204 if (x4Val > 0.0f && x4Val <= x4max) {
1205 nskip = (int) ((x4Val - x4min) / x4step + .5);
1206 for (i = 0; i < nskip; i++) {
1207 for (j = 0; j < step3; j++) {
1208 ireslt = fread (ftmp, 4, nrec, mapin);
1209 if (ireslt <= 0) {
1210 // fprintf (stderr, "error parsing superspace map\n");
1211 Error_Box ("Error seeking section in superspace map.");
1212 free (ftmp);
1213 fclose (mapin);
1214 return;
1215 }
1216 }
1217 }
1218 }
1219 break;
1220 case 465:
1221 if (x5Val > 0.0f && x5Val <= x5max) {
1222 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1223 for (i = 0; i < nskip5; i++) {
1224 for (j = 0; j < step6; j++) {
1225 for (k = 0; k < step4; k++) {
1226 for (l = 0; l < step3; l++) {
1227 ireslt = fread (ftmp, 4, nrec, mapin);
1228 if (ireslt <= 0) {
1229 // fprintf (stderr, "error parsing superspace map\n");
1230 Error_Box ("Error seeking section in superspace map.");
1231 free (ftmp);
1232 fclose (mapin);
1233 return;
1234 }
1235 }
1236 }
1237 }
1238 }
1239 }
1240 if (x6Val > 0.0f && x6Val <= x6max) {
1241 nskip6 = (int) ((x6Val - x6min) / x6step + .5);
1242 for (i = 0; i < nskip6; i++) {
1243 for (j = 0; j < step4; j++) {
1244 for (k = 0; k < step3; k++) {
1245 ireslt = fread (ftmp, 4, nrec, mapin);
1246 if (ireslt <= 0) {
1247 // fprintf (stderr, "error parsing superspace map\n");
1248 Error_Box ("Error seeking section in superspace map.");
1249 free (ftmp);
1250 fclose (mapin);
1251 return;
1252 }
1253 }
1254 }
1255 }
1256 }
1257 if (x4Val > 0.0f && x4Val <= x4max) {
1258 nskip = (int) ((x4Val - x4min) / x4step + .5);
1259 for (i = 0; i < nskip; i++) {
1260 for (j = 0; j < step3; j++) {
1261 ireslt = fread (ftmp, 4, nrec, mapin);
1262 if (ireslt <= 0) {
1263 // fprintf (stderr, "error parsing superspace map\n");
1264 Error_Box ("Error seeking section in superspace map.");
1265 free (ftmp);
1266 fclose (mapin);
1267 return;
1268 }
1269 }
1270 }
1271 }
1272 break;
1273 case 546:
1274 if (x6Val > 0.0f && x6Val <= x6max) {
1275 nskip6 = (int) ((x6Val - x6min) / x6step + .5);
1276 for (i = 0; i < nskip6; i++) {
1277 for (j = 0; j < step4; j++) {
1278 for (k = 0; k < step5; k++) {
1279 for (l = 0; l < step3; l++) {
1280 ireslt = fread (ftmp, 4, nrec, mapin);
1281 if (ireslt <= 0) {
1282 // fprintf (stderr, "error parsing superspace map\n");
1283 Error_Box ("Error seeking section in superspace map.");
1284 free (ftmp);
1285 fclose (mapin);
1286 return;
1287 }
1288 }
1289 }
1290 }
1291 }
1292 }
1293 if (x4Val > 0.0f && x4Val <= x4max) {
1294 nskip = (int) ((x4Val - x4min) / x4step + .5);
1295 for (i = 0; i < nskip; i++) {
1296 for (j = 0; j < step5; j++) {
1297 for (k = 0; k < step3; k++) {
1298 ireslt = fread (ftmp, 4, nrec, mapin);
1299 if (ireslt <= 0) {
1300 // fprintf (stderr, "error parsing superspace map\n");
1301 Error_Box ("Error seeking section in superspace map.");
1302 free (ftmp);
1303 fclose (mapin);
1304 return;
1305 }
1306 }
1307 }
1308 }
1309 }
1310 if (x5Val > 0.0f && x5Val <= x5max) {
1311 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1312 for (i = 0; i < nskip5; i++) {
1313 for (j = 0; j < step3; j++) {
1314 ireslt = fread (ftmp, 4, nrec, mapin);
1315 if (ireslt <= 0) {
1316 // fprintf (stderr, "error parsing superspace map\n");
1317 Error_Box ("Error seeking section in superspace map.");
1318 free (ftmp);
1319 fclose (mapin);
1320 return;
1321 }
1322 }
1323 }
1324 }
1325 break;
1326 case 564:
1327 if (x4Val > 0.0f && x4Val <= x4max) {
1328 nskip = (int) ((x4Val - x4min) / x4step + .5);
1329 for (i = 0; i < nskip; i++) {
1330 for (j = 0; j < step6; j++) {
1331 for (k = 0; k < step5; k++) {
1332 for (l = 0; l < step3; l++) {
1333 ireslt = fread (ftmp, 4, nrec, mapin);
1334 if (ireslt <= 0) {
1335 // fprintf (stderr, "error parsing superspace map\n");
1336 Error_Box ("Error seeking section in superspace map.");
1337 free (ftmp);
1338 fclose (mapin);
1339 return;
1340 }
1341 }
1342 }
1343 }
1344 }
1345 }
1346 if (x6Val > 0.0f && x6Val <= x6max) {
1347 nskip6 = (int) ((x6Val - x6min) / x6step + .5);
1348 for (i = 0; i < nskip6; i++) {
1349 for (j = 0; j < step5; j++) {
1350 for (k = 0; k < step3; k++) {
1351 ireslt = fread (ftmp, 4, nrec, mapin);
1352 if (ireslt <= 0) {
1353 // fprintf (stderr, "error parsing superspace map\n");
1354 Error_Box ("Error seeking section in superspace map.");
1355 free (ftmp);
1356 fclose (mapin);
1357 return;
1358 }
1359 }
1360 }
1361 }
1362 }
1363 if (x5Val > 0.0f && x5Val <= x5max) {
1364 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1365 for (i = 0; i < nskip5; i++) {
1366 for (j = 0; j < step3; j++) {
1367 ireslt = fread (ftmp, 4, nrec, mapin);
1368 if (ireslt <= 0) {
1369 // fprintf (stderr, "error parsing superspace map\n");
1370 Error_Box ("Error seeking section in superspace map.");
1371 free (ftmp);
1372 fclose (mapin);
1373 return;
1374 }
1375 }
1376 }
1377 }
1378 break;
1379 case 654:
1380 if (x4Val > 0.0f && x4Val <= x4max) {
1381 nskip = (int) ((x4Val - x4min) / x4step + .5);
1382 for (i = 0; i < nskip; i++) {
1383 for (j = 0; j < step5; j++) {
1384 for (k = 0; k < step6; k++) {
1385 for (l = 0; l < step3; l++) {
1386 ireslt = fread (ftmp, 4, nrec, mapin);
1387 if (ireslt <= 0) {
1388 // fprintf (stderr, "error parsing superspace map\n");
1389 Error_Box ("Error seeking section in superspace map.");
1390 free (ftmp);
1391 fclose (mapin);
1392 return;
1393 }
1394 }
1395 }
1396 }
1397 }
1398 }
1399 if (x5Val > 0.0f && x5Val <= x5max) {
1400 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1401 for (i = 0; i < nskip5; i++) {
1402 for (j = 0; j < step6; j++) {
1403 for (k = 0; k < step3; k++) {
1404 ireslt = fread (ftmp, 4, nrec, mapin);
1405 if (ireslt <= 0) {
1406 // fprintf (stderr, "error parsing superspace map\n");
1407 Error_Box ("Error seeking section in superspace map.");
1408 free (ftmp);
1409 fclose (mapin);
1410 return;
1411 }
1412 }
1413 }
1414 }
1415 }
1416 if (x6Val > 0.0f && x6Val <= x6max) {
1417 nskip6 = (int) ((x6Val - x6min) / x6step + .5);
1418 for (i = 0; i < nskip6; i++) {
1419 for (j = 0; j < step3; j++) {
1420 ireslt = fread (ftmp, 4, nrec, mapin);
1421 if (ireslt <= 0) {
1422 // fprintf (stderr, "error parsing superspace map\n");
1423 Error_Box ("Error seeking section in superspace map.");
1424 free (ftmp);
1425 fclose (mapin);
1426 return;
1427 }
1428 }
1429 }
1430 }
1431 break;
1432 case 645:
1433 if (x5Val > 0.0f && x5Val <= x5max) {
1434 nskip5 = (int) ((x5Val - x5min) / x5step + .5);
1435 for (i = 0; i < nskip5; i++) {
1436 for (j = 0; j < step4; j++) {
1437 for (k = 0; k < step6; k++) {
1438 for (l = 0; l < step3; l++) {
1439 ireslt = fread (ftmp, 4, nrec, mapin);
1440 if (ireslt <= 0) {
1441 // fprintf (stderr, "error parsing superspace map\n");
1442 Error_Box ("Error seeking section in superspace map.");
1443 free (ftmp);
1444 fclose (mapin);
1445 return;
1446 }
1447 }
1448 }
1449 }
1450 }
1451 }
1452 if (x4Val > 0.0f && x4Val <= x4max) {
1453 nskip = (int) ((x4Val - x4min) / x4step + .5);
1454 for (i = 0; i < nskip; i++) {
1455 for (j = 0; j < step6; j++) {
1456 for (k = 0; k < step3; k++) {
1457 ireslt = fread (ftmp, 4, nrec, mapin);
1458 if (ireslt <= 0) {
1459 // fprintf (stderr, "error parsing superspace map\n");
1460 Error_Box ("Error seeking section in superspace map.");
1461 free (ftmp);
1462 fclose (mapin);
1463 return;
1464 }
1465 }
1466 }
1467 }
1468 }
1469 if (x6Val > 0.0f && x6Val <= x6max) {
1470 nskip6 = (int) ((x6Val - x6min) / x6step + .5);
1471 for (i = 0; i < nskip6; i++) {
1472 for (j = 0; j < step3; j++) {
1473 ireslt = fread (ftmp, 4, nrec, mapin);
1474 if (ireslt <= 0) {
1475 // fprintf (stderr, "error parsing superspace map\n");
1476 Error_Box ("Error seeking section in superspace map.");
1477 free (ftmp);
1478 fclose (mapin);
1479 return;
1480 }
1481 }
1482 }
1483 }
1484 break;
1485
1486 default:
1487 Error_Box
1488 ("ERROR -- Unable to select desired map section (unsupported axis numbering");
1489 free (ftmp);
1490 fclose (mapin);
1491 return;
1492 break;
1493 }
1494
1495 // read the Rho values
1496 rhomin = 100000.0f;
1497 rhomax = 0.000001f;
1498 if (fullcell == 1) { //need to skip redundant information at end
1499 new_mapstep_a--;
1500 new_mapstep_b--;
1501 new_mapstep_c--;
1502 switch (axes) {
1503 case 123: // X Y Z order
1504 for (i = cstart; i < cend; i++) { //c
1505 for (j = bstart; j < bend - 1; j++) { //b
1506 for (k = astart; k < aend - 1; k++) {
1507 ireslt = fread (ftmp, 4, 1, mapin);
1508 ijk = k * new_mapstep_c * new_mapstep_b + j * new_mapstep_c + i;
1509 FourierPt[ijk] = ftmp[0];
1510 if (FourierPt[ijk] < rhomin)
1511 rhomin = FourierPt[ijk];
1512 if (FourierPt[ijk] > rhomax)
1513 rhomax = FourierPt[ijk];
1514 }
1515 fread (ftmp, 4, 1, mapin);
1516 }
1517 for (k = astart; k < aend; k++)
1518 fread (ftmp, 4, 1, mapin);
1519 }
1520 break;
1521 case 312: // ZXY
1522 for (i = bstart; i < bend - 1; i++) {
1523 for (j = astart; j < aend - 1; j++) {
1524 for (k = cstart; k < cend - 1; k++) {
1525 fread (ftmp, 4, 1, mapin);
1526 ijk = j * new_mapstep_c * new_mapstep_b + i * new_mapstep_c + k;
1527 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1528 FourierPt[ijk] = ftmp[0];
1529 if (FourierPt[ijk] < rhomin)
1530 rhomin = FourierPt[ijk];
1531 if (FourierPt[ijk] > rhomax)
1532 rhomax = FourierPt[ijk];
1533 }
1534 fread (ftmp, 4, 1, mapin);
1535 }
1536 for (k = cstart; k < cend; k++)
1537 fread (ftmp, 4, 1, mapin);
1538 }
1539 break;
1540 case 321: //ZYX
1541 for (i = astart; i < aend - 1; i++) {
1542 for (j = cstart; j < cend - 1; j++) {
1543 for (k = bstart; k < bend - 1; k++) {
1544 fread (ftmp, 4, 1, mapin);
1545 ijk = i * new_mapstep_c * new_mapstep_b + j * new_mapstep_b + k;
1546 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1547 FourierPt[ijk] = ftmp[0];
1548 if (FourierPt[ijk] < rhomin)
1549 rhomin = FourierPt[ijk];
1550 if (FourierPt[ijk] > rhomax)
1551 rhomax = FourierPt[ijk];
1552 }
1553 fread (ftmp, 4, 1, mapin);
1554 }
1555 for (k = bstart; k < bend; k++)
1556 fread (ftmp, 4, 1, mapin);
1557 }
1558 break;
1559 case 213: //YXZ
1560 for (i = cstart; i < cend - 1; i++) {
1561 for (j = astart; j < aend - 1; j++) {
1562 for (k = bstart; k < bend - 1; k++) {
1563 fread (ftmp, 4, 1, mapin);
1564 ijk = j * new_mapstep_c * new_mapstep_b + k * new_mapstep_c + i;
1565 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1566 FourierPt[ijk] = ftmp[0];
1567 if (FourierPt[ijk] < rhomin)
1568 rhomin = FourierPt[ijk];
1569 if (FourierPt[ijk] > rhomax)
1570 rhomax = FourierPt[ijk];
1571 }
1572 fread (ftmp, 4, 1, mapin);
1573 }
1574 for (k = bstart; k < bend; k++)
1575 fread (ftmp, 4, 1, mapin);
1576 }
1577 break;
1578 case 231: //YZX
1579 for (i = astart; i < aend - 1; i++) {
1580 for (j = cstart; j < cend - 1; j++) {
1581 for (k = bstart; k < bend - 1; k++) {
1582 fread (ftmp, 4, 1, mapin);
1583 ijk = i * new_mapstep_c * new_mapstep_b + k * new_mapstep_c + j;
1584 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1585 FourierPt[ijk] = ftmp[0];
1586 if (FourierPt[ijk] < rhomin)
1587 rhomin = FourierPt[ijk];
1588 if (FourierPt[ijk] > rhomax)
1589 rhomax = FourierPt[ijk];
1590 }
1591 fread (ftmp, 4, 1, mapin);
1592 }
1593 for (k = bstart; k < bend; k++)
1594 fread (ftmp, 4, 1, mapin);
1595 }
1596 break;
1597 case 132: //XZY
1598 for (i = bstart; i < bend - 1; i++) {
1599 for (j = cstart; j < cend - 1; j++) {
1600 for (k = astart; k < aend - 1; k++) {
1601 fread (ftmp, 4, 1, mapin);
1602 ijk = k * new_mapstep_c * new_mapstep_b + i * new_mapstep_c + j;
1603 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1604 FourierPt[ijk] = ftmp[0];
1605 if (FourierPt[ijk] < rhomin)
1606 rhomin = FourierPt[ijk];
1607 if (FourierPt[ijk] > rhomax)
1608 rhomax = FourierPt[ijk];
1609 }
1610 fread (ftmp, 4, 1, mapin);
1611 }
1612 for (k = astart; k < aend; k++)
1613 fread (ftmp, 4, 1, mapin);
1614 }
1615 break;
1616 default:
1617 Error_Box ("Primary axes must be some permutation of X,Y,Z .\n");
1618 free (ftmp);
1619 fclose (mapin);
1620 return;
1621 break;
1622 }
1623 // new_mapstep_a-=1;
1624 // new_mapstep_b-=1;
1625 // new_mapstep_c-=1;
1626 } else {
1627 switch (axes) {
1628 case 123: // X Y Z order
1629 for (i = cstart; i < cend; i++) { //c
1630 for (j = bstart; j < bend; j++) { //b
1631 for (k = astart; k < aend; k++) {
1632 ireslt = fread (ftmp, 4, 1, mapin);
1633 ijk = k * new_mapstep_c * new_mapstep_b + j * new_mapstep_c + i;
1634 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1635 FourierPt[ijk] = ftmp[0];
1636 if (FourierPt[ijk] < rhomin)
1637 rhomin = FourierPt[ijk];
1638 if (FourierPt[ijk] > rhomax)
1639 rhomax = FourierPt[ijk];
1640 }
1641 }
1642 }
1643 break;
1644 case 312: // ZXY
1645 for (i = bstart; i < bend; i++) {
1646 for (j = astart; j < aend; j++) {
1647 for (k = cstart; k < cend; k++) {
1648 fread (ftmp, 4, 1, mapin);
1649 ijk = j * new_mapstep_c * new_mapstep_b + i * new_mapstep_c + k;
1650 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1651 FourierPt[ijk] = ftmp[0];
1652 if (FourierPt[ijk] < rhomin)
1653 rhomin = FourierPt[ijk];
1654 if (FourierPt[ijk] > rhomax)
1655 rhomax = FourierPt[ijk];
1656 }
1657 }
1658 }
1659 break;
1660 case 321: //ZYX
1661 for (i = astart; i < aend; i++) {
1662 for (j = cstart; j < cend; j++) {
1663 for (k = bstart; k < bend; k++) {
1664 fread (ftmp, 4, 1, mapin);
1665 ijk = i * new_mapstep_c * new_mapstep_b + j * new_mapstep_b + k;
1666 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1667 FourierPt[ijk] = ftmp[0];
1668 if (FourierPt[ijk] < rhomin)
1669 rhomin = FourierPt[ijk];
1670 if (FourierPt[ijk] > rhomax)
1671 rhomax = FourierPt[ijk];
1672 }
1673 }
1674 }
1675 break;
1676 case 213: //YXZ
1677 for (i = cstart; i < cend; i++) {
1678 for (j = astart; j < aend; j++) {
1679 for (k = bstart; k < bend; k++) {
1680 fread (ftmp, 4, 1, mapin);
1681 ijk = j * new_mapstep_c * new_mapstep_b + k * new_mapstep_c + i;
1682 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1683 FourierPt[ijk] = ftmp[0];
1684 if (FourierPt[ijk] < rhomin)
1685 rhomin = FourierPt[ijk];
1686 if (FourierPt[ijk] > rhomax)
1687 rhomax = FourierPt[ijk];
1688 }
1689 }
1690 }
1691 break;
1692 case 231: //YZX
1693 for (i = astart; i < aend; i++) {
1694 for (j = cstart; j < cend; j++) {
1695 for (k = bstart; k < bend; k++) {
1696 fread (ftmp, 4, 1, mapin);
1697 ijk = i * new_mapstep_c * new_mapstep_b + k * new_mapstep_c + j;
1698 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1699 FourierPt[ijk] = ftmp[0];
1700 if (FourierPt[ijk] < rhomin)
1701 rhomin = FourierPt[ijk];
1702 if (FourierPt[ijk] > rhomax)
1703 rhomax = FourierPt[ijk];
1704 }
1705 }
1706 }
1707 break;
1708 case 132: //XZY
1709 for (i = bstart; i < bend; i++) {
1710 for (j = cstart; j < cend; j++) {
1711 for (k = astart; k < aend; k++) {
1712 fread (ftmp, 4, 1, mapin);
1713 ijk = k * new_mapstep_c * new_mapstep_b + i * new_mapstep_c + j;
1714 // if (wrong_end) dens[k] = end_flip_real(dens[k]);
1715 FourierPt[ijk] = ftmp[0];
1716 if (FourierPt[ijk] < rhomin)
1717 rhomin = FourierPt[ijk];
1718 if (FourierPt[ijk] > rhomax)
1719 rhomax = FourierPt[ijk];
1720 }
1721 }
1722 }
1723 break;
1724 default:
1725 Error_Box ("Primary axes must be some permutation of X,Y,Z .\n");
1726 free (ftmp);
1727 fclose (mapin);
1728 return;
1729 break;
1730 }
1731 }
1732
1733 mapstep_a = new_mapstep_a;
1734 mapstep_b = new_mapstep_b;
1735 mapstep_c = new_mapstep_c;
1736 fclose (mapin);
1737 free (ftmp);
1738
1739 ReadFourMap = 1;
1740 sprintf (Map_Info.title, "%s map with axes %s %s %s %s %s %s\n", maptypes[maptype],
1741 theaxes[axis[0]], theaxes[axis[1]], theaxes[axis[2]],
1742 theaxes[axis[3]], theaxes[axis[4]], theaxes[axis[5]]);
1743 Map_Info.rhomn = rhomin;
1744 Map_Info.rhomx = rhomax;
1745 Map_Info.map_int[0] = mapstep_a - 1;
1746 Map_Info.map_int[1] = mapstep_b - 1;
1747 Map_Info.map_int[2] = mapstep_c - 1;
1748 Map_Info.lat_con[0] = drvui->lat_con[0];
1749 Map_Info.lat_con[1] = drvui->lat_con[1];
1750 Map_Info.lat_con[2] = drvui->lat_con[2];
1751 Map_Info.lat_con[3] = drvui->lat_con[3];
1752 Map_Info.lat_con[4] = drvui->lat_con[4];
1753 Map_Info.lat_con[5] = drvui->lat_con[5];
1754 Map_Info.info_valid = 1;
1755
1756 Map_Info.x4lim[0] = x4min;
1757 Map_Info.x4lim[1] = x4max;
1758 Map_Info.x5lim[0] = x5min;
1759 Map_Info.x5lim[1] = x5max;
1760 Map_Info.x6lim[0] = x6min;
1761 Map_Info.x6lim[1] = x6max;
1762 #if 1
1763 xMin = Map_Info.xlim[0] = xmin;
1764 xMax = Map_Info.xlim[1] = xmax;
1765 yMin = Map_Info.ylim[0] = ymin;
1766 yMax = Map_Info.ylim[1] = ymax;
1767 zMin = Map_Info.zlim[0] = zmin;
1768 zMax = Map_Info.zlim[1] = zmax;
1769 #endif
1770
1771 if (!Quick) {
1772 fprintf (drvui->flout, "Reading %s map file %s\n"
1773 "axis sequence = %s %s %s %s %s %s\n map steps: a=%d b=%d c=%d\n "
1774 "Range of rho values: %f to %f\n", maptypes[maptype], infile,
1775 theaxes[axis[0]], theaxes[axis[1]], theaxes[axis[2]],
1776 theaxes[axis[3]], theaxes[axis[4]], theaxes[axis[5]],
1777 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax);
1778 if (axis[3] != 0)
1779 fprintf (drvui->flout, "X4 section %2f to %2f stepsize %2f\n",
1780 x4min, x4max, x4step);
1781 if (axis[4] != 0)
1782 fprintf (drvui->flout, "X5 section %2f to %2f stepsize %2f\n",
1783 x5min, x5max, x5step);
1784 if (axis[5] != 0)
1785 fprintf (drvui->flout, "X6 section %2f to %2f stepsize %2f\n",
1786 x6min, x6max, x6step);
1787 }
1788
1789 }
1790
1791 void
read_m80(char * infile,int Quick)1792 read_m80 (char *infile, int Quick)
1793 {
1794 // read JANA .m80 Fo/Fc file and compute density
1795
1796 FILE *mapin;
1797
1798 float *fo, *fc, *f2, *ftmp1, *ftmp2, *ftmp3;
1799
1800 short *ih, *ik, *il, *itmp1, *itmp2, *itmp3;
1801
1802 float *Cklx;
1803
1804 float *Dklx;
1805
1806 float *Elxy;
1807
1808 float *Flxy;
1809
1810 double x, y, z;
1811
1812 char string[200];
1813
1814 int nval = 200;
1815
1816 int nr;
1817
1818 int i, j, k, l;
1819
1820 int ijk = 0;
1821
1822 float rhomin, rhomax;
1823
1824 float factor;
1825
1826 short ih0, ik0, il0;
1827
1828 short ih1, ik1, il1;
1829
1830 short im, in, io, is;
1831
1832 float fo0, fc0, fc1, f20, f21;
1833
1834 int kmin, kmax, nk;
1835
1836 int lmin, lmax, nl;
1837
1838 float phase;
1839
1840 int skip;
1841
1842 int progress;
1843
1844 float foj;
1845
1846 int warnings;
1847
1848 int dimension;
1849
1850 if (FourierPt != NULL) {
1851 return;
1852 }
1853
1854 if ((mapin = fopen (infile, "r")) == NULL) {
1855 sprintf (string, "Cannot open Fo/Fc (JANA m80) file %s", infile);
1856 Error_Box (string);
1857 return;
1858 }
1859
1860 fo = (float *) zalloc (200 * sizeof (float));
1861 if (fo == NULL) {
1862 sprintf (string, "Unable to obtain storage for Fo data");
1863 Error_Box (string);
1864 fclose (mapin);
1865 return;
1866 }
1867 fc = (float *) zalloc (200 * sizeof (float));
1868 if (fc == NULL) {
1869 sprintf (string, "Unable to obtain storage for Fc data");
1870 Error_Box (string);
1871 free (fo);
1872 fclose (mapin);
1873 return;
1874 }
1875 f2 = (float *) zalloc (200 * sizeof (float));
1876 if (f2 == NULL) {
1877 sprintf (string, "Unable to obtain storage for B or phase data");
1878 Error_Box (string);
1879 free (fo);
1880 free (fc);
1881 fclose (mapin);
1882 return;
1883 }
1884 ih = (short *) zalloc (200 * sizeof (short));
1885 if (ih == NULL) {
1886 sprintf (string, "Unable to obtain storage for h data");
1887 Error_Box (string);
1888 free (fo);
1889 free (fc);
1890 free (f2);
1891 fclose (mapin);
1892 return;
1893 }
1894 ik = (short *) zalloc (200 * sizeof (short));
1895 if (ik == NULL) {
1896 sprintf (string, "Unable to obtain storage for k data");
1897 Error_Box (string);
1898 free (fo);
1899 free (fc);
1900 free (f2);
1901 free (ih);
1902 fclose (mapin);
1903 return;
1904 }
1905 il = (short *) zalloc (200 * sizeof (short));
1906 if (il == NULL) {
1907 sprintf (string, "Unable to obtain storage for l data");
1908 Error_Box (string);
1909 free (fo);
1910 free (fc);
1911 free (f2);
1912 free (ih);
1913 free (ik);
1914 fclose (mapin);
1915 return;
1916 }
1917
1918 if (fgets (string, 199, mapin)) {
1919 i = sscanf (string, "%hd %hd %hd %hd %hd %hd %hd %f", &ih0, &ik0, &il0, &im, &in,
1920 &io, &is, &fo0);
1921 if (i > 4)
1922 Map_Info.info_valid = 1;
1923 dimension = i - 2;
1924 } else {
1925 sprintf (string, "Error reading M80 file");
1926 Error_Box (string);
1927 free (fo);
1928 free (fc);
1929 free (f2);
1930 free (ih);
1931 free (ik);
1932 free (il);
1933 fclose (mapin);
1934 return;
1935 }
1936
1937 if (dimension < 3) {
1938 sprintf (string, "Error reading M80 file");
1939 Error_Box (string);
1940 free (fo);
1941 free (fc);
1942 free (f2);
1943 free (ih);
1944 free (ik);
1945 free (il);
1946 fclose (mapin);
1947 return;
1948 }
1949
1950 nr = 0;
1951 warnings = 0;
1952
1953 do {
1954 skip = 0;
1955 if (dimension > 3) {
1956 switch (dimension) {
1957 case 4:
1958 im = 0;
1959 i = sscanf (string, "%hd %hd %hd %hd %*d %f %*f %*f %f %f", &ih0, &ik0,
1960 &il0, &im, &fo0, &fc0, &f20);
1961 if (im != 0)
1962 skip = 1;
1963 break;
1964 case 5:
1965 im = in = 0;
1966 i = sscanf (string, "%hd %hd %hd %hd %hd %*d %f %*f %*f %f %f", &ih0,
1967 &ik0, &il0, &im, &in, &fo0, &fc0, &f20);
1968 if (im != 0 || in != 0)
1969 skip = 1;
1970 break;
1971 case 6:
1972 im = in = io = 0;
1973 i = sscanf (string, "%hd %hd %hd %hd %hd %hd %*d %f %*f %*f %f %f",
1974 &ih0, &ik0, &il0, &im, &in, &io, &fo0, &fc0, &f20);
1975 if (im != 0 || in != 0 || io != 0)
1976 skip = 1;
1977 break;
1978 }
1979 } else {
1980 i = sscanf (string, "%hd %hd %hd %*d %f %*f %*f %f %f", &ih0, &ik0, &il0,
1981 &fo0, &fc0, &f20);
1982 }
1983 if (i < dimension + 3)
1984 break;
1985
1986 if (skip == 1) { // satellite reflection of modulated structure
1987 fgets (string, 199, mapin);
1988 continue;
1989 }
1990 // change all data to 'phase' type
1991 fc1 = (float) sqrt (fc0 * fc0 + f20 * f20);
1992 if (fc1 == 0.) // happens with data generated by superflip
1993 f20 = 0.;
1994 else
1995 f20 = (float) atan2 (f20 / fc1, fc0 / fc1);
1996 fc0 = fc1;
1997
1998 // apply spacegroup symmetry (including xyz - catch duplicates in input list)
1999 for (k = 0; k < drvui->ng; ++k) {
2000 skip = 0;
2001
2002 ih1 =
2003 drvui->ss[k][0][0] * ih0 + drvui->ss[k][1][0] * ik0 +
2004 drvui->ss[k][2][0] * il0;
2005 ik1 =
2006 drvui->ss[k][0][1] * ih0 + drvui->ss[k][1][1] * ik0 +
2007 drvui->ss[k][2][1] * il0;
2008 il1 =
2009 drvui->ss[k][0][2] * ih0 + drvui->ss[k][1][2] * ik0 +
2010 drvui->ss[k][2][2] * il0;
2011 phase =
2012 2.0f * (float) PI *(drvui->ts[k][0] * ih0 + drvui->ts[k][1] * ik0 +
2013 drvui->ts[k][2] * il0);
2014 f21 = f20 + phase;
2015 f21 = (float) atan2 (sin (f21), cos (f21)); // get f21 in range -PI to PI
2016
2017 for (j = 0; j < nr; j++) {
2018 if (ih[j] == ih1 && ik[j] == ik1 && il[j] == il1) {
2019 skip = 1;
2020 foj = fo[j];
2021 break;
2022 }
2023 }
2024 if (skip == 0) {
2025 ih[nr] = ih1;
2026 ik[nr] = ik1;
2027 il[nr] = il1;
2028 fo[nr] = fo0;
2029 fc[nr] = fc0;
2030 f2[nr] = f21;
2031 nr++;
2032 } else {
2033 if (fabs (foj - fo0) > 1.e-3) {
2034 if (!Quick && warnings < 100)
2035 fprintf (drvui->flout,
2036 "Warning - different Fo for %d %d %d and %d %d %d\n",
2037 ih0, ik0, il0, ih1, ik1, il1);
2038 warnings++;
2039 }
2040 }
2041
2042 ih1 *= -1; // add Friedel opposite
2043 ik1 *= -1;
2044 il1 *= -1;
2045 f21 *= -1.0;
2046 skip = 0;
2047
2048 for (j = 0; j < nr; j++) {
2049 if (ih[j] == ih1 && ik[j] == ik1 && il[j] == il1) {
2050 skip = 1;
2051 foj = fo[j];
2052 break;
2053 }
2054 }
2055 if (skip == 0) {
2056 ih[nr] = ih1;
2057 ik[nr] = ik1;
2058 il[nr] = il1;
2059 fo[nr] = fo0;
2060 fc[nr] = fc0;
2061 f2[nr] = f21;
2062 nr++;
2063 } else {
2064 if (fabs (foj - fo0) > 1.e-3) {
2065 if (!Quick && warnings < 100)
2066 fprintf (drvui->flout,
2067 "Warning - different Fo for %d %d %d and %d %d %d\n",
2068 ih0, ik0, il0, ih1, ik1, il1);
2069 warnings++;
2070 }
2071 }
2072
2073 if (nr > nval - 10) {
2074 nval += 100;
2075 ftmp1 = (float *) realloc (fo, nval * sizeof (float));
2076 ftmp2 = (float *) realloc (fc, nval * sizeof (float));
2077 ftmp3 = (float *) realloc (f2, nval * sizeof (float));
2078 itmp1 = (short *) realloc (ih, nval * sizeof (short));
2079 itmp2 = (short *) realloc (ik, nval * sizeof (short));
2080 itmp3 = (short *) realloc (il, nval * sizeof (short));
2081 if (!ftmp1 || !ftmp2 || !ftmp3 || !itmp1 || !itmp2 || !itmp3) {
2082 sprintf (string,
2083 "Unable to expand storage for h, k, l, fo, fc or phase data");
2084 free (fo);
2085 free (fc);
2086 free (f2);
2087 free (ih);
2088 free (ik);
2089 free (il);
2090 if (ftmp1)
2091 free (ftmp1);
2092 if (ftmp2)
2093 free (ftmp2);
2094 if (ftmp3)
2095 free (ftmp3);
2096 if (itmp1)
2097 free (itmp1);
2098 if (itmp2)
2099 free (itmp2);
2100 if (itmp3)
2101 free (itmp3);
2102 Error_Box (string);
2103 fclose (mapin);
2104 return;
2105 } else {
2106 fo = ftmp1;
2107 fc = ftmp2;
2108 f2 = ftmp3;
2109 ih = itmp1;
2110 ik = itmp2;
2111 il = itmp3;
2112 }
2113 } // if we need to realloc
2114 } // loop over all symmetry operators
2115 fgets (string, 199, mapin);
2116 } while (i > 0 && !feof (mapin));
2117
2118 fclose (mapin);
2119
2120 if (dimension > 3) {
2121 sprintf (string,
2122 "%d-dimensional data encountered - only main reflections will be used",
2123 dimension);
2124 Error_Box (string);
2125 }
2126 if (warnings > 0) {
2127 sprintf (string,
2128 "Symmetry problem - %d mismatches in symmetry-equivalent Fo values\n",
2129 warnings);
2130 Error_Box (string);
2131 }
2132 factor = 1.0f / (drvui->lat_con[0] * drvui->lat_con[1] * drvui->lat_con[2]
2133 * (float) sqrt (1.0 - cos (drvui->lat_con[3] * PI / 180.0) *
2134 cos (drvui->lat_con[3] * PI / 180.0) *
2135 cos (drvui->lat_con[4] * PI / 180.0) *
2136 cos (drvui->lat_con[4] * PI / 180.0) *
2137 cos (drvui->lat_con[5] * PI / 180.0) *
2138 cos (drvui->lat_con[5] * PI / 180.0)));
2139 kmin = kmax = ik[0];
2140 lmin = lmax = il[0];
2141
2142 for (l = 0; l < nr; l++) {
2143 if (ik[l] < kmin)
2144 kmin = ik[l]; // find range of k and l
2145 if (ik[l] > kmax)
2146 kmax = ik[l];
2147 if (il[l] < lmin)
2148 lmin = il[l];
2149 if (il[l] > lmax)
2150 lmax = il[l];
2151 switch (Map_Info.map_type) { // get map coefficients/V
2152 default:
2153 case 0: // Fo map
2154 break;
2155 case 1: // Fc map
2156 fo[l] = fc[l];
2157 break;
2158 case 2: // Fo - Fc map
2159 fo[l] = fo[l] - fc[l];
2160 break;
2161 case 3: // 2Fo - Fc map
2162 fo[l] = 2.0f * fo[l] - fc[l];
2163 break;
2164 case 4: // Fo2 (Patterson)
2165 fo[l] = fo[l] * fo[l];
2166 f2[l] = 0.0f;
2167 break;
2168 }
2169 fo[l] *= factor; // update coefficients times 1/V
2170 }
2171
2172 nk = kmax - kmin + 1;
2173 nl = lmax - lmin + 1;
2174 Cklx = (float *) zalloc (nk * nl * sizeof (float));
2175 Dklx = (float *) zalloc (nk * nl * sizeof (float));
2176 Elxy = (float *) zalloc (nl * sizeof (float));
2177 Flxy = (float *) zalloc (nl * sizeof (float));
2178
2179 if (!Cklx || !Dklx || !Elxy || !Flxy) {
2180 Error_Box ("Unable to allocate Beevers-Lipson arrays.");
2181 if (Cklx) free (Cklx);
2182 if (Dklx) free (Dklx);
2183 if (Elxy) free (Elxy);
2184 if (Flxy) free (Flxy);
2185 free (fo);
2186 free (fc);
2187 free (f2);
2188 free (ih);
2189 free (ik);
2190 free (il);
2191 return;
2192 }
2193 mapstep_a = (int) (4.0f * drvui->lat_con[0] + 0.5f);
2194 mapstep_b = (int) (4.0f * drvui->lat_con[1] + 0.5f);
2195 mapstep_c = (int) (4.0f * drvui->lat_con[2] + 0.5f);
2196
2197 FourierPt = (float *) zalloc (mapstep_a * mapstep_b * mapstep_c * sizeof (float));
2198 if (FourierPt == NULL) {
2199 Error_Box ("ERROR -- Unable to allocate space for map.");
2200 return;
2201 }
2202
2203 rhomin = 1.0e15f, rhomax = -1.0e15f;
2204
2205 // start the Beevers-Lipson expansion
2206
2207 Progress_Window (-1, "Computing Density", (float) mapstep_a * mapstep_b);
2208
2209 progress = 0;
2210 for (i = 0; i < mapstep_a; i++) {
2211 x = (float) i / (float) mapstep_a;
2212 memset (Cklx, 0, nk * nl * sizeof (float));
2213 memset (Dklx, 0, nk * nl * sizeof (float));
2214 for (l = 0; l < nr; l++) {
2215 ijk = (il[l] - lmin) * nk + ik[l] - kmin;
2216 Cklx[ijk] += fo[l] * (float) cos (2.0 * PI * ih[l] * x - f2[l]);
2217 Dklx[ijk] -= fo[l] * (float) sin (2.0 * PI * ih[l] * x - f2[l]);
2218 }
2219 for (j = 0; j < mapstep_b; j++) {
2220 Progress_Window (0, NULL, (float) ++progress);
2221 y = (float) j / (float) mapstep_b;
2222 memset (Elxy, 0, nl * sizeof (float));
2223 memset (Flxy, 0, nl * sizeof (float));
2224 for (l = 0; l < nl; l++) {
2225 for (k = 0; k < nk; k++) {
2226 ijk = l * nk + k;
2227 Elxy[l] += Cklx[ijk] * (float) cos (2.0 * PI * (k + kmin) * y)
2228 + Dklx[ijk] * (float) sin (2.0 * PI * (k + kmin) * y);
2229 Flxy[l] += Dklx[ijk] * (float) cos (2.0 * PI * (k + kmin) * y)
2230 - Cklx[ijk] * (float) sin (2.0 * PI * (k + kmin) * y);
2231 }
2232 }
2233 for (k = 0; k < mapstep_c; k++) {
2234 z = (float) k / (float) mapstep_c;
2235 ijk = i * mapstep_b * mapstep_c + j * mapstep_c + k;
2236
2237 for (l = 0; l < nl; l++) {
2238 FourierPt[ijk] += Elxy[l] * (float) cos (2.0 * PI * (l + lmin) * z)
2239 + Flxy[l] * (float) sin (2.0 * PI * (l + lmin) * z);
2240 }
2241 if (FourierPt[ijk] < rhomin)
2242 rhomin = FourierPt[ijk];
2243 if (FourierPt[ijk] > rhomax)
2244 rhomax = FourierPt[ijk];
2245 }
2246 }
2247 }
2248 if (!Quick)
2249 fprintf (drvui->fcns, "Reading Fourier map file %s\n "
2250 "map steps: a=%d b=%d c=%d\n "
2251 "Range of rho values: %f to %f, gridsize %d\n", infile,
2252 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax, ijk);
2253
2254
2255 Progress_Window (-2, NULL, 0.0f);
2256 ReadFourMap = 1;
2257 Map_Info.rhomn = rhomin;
2258 Map_Info.rhomx = rhomax;
2259 Map_Info.map_int[0] = mapstep_a;
2260 Map_Info.map_int[1] = mapstep_b;
2261 Map_Info.map_int[2] = mapstep_c;
2262 Map_Info.lat_con[0] = drvui->lat_con[0];
2263 Map_Info.lat_con[1] = drvui->lat_con[1];
2264 Map_Info.lat_con[2] = drvui->lat_con[2];
2265 Map_Info.lat_con[3] = drvui->lat_con[3];
2266 Map_Info.lat_con[4] = drvui->lat_con[4];
2267 Map_Info.lat_con[5] = drvui->lat_con[5];
2268 Map_Info.xlim[0] = 0.0f;
2269 Map_Info.xlim[1] = 1.0f;
2270 Map_Info.ylim[0] = 0.0f;
2271 Map_Info.ylim[1] = 1.0f;
2272 Map_Info.zlim[0] = 0.0f;
2273 Map_Info.zlim[1] = 1.0f;
2274
2275 free (Cklx);
2276 free (Dklx);
2277 free (Elxy);
2278 free (Flxy);
2279 free (fo);
2280 free (fc);
2281 free (f2);
2282 free (ih);
2283 free (ik);
2284 free (il);
2285 }
2286
2287 /* ************************************************************** */
2288 /* ************************************************************** */
2289
2290 void
read_flp(char * infile,int Quick)2291 read_flp (char *infile, int Quick)
2292 {
2293 // read FullProf binary map file
2294
2295 #include "read_flp.h" // define all the struct's
2296
2297 char titulo[81];
2298
2299 char version[6];
2300
2301 struct crystal_cell_type celda;
2302
2303 FILE *mapin;
2304
2305 float rhomin = 1.0e+15f, rhomax = -1.0e+15f;
2306
2307 int ireslt = 0;
2308
2309 char string[200];
2310
2311 int icount;
2312
2313 struct space_group spg;
2314
2315 char temp[16384];
2316
2317 int i, j, k, ijk = 0;
2318
2319 int natoms;
2320
2321 struct FFT_param map_param;
2322
2323 struct Atom_Type atm[200];
2324
2325 float dens[200];
2326
2327 int new_mapstep_a, astart, aend;
2328
2329 int new_mapstep_b, bstart, bend;
2330
2331 int new_mapstep_c, cstart, cend;
2332
2333 int wrong_end = 0;
2334
2335 memset (version, 0, 6);
2336 memset (titulo, 0, 81);
2337 if ((mapin = fopen (infile, "r+b")) == NULL) {
2338 sprintf (string, "Cannot open Fourier map (FulProf) file %s", infile);
2339 Error_Box (string);
2340 return;
2341 }
2342 ireslt = fread (&icount, sizeof (int), 1, mapin); // read record length
2343 if (!ireslt) {
2344 char string[200];
2345
2346 sprintf (string, "Error reading Fourier map file %s", infile);
2347 Error_Box (string);
2348 fclose (mapin);
2349 return;
2350 }
2351 if (icount != 80) {
2352 icount = end_flip (icount);
2353 wrong_end = 1;
2354 }
2355 if (icount != 80) {
2356 Error_Box ("Map does not have the correct format\nJob aborted.");
2357 fclose (mapin);
2358 return;
2359 }
2360 fread (titulo, sizeof (char), icount, mapin); // read title
2361 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2362 trim_string (titulo, 80);
2363 strcpy (Map_Info.title, titulo);
2364
2365 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2366 if (wrong_end)
2367 icount = end_flip (icount);
2368 fread (version, sizeof (char), icount, mapin);
2369 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2370
2371 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2372 if (wrong_end)
2373 icount = end_flip (icount);
2374 fread (&celda, sizeof (char), icount, mapin);
2375 if (wrong_end) {
2376 for (i = 0; i < 3; i++) {
2377 celda.ang[i] = end_flip_real (celda.ang[i]);
2378 celda.cell[i] = end_flip_real (celda.cell[i]);
2379 }
2380 }
2381 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2382 for (i = 0; i < 3; i++) {
2383 Map_Info.lat_con[i] = celda.cell[i];
2384 Map_Info.lat_con[i + 3] = celda.ang[i];
2385 }
2386
2387 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2388 if (wrong_end)
2389 icount = end_flip (icount);
2390 fread (&spg, sizeof (char), icount, mapin);
2391 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2392
2393 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2394 if (wrong_end)
2395 icount = end_flip (icount);
2396 fread (&natoms, sizeof (char), icount, mapin);
2397 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2398
2399 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2400 if (wrong_end)
2401 icount = end_flip (icount);
2402 fread (&atm, sizeof (char), icount, mapin);
2403 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2404
2405 for (i = 0; i < 2; i++) { // skip next 2 records
2406 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2407 if (wrong_end)
2408 icount = end_flip (icount);
2409 fread (&temp, sizeof (char), icount, mapin);
2410 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2411 }
2412
2413 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2414 if (wrong_end)
2415 icount = end_flip (icount);
2416 fread (&map_param, sizeof (char), icount, mapin);
2417 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2418 if (wrong_end) {
2419 for (i = 0; i < 3; i++) {
2420 map_param.ngrid[i] = end_flip (map_param.ngrid[i]);
2421 }
2422 }
2423 Map_Info.map_int[0] = mapstep_a = map_param.ngrid[0];
2424 Map_Info.map_int[1] = mapstep_b = map_param.ngrid[1];
2425 Map_Info.map_int[2] = mapstep_c = map_param.ngrid[2];
2426 for (i = 0; i < 2; i++) {
2427 if (wrong_end) {
2428 map_param.xlim[i] = end_flip_real (map_param.xlim[i]);
2429 map_param.ylim[i] = end_flip_real (map_param.ylim[i]);
2430 map_param.zlim[i] = end_flip_real (map_param.zlim[i]);
2431 }
2432 Map_Info.xlim[i] = map_param.xlim[i];
2433 Map_Info.ylim[i] = map_param.ylim[i];
2434 Map_Info.zlim[i] = map_param.zlim[i];
2435 }
2436 Map_Info.info_valid = 1;
2437
2438 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2439 if (wrong_end)
2440 icount = end_flip (icount);
2441 fread (&temp, sizeof (char), icount, mapin);
2442 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2443
2444 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
2445 // not full cell along x
2446 new_mapstep_a =
2447 (int) (1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0]) + 0.5f) * mapstep_a;
2448 astart = (int) (Map_Info.xlim[0] * new_mapstep_a);
2449 aend = (int) (Map_Info.xlim[1] * new_mapstep_a);
2450 } else {
2451 new_mapstep_a = mapstep_a;
2452 astart = 0;
2453 aend = mapstep_a;
2454 }
2455 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
2456 // not full cell along y
2457 new_mapstep_b =
2458 (int) (1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0]) + 0.5f) * mapstep_b;
2459 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b);
2460 bend = (int) (Map_Info.ylim[1] * new_mapstep_b);
2461 } else {
2462 new_mapstep_b = mapstep_b;
2463 bstart = 0;
2464 bend = mapstep_b;
2465 }
2466 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
2467 // not full cell along z
2468 new_mapstep_c =
2469 (int) (1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0]) + 0.5f) * mapstep_c;
2470 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c);
2471 cend = (int) (Map_Info.zlim[1] * new_mapstep_c);
2472 } else {
2473 new_mapstep_c = mapstep_c;
2474 cstart = 0;
2475 cend = mapstep_c;
2476 }
2477 if (FourierPt)
2478 free (FourierPt);
2479 FourierPt =
2480 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c * sizeof (float));
2481 if (FourierPt == NULL) {
2482 ireslt = 0;
2483 Error_Box ("ERROR -- Unable to allocate space for map.");
2484 fclose (mapin);
2485 return;
2486 }
2487 // read the Rho values
2488
2489 for (i = astart; i < aend; i++) {
2490 for (j = bstart; j < bend; j++) {
2491 fread (&icount, sizeof (int), 1, mapin); // read header for next record
2492 if (wrong_end)
2493 icount = end_flip (icount);
2494 fread (&dens, sizeof (char), icount, mapin);
2495 fread (&icount, sizeof (int), 1, mapin); // skip trailer
2496 for (k = cstart; k < cend; k++) {
2497 ijk = i * new_mapstep_c * new_mapstep_b + j * new_mapstep_c + k;
2498 if (wrong_end)
2499 dens[k] = end_flip_real (dens[k]);
2500 FourierPt[ijk] = dens[k];
2501 if (FourierPt[ijk] < rhomin)
2502 rhomin = FourierPt[ijk];
2503 if (FourierPt[ijk] > rhomax)
2504 rhomax = FourierPt[ijk];
2505 }
2506 }
2507 }
2508 fclose (mapin);
2509 if (!Quick)
2510 fprintf (drvui->fcns, "Reading Fourier map file %s\n "
2511 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
2512 "Range of rho values: %f to %f, gridsize %d\n", infile,
2513 celda.cell[0], celda.cell[1], celda.cell[2],
2514 celda.ang[0], celda.ang[1], celda.ang[2],
2515 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax, ijk);
2516 ReadFourMap = 1;
2517 Map_Info.rhomn = rhomin;
2518 Map_Info.rhomx = rhomax;
2519 mapstep_a = new_mapstep_a;
2520 mapstep_b = new_mapstep_b;
2521 mapstep_c = new_mapstep_c;
2522 Map_Info.map_int[0] = mapstep_a;
2523 Map_Info.map_int[1] = mapstep_b;
2524 Map_Info.map_int[2] = mapstep_c;
2525 }
2526
2527 /* ************************************************************** */
2528 /* ************************************************************** */
2529
2530 void
read_grd(char * infile,int Quick)2531 read_grd (char *infile, int Quick)
2532 {
2533 FILE *mapin;
2534
2535 char line[81];
2536
2537 char *reslt;
2538
2539 float rhomin = 1.0e15f, rhomax = -1.0e15f;
2540
2541 int ireslt = 0;
2542
2543 char string[200];
2544
2545 int new_mapstep_a = 0, astart, aend;
2546
2547 int new_mapstep_b = 0, bstart, bend;
2548
2549 int new_mapstep_c = 0, cstart, cend;
2550
2551 if ((mapin = fopen (infile, "r")) == NULL) {
2552 sprintf (string, "Cannot open Fourier map (.grd) file %s", infile);
2553 Error_Box (string);
2554 return;
2555 }
2556 reslt = fgets (line, sizeof (line), mapin);
2557 trim_string (line, 81);
2558 strcpy (Map_Info.title, line);
2559 Map_Info.info_valid = 1;
2560 if (reslt != NULL) {
2561 ireslt = fscanf (mapin, "%f %f %f %f %f %f",
2562 &map_a, &map_b, &map_c, &map_alpha, &map_beta, &map_gamma);
2563 Map_Info.lat_con[0] = map_a;
2564 Map_Info.lat_con[1] = map_b;
2565 Map_Info.lat_con[2] = map_c;
2566 Map_Info.lat_con[3] = map_alpha;
2567 Map_Info.lat_con[4] = map_beta;
2568 Map_Info.lat_con[5] = map_gamma;
2569 }
2570 if (ireslt > 0)
2571 ireslt = fscanf (mapin, "%d %d %d", &mapstep_a, &mapstep_b, &mapstep_c);
2572 if (ireslt > 0) {
2573 int i, j, k, ijk;
2574
2575 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
2576 // not full cell along x
2577 new_mapstep_a =
2578 (int) (1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0]) + 0.5f) * mapstep_a;
2579 astart = (int) (Map_Info.xlim[0] * new_mapstep_a);
2580 aend = (int) (Map_Info.xlim[1] * new_mapstep_a);
2581 } else {
2582 new_mapstep_a = mapstep_a;
2583 astart = 0;
2584 aend = mapstep_a;
2585 }
2586 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
2587 // not full cell along y
2588 new_mapstep_b =
2589 (int) (1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0]) + 0.5f) * mapstep_b;
2590 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b);
2591 bend = (int) (Map_Info.ylim[1] * new_mapstep_b);
2592 } else {
2593 new_mapstep_b = mapstep_b;
2594 bstart = 0;
2595 bend = mapstep_b;
2596 }
2597 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
2598 // not full cell along z
2599 new_mapstep_c =
2600 (int) (1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0]) + 0.5f) * mapstep_c;
2601 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c);
2602 cend = (int) (Map_Info.zlim[1] * new_mapstep_c);
2603 } else {
2604 new_mapstep_c = mapstep_c;
2605 cstart = 0;
2606 cend = mapstep_c;
2607 }
2608
2609 if (FourierPt)
2610 free (FourierPt);
2611
2612 FourierPt =
2613 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c *
2614 sizeof (float));
2615 if (FourierPt == NULL) {
2616 ireslt = 0;
2617 Error_Box ("ERROR -- Unable to allocate space for map.\n");
2618 fclose (mapin);
2619 return;
2620 }
2621
2622 /* read the Rho values */
2623 for (i = astart; i < aend; i++) {
2624 for (j = bstart; j < bend; j++) {
2625 for (k = cstart; k < cend; k++) {
2626 ijk = i * new_mapstep_b * new_mapstep_c + j * new_mapstep_c + k;
2627 ireslt = fscanf (mapin, "%f", &FourierPt[ijk]);
2628 if (FourierPt[ijk] < rhomin)
2629 rhomin = FourierPt[ijk];
2630 if (FourierPt[ijk] > rhomax)
2631 rhomax = FourierPt[ijk];
2632 }
2633 }
2634 }
2635 }
2636 Map_Info.rhomn = rhomin;
2637 Map_Info.rhomx = rhomax;
2638 Map_Info.map_int[0] = mapstep_a = new_mapstep_a;
2639 Map_Info.map_int[1] = mapstep_b = new_mapstep_b;
2640 Map_Info.map_int[2] = mapstep_c = new_mapstep_c;
2641
2642 if (ireslt > 0) {
2643 if (!Quick) {
2644 fprintf (drvui->flout, "Reading Fourier map file %s\n\tTitle= %s\n "
2645 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
2646 "Range of rho values: %f to %f\n", infile, line,
2647 map_a, map_b, map_c, map_alpha, map_beta, map_gamma,
2648 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax);
2649 }
2650 ReadFourMap = 1;
2651 } else {
2652 char string[200];
2653
2654 sprintf (string, "Error reading Fourier map file %s", infile);
2655 Error_Box (string);
2656 }
2657 fclose (mapin);
2658 }
2659
2660 /* ************************************************************** */
2661 /* ************************************************************** */
2662
2663 void
read_stf(char * infile,int Quick)2664 read_stf (char *infile, int Quick)
2665 {
2666 FILE *mapin;
2667
2668 char line[81];
2669
2670 char *reslt;
2671
2672 float rhomin = 1.0e15f, rhomax = -1.0e15f;
2673
2674 int ireslt = 0;
2675
2676 int ijk = 0;
2677
2678 int i, j, k;
2679
2680 int ii, jj, kk;
2681
2682 int fullcell = 1;
2683
2684 float bx1, bx2, by1, by2, bz1, bz2;
2685
2686 char string[200];
2687
2688 int new_mapstep_a, astart, aend;
2689
2690 int new_mapstep_b, bstart, bend;
2691
2692 int new_mapstep_c, cstart, cend;
2693
2694 if ((mapin = fopen (infile, "r")) == NULL) {
2695 sprintf (string, "Cannot open Fourier map (.stf) file %s\n", infile);
2696 Error_Box (string);
2697 return;
2698 }
2699 reslt = fgets (line, sizeof (line), mapin);
2700 reslt = fgets (line, sizeof (line), mapin);
2701 reslt = fgets (line, sizeof (line), mapin);
2702 if (reslt) {
2703 ireslt = sscanf (line, "%*s %d %d %d", &mapstep_a, &mapstep_b, &mapstep_c);
2704 }
2705 if (ireslt > 0) {
2706 reslt = fgets (line, sizeof (line), mapin);
2707 sscanf (line, "%*s %f %f %f %f %f %f", &bx1, &bx2, &by1, &by2, &bz1, &bz2);
2708 Map_Info.xlim[0] = bx1;
2709 Map_Info.xlim[1] = bx2;
2710 Map_Info.ylim[0] = by1;
2711 Map_Info.ylim[1] = by2;
2712 Map_Info.zlim[0] = bz1;
2713 Map_Info.zlim[1] = bz2;
2714 Map_Info.lat_con[0] = map_a = drvui->lat_con[0];
2715 Map_Info.lat_con[1] = map_b = drvui->lat_con[1];
2716 Map_Info.lat_con[2] = map_c = drvui->lat_con[2];
2717 Map_Info.lat_con[3] = map_alpha = drvui->lat_con[3];
2718 Map_Info.lat_con[4] = map_beta = drvui->lat_con[4];
2719 Map_Info.lat_con[5] = map_gamma = drvui->lat_con[5];
2720
2721 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
2722 // not full cell along x
2723 new_mapstep_a =
2724 (int) ((1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0])) *
2725 (float) mapstep_a + 0.5f);
2726 astart = (int) (Map_Info.xlim[0] * new_mapstep_a + 0.5f);
2727 aend = (int) (Map_Info.xlim[1] * new_mapstep_a + 0.5f);
2728 fullcell = 0;
2729 if (astart < 0 && aend - astart != mapstep_a)
2730 aend++;
2731 if (astart > 0 && aend - astart != mapstep_a)
2732 aend--;
2733 if (aend - astart != mapstep_a)
2734 Error_Box ("Error recalculating map a limits for full cell");
2735 } else {
2736 new_mapstep_a = mapstep_a;
2737 astart = 0;
2738 aend = mapstep_a;
2739 }
2740 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
2741 // not full cell along y
2742 new_mapstep_b =
2743 (int) ((1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0])) *
2744 (float) mapstep_b + 0.5f);
2745 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b + 0.5f);
2746 bend = (int) (Map_Info.ylim[1] * new_mapstep_b + 0.5f);
2747 fullcell = 0;
2748 if (bstart < 0 && bend - bstart != mapstep_b)
2749 bend++;
2750 if (bstart > 0 && bend - bstart != mapstep_b)
2751 bend--;
2752 if (bend - bstart != mapstep_b)
2753 Error_Box ("Error recalculating map b limits for full cell");
2754 } else {
2755 new_mapstep_b = mapstep_b;
2756 bstart = 0;
2757 bend = mapstep_b;
2758 }
2759 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
2760 // not full cell along z
2761 new_mapstep_c =
2762 (int) ((1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0])) *
2763 (float) mapstep_c + 0.5f);
2764 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c + 0.5f);
2765 cend = (int) (Map_Info.zlim[1] * new_mapstep_c + 0.5f);
2766 fullcell = 0;
2767 if (cstart < 0 && cend - cstart != mapstep_c)
2768 cend++;
2769 if (cstart > 0 && cend - cstart != mapstep_c)
2770 cend--;
2771 if (cend - cstart != mapstep_c)
2772 Error_Box ("Error recalculating map c limits for full cell");
2773 } else {
2774 new_mapstep_c = mapstep_c;
2775 cstart = 0;
2776 cend = mapstep_c;
2777 }
2778
2779 if (FourierPt)
2780 free (FourierPt);
2781
2782 FourierPt =
2783 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c *
2784 sizeof (float));
2785 if (FourierPt == NULL) {
2786 ireslt = 0;
2787 Error_Box ("ERROR -- Unable to allocate space for map.\n");
2788 fclose (mapin);
2789 return;
2790 }
2791 reslt = fgets (line, sizeof (line), mapin);
2792 reslt = fgets (line, sizeof (line), mapin);
2793 reslt = fgets (line, sizeof (line), mapin);
2794
2795 for (i = astart; i < aend; i++) {
2796 for (j = bstart; j < bend - 1; j++) {
2797 for (k = cstart; k < cend - 1; k++) {
2798 ii = i + (i < 0 ? new_mapstep_a : 0);
2799 jj = j + (j < 0 ? new_mapstep_b : 0);
2800 kk = k + (k < 0 ? new_mapstep_c : 0);
2801 ijk =
2802 kk * (new_mapstep_a - 1) * (new_mapstep_b - 1) +
2803 jj * (new_mapstep_a - 1) + ii;
2804 ireslt = fscanf (mapin, "%E", &FourierPt[ijk]);
2805 if (FourierPt[ijk] < rhomin)
2806 rhomin = FourierPt[ijk];
2807 if (FourierPt[ijk] > rhomax)
2808 rhomax = FourierPt[ijk];
2809 }
2810 fscanf (mapin, "%*f");
2811 }
2812 for (k = 0; k < mapstep_c; k++)
2813 fscanf (mapin, "%*f");
2814 }
2815 new_mapstep_a--;
2816 new_mapstep_b--;
2817 new_mapstep_c--;
2818
2819 Map_Info.map_int[0] = mapstep_a = new_mapstep_a;
2820 Map_Info.map_int[1] = mapstep_b = new_mapstep_b;
2821 Map_Info.map_int[2] = mapstep_c = new_mapstep_c;
2822 Map_Info.rhomn = rhomin;
2823 Map_Info.rhomx = rhomax;
2824 Map_Info.info_valid = 1;
2825
2826 }
2827 /* if mapstep line present */
2828 if (ijk > 0) {
2829 if (!Quick)
2830 fprintf (drvui->flout, "Reading Fourier map file %s\n "
2831 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
2832 "Range of rho values: %f to %f, gridsize %d\n", infile,
2833 map_a, map_b, map_c, map_alpha, map_beta, map_gamma,
2834 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax, ijk);
2835 ReadFourMap = 1;
2836 } else {
2837 sprintf (string, "Error reading Fourier map file %s", infile);
2838 Error_Box (string);
2839 }
2840 fclose (mapin);
2841 }
2842
2843 /* ************************************************************** */
2844 /* ************************************************************** */
2845
2846 void
read_vasp(char * infile,int Quick)2847 read_vasp (char *infile, int Quick)
2848 {
2849 FILE *mapin;
2850
2851 char line[200];
2852
2853 char *reslt;
2854
2855 float rhomin = 1.0e15f, rhomax = -1.0e15f;
2856
2857 int ireslt = 0;
2858
2859 int ijk = 0;
2860
2861 int i, j, k, l, m;
2862
2863 float tmp_rho[10];
2864
2865 char string[200];
2866
2867 float scale;
2868
2869 float cmat[3][3];
2870
2871 int new_mapstep_a = 0, astart, aend;
2872
2873 int new_mapstep_b = 0, bstart, bend;
2874
2875 int new_mapstep_c = 0, cstart, cend;
2876
2877 map_a = drvui->lat_con[0];
2878 map_b = drvui->lat_con[1];
2879 map_c = drvui->lat_con[2];
2880 map_alpha = drvui->lat_con[3];
2881 map_beta = drvui->lat_con[4];
2882 map_gamma = drvui->lat_con[5];
2883 if ((mapin = fopen (infile, "r")) == NULL) {
2884 sprintf (string, "Cannot open VASP file %s\n", infile);
2885 Error_Box (string);
2886 return;
2887 }
2888 fgets (line, sizeof (line), mapin);
2889 trim_string (line, 80);
2890 strcpy (Map_Info.title, line);
2891 fgets (line, sizeof (line), mapin);
2892 sscanf (line, "%f", &scale);
2893 for (i = 0; i < 3; i++) { // read the cell matrix
2894 if (!(reslt = fgets (line, sizeof (line), mapin))) {
2895 sprintf (string, "Error reading VASP file.");
2896 Error_Box (string);
2897 fclose (mapin);
2898 return;
2899 }
2900 sscanf (line, "%f %f %f", &cmat[i][0], &cmat[i][1], &cmat[i][2]);
2901 }
2902 for (i = 0; i < 3; i++)
2903 for (j = 0; j < 3; j++)
2904 cmat[i][j] *= scale;
2905 fgets (line, sizeof (line), mapin);
2906 Map_Info.lat_con[0] = (float) sqrt (cmat[0][0]); // this section needs fixing
2907 Map_Info.lat_con[1] = (float) sqrt (cmat[1][1]);
2908 Map_Info.lat_con[2] = (float) sqrt (cmat[2][2]);
2909 Map_Info.lat_con[3] = 90.0f;
2910 Map_Info.lat_con[4] = 90.0f;
2911 Map_Info.lat_con[5] = 90.0f;
2912 j = 0;
2913 i = strlen (line) / 4;
2914 for (k = 0; k < i; k++) {
2915 sscanf (line + 4 * k, "%4d", &l);
2916 j += l;
2917 }
2918 for (i = 0; i < j + 2; i++) { // set to skip k-point stuff
2919 if (!(reslt = fgets (line, sizeof (line), mapin))) {
2920 sprintf (string, "Error reading VASP file.");
2921 Error_Box (string);
2922 fclose (mapin);
2923 return;
2924 }
2925 }
2926 fgets (line, sizeof (line), mapin);
2927 ireslt = sscanf (line, "%d %d %d", &mapstep_a, &mapstep_b, &mapstep_c);
2928 if (ireslt > 0) {
2929 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
2930 // not full cell along x
2931 new_mapstep_a =
2932 (int) (1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0]) + 0.5f) * mapstep_a;
2933 astart = (int) (Map_Info.xlim[0] * new_mapstep_a);
2934 aend = (int) (Map_Info.xlim[1] * new_mapstep_a);
2935 } else {
2936 new_mapstep_a = mapstep_a;
2937 astart = 0;
2938 aend = mapstep_a;
2939 }
2940 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
2941 // not full cell along y
2942 new_mapstep_b =
2943 (int) (1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0]) + 0.5f) * mapstep_b;
2944 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b);
2945 bend = (int) (Map_Info.ylim[1] * new_mapstep_b);
2946 } else {
2947 new_mapstep_b = mapstep_b;
2948 bstart = 0;
2949 bend = mapstep_b;
2950 }
2951 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
2952 // not full cell along z
2953 new_mapstep_c =
2954 (int) (1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0]) + 0.5f) * mapstep_c;
2955 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c);
2956 cend = (int) (Map_Info.zlim[1] * new_mapstep_c);
2957 } else {
2958 new_mapstep_c = mapstep_c;
2959 cstart = 0;
2960 cend = mapstep_c;
2961 }
2962
2963 if (FourierPt)
2964 free (FourierPt);
2965
2966 FourierPt =
2967 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c *
2968 sizeof (float));
2969 if (FourierPt == NULL) {
2970 ireslt = 0;
2971 Error_Box ("ERROR -- Unable to allocate space for map.\n");
2972 fclose (mapin);
2973 return;
2974 }
2975 // read the Rho values
2976
2977 m = 11;
2978 for (i = cstart; i < cend; i++) {
2979 for (j = bstart; j < bend; j++) {
2980 for (k = astart; k < aend; k++) {
2981 if (m > 9) {
2982 reslt = fgets (line, sizeof (line), mapin);
2983 sscanf (line, "%f %f %f %f %f %f %f %f %f %f", &tmp_rho[0],
2984 &tmp_rho[1], &tmp_rho[2], &tmp_rho[3], &tmp_rho[4],
2985 &tmp_rho[5], &tmp_rho[6], &tmp_rho[7], &tmp_rho[8],
2986 &tmp_rho[9]);
2987 m = 0;
2988 }
2989 ijk = k * new_mapstep_b * new_mapstep_c + j * new_mapstep_c + i;
2990 FourierPt[ijk] = tmp_rho[m++];
2991 if (FourierPt[ijk] < rhomin)
2992 rhomin = FourierPt[ijk];
2993 if (FourierPt[ijk] > rhomax)
2994 rhomax = FourierPt[ijk];
2995 }
2996 }
2997 }
2998 }
2999
3000 if (reslt) {
3001 ReadFourMap = 1;
3002 } else {
3003 char string[200];
3004
3005 sprintf (string, "Error reading Fourier map file %s", infile);
3006 Error_Box (string);
3007 }
3008 Map_Info.rhomn = rhomin;
3009 Map_Info.rhomx = rhomax;
3010 Map_Info.map_int[0] = mapstep_a = new_mapstep_a;
3011 Map_Info.map_int[1] = mapstep_b = new_mapstep_b;
3012 Map_Info.map_int[2] = mapstep_c = new_mapstep_c;
3013 if (!Quick) {
3014 fprintf (drvui->flout, "Reading Fourier map file %s\n"
3015 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
3016 "Range of rho values: %f to %f\n", infile,
3017 map_a, map_b, map_c, map_alpha, map_beta, map_gamma,
3018 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax);
3019 }
3020 Map_Info.info_valid = 1;
3021 fclose (mapin);
3022
3023 }
3024
3025 /* ************************************************************** */
3026 /* ************************************************************** */
3027
3028 void
read_w2k(char * infile,int Quick)3029 read_w2k (char *infile, int Quick)
3030 {
3031 FILE *mapin;
3032
3033 char line[81];
3034
3035 char *reslt;
3036
3037 float rhomin = 1.0e15f, rhomax = -1.0e15f;
3038
3039 int ireslt = 0;
3040
3041 int dimen = 0;
3042
3043 int ijk = 0;
3044
3045 int i, j, k;
3046
3047 int ii, jj, kk;
3048
3049 char string[200];
3050
3051 int new_mapstep_a = 0, astart, aend;
3052
3053 int new_mapstep_b = 0, bstart, bend;
3054
3055 int new_mapstep_c = 0, cstart, cend;
3056
3057 Map_Info.lat_con[0] = map_a = drvui->lat_con[0];
3058 Map_Info.lat_con[1] = map_b = drvui->lat_con[1];
3059 Map_Info.lat_con[2] = map_c = drvui->lat_con[2];
3060 Map_Info.lat_con[3] = map_alpha = drvui->lat_con[3];
3061 Map_Info.lat_con[4] = map_beta = drvui->lat_con[4];
3062 Map_Info.lat_con[5] = map_gamma = drvui->lat_con[5];
3063
3064 if ((mapin = fopen (infile, "r")) == NULL) {
3065 sprintf (string, "Cannot open Fourier map (.w2k) file %s\n", infile);
3066 Error_Box (string);
3067 return;
3068 }
3069
3070 do {
3071 reslt = fgets (line, sizeof (line), mapin);
3072 if (!strncmp (line, "(@3 [3]", 7))
3073 dimen = 1;
3074 } while (!feof (mapin) && dimen == 0);
3075
3076 if (dimen == 0) { // not in lapw5_3d format, probably from wien2venus.py
3077 rewind (mapin);
3078 reslt = fgets (line, sizeof (line), mapin);
3079 if (!strncmp (line, "cell", 4))
3080 dimen = 2;
3081 }
3082
3083 if (dimen == 0) {
3084 sprintf (string, "Cannot read Fourier map (.w2k) file %s\n", infile);
3085 Error_Box (string);
3086 fclose (mapin);
3087 return;
3088 }
3089
3090 if (dimen == 2) { // file generated by wien2venus.py
3091 reslt = fgets (line, sizeof (line), mapin);
3092 reslt = fgets (line, sizeof (line), mapin);
3093 reslt = fgets (line, sizeof (line), mapin);
3094 ireslt = sscanf (line, "%d %d %d", &mapstep_a, &mapstep_b, &mapstep_c);
3095 } else {
3096 ireslt = sscanf (line, "%*s %*s %d %d %d", &mapstep_a, &mapstep_b, &mapstep_c);
3097 }
3098
3099 if (ireslt > 0) {
3100 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
3101 // not full cell along x
3102 new_mapstep_a =
3103 (int) (1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0]) + 0.5f) * mapstep_a;
3104 astart = (int) (Map_Info.xlim[0] * new_mapstep_a);
3105 aend = (int) (Map_Info.xlim[1] * new_mapstep_a);
3106 } else {
3107 new_mapstep_a = mapstep_a;
3108 astart = 0;
3109 aend = mapstep_a;
3110 }
3111 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
3112 // not full cell along y
3113 new_mapstep_b =
3114 (int) (1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0]) + 0.5f) * mapstep_b;
3115 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b);
3116 bend = (int) (Map_Info.ylim[1] * new_mapstep_b);
3117 } else {
3118 new_mapstep_b = mapstep_b;
3119 bstart = 0;
3120 bend = mapstep_b;
3121 }
3122 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
3123 // not full cell along z
3124 new_mapstep_c =
3125 (int) (1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0]) + 0.5f) * mapstep_c;
3126 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c);
3127 cend = (int) (Map_Info.zlim[1] * new_mapstep_c);
3128 } else {
3129 new_mapstep_c = mapstep_c;
3130 cstart = 0;
3131 cend = mapstep_c;
3132 }
3133
3134 if (FourierPt)
3135 free (FourierPt);
3136
3137 FourierPt =
3138 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c *
3139 sizeof (float));
3140 if (FourierPt == NULL) {
3141 ireslt = 0;
3142 Error_Box ("ERROR -- Unable to allocate space for map.\n");
3143 fclose (mapin);
3144 return;
3145 }
3146
3147 if (dimen == 1) {
3148 reslt = fgets (line, sizeof (line), mapin);
3149
3150 for (k = cstart; k < cend; k++) {
3151 for (j = bstart; j < bend - 1; j++) {
3152 for (i = astart; i < aend - 1; i++) {
3153 ijk =
3154 j * (new_mapstep_a - 1) * (new_mapstep_b - 1) +
3155 k * (new_mapstep_a - 1) + i;
3156 ireslt = fscanf (mapin, "%f", &FourierPt[ijk]);
3157 if (FourierPt[ijk] < rhomin)
3158 rhomin = FourierPt[ijk];
3159 if (FourierPt[ijk] > rhomax)
3160 rhomax = FourierPt[ijk];
3161 }
3162 fscanf (mapin, "%*f");
3163 }
3164 for (i = 0; i < mapstep_a; i++)
3165 fscanf (mapin, "%*f");
3166 }
3167 } else {
3168 ijk = -1;
3169 for (k = astart; k < aend; k++) {
3170 for (j = bstart; j < bend - 1; j++) {
3171 for (i = cstart; i < cend - 1; i++) {
3172 ijk++;
3173
3174 if (drvui->sys == 5) { /* convert from rhombohedral setting */
3175 float xi = (float) k / (float) aend;
3176
3177 float yi = (float) j / (float) bend;
3178
3179 float zi = (float) i / (float) cend;
3180
3181 float xn = xi * 2.f / 3.f - yi * 1.f / 3.f - zi * 1.f / 3.f;
3182
3183 if (xn < 0.)
3184 xn += 1.;
3185 if (xn > 1.)
3186 xn -= 1.;
3187 float yn = xi * 1.f / 3.f + yi * 1.f / 3.f - zi * 2.f / 3.f;
3188
3189 if (yn < 0.)
3190 yn += 1.;
3191 if (yn > 1.)
3192 yn -= 1.;
3193 float zn = xi * 1.f / 3.f + yi * 1.f / 3.f + zi * 1.f / 3.f;
3194
3195 if (zn < 0.)
3196 zn += 1.;
3197 if (zn > 1.)
3198 zn -= 1.;
3199 kk = (int) (xn * aend + .5);
3200 jj = (int) (yn * bend + .5);
3201 ii = (int) (zn * cend + .5);
3202
3203 ijk =
3204 kk * (new_mapstep_c - 1) * (new_mapstep_b - 1) +
3205 jj * (new_mapstep_c - 1) + ii;
3206 }
3207
3208 ireslt = fscanf (mapin, "%f", &FourierPt[ijk]);
3209 if (FourierPt[ijk] < rhomin)
3210 rhomin = FourierPt[ijk];
3211 if (FourierPt[ijk] > rhomax)
3212 rhomax = FourierPt[ijk];
3213 }
3214 fscanf (mapin, "%*f");
3215 }
3216 for (i = 0; i < mapstep_a; i++)
3217 fscanf (mapin, "%*f");
3218 }
3219 }
3220
3221 }
3222
3223 if (ijk > 0) {
3224 if (!Quick)
3225 fprintf (drvui->flout, "Reading Fourier map file %s\n "
3226 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
3227 "Range of rho values: %f to %f, gridsize %d\n", infile,
3228 map_a, map_b, map_c, map_alpha, map_beta, map_gamma,
3229 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax, ijk);
3230 ReadFourMap = 1;
3231 } else {
3232 sprintf (string, "Error reading Fourier map file %s", infile);
3233 Error_Box (string);
3234 }
3235 fclose (mapin);
3236 Map_Info.rhomn = rhomin;
3237 Map_Info.rhomx = rhomax;
3238 Map_Info.map_int[0] = mapstep_a = --new_mapstep_a;
3239 Map_Info.map_int[1] = mapstep_b = --new_mapstep_b;
3240 Map_Info.map_int[2] = mapstep_c = --new_mapstep_c;
3241 Map_Info.info_valid = 1;
3242 }
3243
3244 /* ************************************************************** */
3245 /* ************************************************************** */
3246
3247 void
read_exc(char * infile,int Quick)3248 read_exc (char *infile, int Quick)
3249 {
3250 FILE *mapin;
3251
3252 char line[81];
3253
3254 char *reslt;
3255
3256 float rhomin, rhomax;
3257
3258 int ireslt = 0;
3259
3260 int ijk = 0;
3261
3262 int i, j, k;
3263
3264 char string[200];
3265
3266 int new_mapstep_a = 0, astart, aend;
3267
3268 int new_mapstep_b = 0, bstart, bend;
3269
3270 int new_mapstep_c = 0, cstart, cend;
3271
3272 Map_Info.lat_con[0] = map_a = drvui->lat_con[0];
3273 Map_Info.lat_con[1] = map_b = drvui->lat_con[1];
3274 Map_Info.lat_con[2] = map_c = drvui->lat_con[2];
3275 Map_Info.lat_con[3] = map_alpha = drvui->lat_con[3];
3276 Map_Info.lat_con[4] = map_beta = drvui->lat_con[4];
3277 Map_Info.lat_con[5] = map_gamma = drvui->lat_con[5];
3278
3279 rhomin = 1.0e15f;
3280 rhomax = -1.0e15f;
3281
3282 if ((mapin = fopen (infile, "r")) == NULL) {
3283 sprintf (string, "Cannot open Fourier map (.OUT) file %s\n", infile);
3284 Error_Box (string);
3285 return;
3286 }
3287
3288 ireslt = fscanf (mapin, "%d %d %d", &mapstep_a, &mapstep_b, &mapstep_c);
3289
3290 if (ireslt > 0) {
3291 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
3292 // not full cell along x
3293 new_mapstep_a =
3294 (int) (1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0]) + 0.5f) * mapstep_a;
3295 astart = (int) (Map_Info.xlim[0] * new_mapstep_a);
3296 aend = (int) (Map_Info.xlim[1] * new_mapstep_a);
3297 } else {
3298 new_mapstep_a = mapstep_a;
3299 astart = 0;
3300 aend = mapstep_a;
3301 }
3302 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
3303 // not full cell along y
3304 new_mapstep_b =
3305 (int) (1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0]) + 0.5f) * mapstep_b;
3306 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b);
3307 bend = (int) (Map_Info.ylim[1] * new_mapstep_b);
3308 } else {
3309 new_mapstep_b = mapstep_b;
3310 bstart = 0;
3311 bend = mapstep_b;
3312 }
3313 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
3314 // not full cell along z
3315 new_mapstep_c =
3316 (int) (1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0]) + 0.5f) * mapstep_c;
3317 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c);
3318 cend = (int) (Map_Info.zlim[1] * new_mapstep_c);
3319 } else {
3320 new_mapstep_c = mapstep_c;
3321 cstart = 0;
3322 cend = mapstep_c;
3323 }
3324
3325 if (FourierPt)
3326 free (FourierPt);
3327
3328 FourierPt =
3329 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c *
3330 sizeof (float));
3331 if (FourierPt == NULL) {
3332 ireslt = 0;
3333 Error_Box ("ERROR -- Unable to allocate space for map.\n");
3334 fclose (mapin);
3335 return;
3336 }
3337
3338 reslt = fgets (line, sizeof (line), mapin); // trailing comments from mapstep line
3339 fgetc (mapin);
3340
3341 for (i = cstart; i < cend; i++) {
3342 for (j = bstart; j < bend; j++) {
3343 for (k = astart; k < aend; k++) {
3344 ijk = k * new_mapstep_b * new_mapstep_c + j * new_mapstep_b + i;
3345 reslt = fgets (line, sizeof (line), mapin);
3346 ireslt = sscanf (line, "%*f %*f %*f %f", &FourierPt[ijk]);
3347 if (FourierPt[ijk] < rhomin)
3348 rhomin = FourierPt[ijk];
3349 if (FourierPt[ijk] > rhomax)
3350 rhomax = FourierPt[ijk];
3351 }
3352 }
3353 }
3354 }
3355
3356 if (ijk > 0) {
3357 if (!Quick)
3358 fprintf (drvui->flout, "Reading Fourier map file %s\n "
3359 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
3360 "Range of rho values: %f to %f, gridsize %d\n", infile,
3361 map_a, map_b, map_c, map_alpha, map_beta, map_gamma,
3362 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax, ijk);
3363 ReadFourMap = 1;
3364 } else {
3365 sprintf (string, "Error reading Fourier map file %s", infile);
3366 Error_Box (string);
3367 }
3368 fclose (mapin);
3369 Map_Info.rhomn = rhomin;
3370 Map_Info.rhomx = rhomax;
3371 Map_Info.map_int[0] = mapstep_a = new_mapstep_a;
3372 Map_Info.map_int[1] = mapstep_b = new_mapstep_b;
3373 Map_Info.map_int[2] = mapstep_c = new_mapstep_c;
3374 Map_Info.info_valid = 1;
3375 }
3376
3377 void
read_aim(char * infile,int Quick)3378 read_aim (char *infile, int Quick)
3379 {
3380 FILE *surf;
3381
3382 char line[81];
3383
3384 char *reslt;
3385
3386 int ireslt = 0;
3387
3388 int i, j, k, m;
3389
3390 int nphi, ntet;
3391
3392 int nvrt;
3393
3394 float phi, theta, r, e;
3395
3396 double sinphi, sinthe, cosphi, costhe;
3397
3398 char string[200];
3399
3400 if ((surf = fopen (infile, "r")) == NULL) {
3401 sprintf (string, "Cannot open AIM surface (.surf) file %s\n", infile);
3402 Error_Box (string);
3403 return;
3404 }
3405
3406 m = drvui->nsurf;
3407
3408 reslt = fgets (line, 80, surf);
3409 reslt = fgets (line, 80, surf);
3410 ireslt = sscanf (line, "%d", &ntet);
3411 drvui->ntet[m] = ntet;
3412 reslt = fgets (line, 80, surf);
3413 ireslt = sscanf (line, "%d", &nphi);
3414 drvui->nphi[m] = nphi;
3415 if (ireslt <= 0) {
3416 sprintf (string, "Error reading AIM surface (.surf) file %s\n", infile);
3417 Error_Box (string);
3418 fclose (surf);
3419 return;
3420 }
3421
3422 nvrt = ntet * nphi;
3423 drvui->surfx[m] = (float *) malloc (nvrt * sizeof (float));
3424 drvui->surfy[m] = (float *) malloc (nvrt * sizeof (float));
3425 drvui->surfz[m] = (float *) malloc (nvrt * sizeof (float));
3426
3427 k = 0;
3428 for (i = 0; i < ntet; ++i) {
3429 for (j = 0; j < nphi; ++j) {
3430 reslt = fgets (line, 80, surf);
3431 ireslt = sscanf (line, "%f %f %gD%f %*f", &phi, &theta, &r, &e);
3432 //if (ireslt<3) break;
3433 r *= powf (10., e);
3434 r *= 0.529f;
3435 sinphi = sin (phi);
3436 cosphi = cos (phi);
3437 sinthe = sin (theta);
3438 costhe = cos (theta);
3439
3440 drvui->surfx[m][k] = (float) (sinphi * costhe * r);
3441 drvui->surfy[m][k] = (float) (sinphi * sinthe * r);
3442 drvui->surfz[m][k] = (float) (cosphi * r);
3443 //fprintf(stderr,"%f %f %f\n",drvui->surfx[m][k],drvui->surfy[m][k],drvui->surfz[m][k]);
3444 k++;
3445 }
3446 }
3447
3448 if (!Quick)
3449 fprintf (drvui->flout, "Read %d vertices from AIM surface file %s\n",
3450 nvrt, infile);
3451
3452 fclose (surf);
3453 return;
3454
3455 }
3456
3457 /* ************************************************************** */
3458 /* ************************************************************** */
3459
3460 void
read_xsf(char * infile,int Quick)3461 read_xsf (char *infile, int Quick)
3462 {
3463 FILE *mapin;
3464
3465 char line[81];
3466
3467 char *reslt = NULL;
3468
3469 float rhomin, rhomax;
3470
3471 int ireslt = 0;
3472
3473 int ijk = 0;
3474
3475 int i, j, k;
3476
3477 char string[200];
3478
3479 int new_mapstep_a = 0, astart, aend;
3480
3481 int new_mapstep_b = 0, bstart, bend;
3482
3483 int new_mapstep_c = 0, cstart, cend;
3484
3485 Map_Info.lat_con[0] = map_a = drvui->lat_con[0];
3486 Map_Info.lat_con[1] = map_b = drvui->lat_con[1];
3487 Map_Info.lat_con[2] = map_c = drvui->lat_con[2];
3488 Map_Info.lat_con[3] = map_alpha = drvui->lat_con[3];
3489 Map_Info.lat_con[4] = map_beta = drvui->lat_con[4];
3490 Map_Info.lat_con[5] = map_gamma = drvui->lat_con[5];
3491
3492 rhomin = 1.0e15f;
3493 rhomax = -1.0e15f;
3494
3495 if ((mapin = fopen (infile, "r")) == NULL) {
3496 sprintf (string, "Cannot open Fourier map (.XSF) file %s\n", infile);
3497 Error_Box (string);
3498 return;
3499 }
3500
3501 while ( !strstr(line,"BEGIN_BLOCK_DATAGRID_3D")) {
3502 reslt = fgets (line, sizeof (line), mapin);
3503 if (feof(mapin)) {
3504 sprintf (string, "Error reading Fourier map (.XSF) file %s\n", infile);
3505 Error_Box (string);
3506 fclose (mapin);
3507 return;
3508 }
3509 }
3510 reslt = fgets (line, sizeof (line), mapin); // skip comment
3511 reslt = fgets (line, sizeof (line), mapin);
3512
3513 if ( !strstr(line,"BEGIN_DATAGRID_3D_") ) {
3514 sprintf (string, "Error reading Fourier map (.XSF) file %s\n", infile);
3515 Error_Box (string);
3516 fclose (mapin);
3517 return;
3518 }
3519
3520 reslt = fgets (line, sizeof (line), mapin); //
3521 ireslt=sscanf (line, "%d %d %d", &mapstep_a, &mapstep_b, &mapstep_c);
3522
3523 if (ireslt > 0) {
3524 if ((fabs (Map_Info.xlim[0]) > 0.001) || (fabs (Map_Info.xlim[1] - 1.0f) > 0.001)) {
3525 // not full cell along x
3526 new_mapstep_a =
3527 (int) (1.0f / (Map_Info.xlim[1] - Map_Info.xlim[0]) + 0.5f) * mapstep_a;
3528 astart = (int) (Map_Info.xlim[0] * new_mapstep_a);
3529 aend = (int) (Map_Info.xlim[1] * new_mapstep_a);
3530 } else {
3531 new_mapstep_a = mapstep_a;
3532 astart = 0;
3533 aend = mapstep_a;
3534 }
3535 if ((fabs (Map_Info.ylim[0]) > 0.001) || (fabs (Map_Info.ylim[1] - 1.0f) > 0.001)) {
3536 // not full cell along y
3537 new_mapstep_b =
3538 (int) (1.0f / (Map_Info.ylim[1] - Map_Info.ylim[0]) + 0.5f) * mapstep_b;
3539 bstart = (int) (Map_Info.ylim[0] * new_mapstep_b);
3540 bend = (int) (Map_Info.ylim[1] * new_mapstep_b);
3541 } else {
3542 new_mapstep_b = mapstep_b;
3543 bstart = 0;
3544 bend = mapstep_b;
3545 }
3546 if ((fabs (Map_Info.zlim[0]) > 0.001) || (fabs (Map_Info.zlim[1] - 1.0f) > 0.001)) {
3547 // not full cell along z
3548 new_mapstep_c =
3549 (int) (1.0f / (Map_Info.zlim[1] - Map_Info.zlim[0]) + 0.5f) * mapstep_c;
3550 cstart = (int) (Map_Info.zlim[0] * new_mapstep_c);
3551 cend = (int) (Map_Info.zlim[1] * new_mapstep_c);
3552 } else {
3553 new_mapstep_c = mapstep_c;
3554 cstart = 0;
3555 cend = mapstep_c;
3556 }
3557
3558 if (FourierPt)
3559 free (FourierPt);
3560
3561 FourierPt =
3562 (float *) zalloc (new_mapstep_a * new_mapstep_b * new_mapstep_c *
3563 sizeof (float));
3564 if (FourierPt == NULL) {
3565 ireslt = 0;
3566 Error_Box ("ERROR -- Unable to allocate space for map.\n");
3567 fclose (mapin);
3568 return;
3569 }
3570
3571 reslt = fgets (line, sizeof (line), mapin); // ignore cell vectors
3572 reslt = fgets (line, sizeof (line), mapin); // for now
3573 reslt = fgets (line, sizeof (line), mapin); //
3574 reslt = fgets (line, sizeof (line), mapin); //
3575
3576 for (i = cstart; i < cend; i++) {
3577 for (j = bstart; j < bend; j++) {
3578 for (k = astart; k < aend; k++) {
3579 ijk = k * new_mapstep_b * new_mapstep_c + j * new_mapstep_c + i;
3580 ireslt = fscanf (mapin, "%E", &FourierPt[ijk]);
3581 if (FourierPt[ijk] < rhomin)
3582 rhomin = FourierPt[ijk];
3583 if (FourierPt[ijk] > rhomax)
3584 rhomax = FourierPt[ijk];
3585 }
3586 }
3587 }
3588 }
3589
3590 if (ijk > 0) {
3591 if (!Quick)
3592 fprintf (drvui->flout, "Reading Fourier map file %s\n "
3593 "cell = %f %f %f %f %f %f\n map steps: a=%d b=%d c=%d\n "
3594 "Range of rho values: %f to %f, gridsize %d\n", infile,
3595 map_a, map_b, map_c, map_alpha, map_beta, map_gamma,
3596 mapstep_a, mapstep_b, mapstep_c, rhomin, rhomax, ijk);
3597 ReadFourMap = 1;
3598 } else {
3599 sprintf (string, "Error reading Fourier map file %s", infile);
3600 Error_Box (string);
3601 }
3602 fclose (mapin);
3603 Map_Info.rhomn = rhomin;
3604 Map_Info.rhomx = rhomax;
3605 Map_Info.map_int[0] = mapstep_a = new_mapstep_a;
3606 Map_Info.map_int[1] = mapstep_b = new_mapstep_b;
3607 Map_Info.map_int[2] = mapstep_c = new_mapstep_c;
3608 Map_Info.info_valid = 1;
3609 }
3610