| 2552 | } |
| 2553 | |
| 2554 | void map_data(mpb_real *d_in_re, int size_in_re, mpb_real *d_in_im, int size_in_im, int n_in[3], |
| 2555 | mpb_real *d_out_re, int size_out_re, mpb_real *d_out_im, int size_out_im, |
| 2556 | int n_out[3], matrix3x3 coord_map, mpb_real *kvector, bool pick_nearest, bool verbose, |
| 2557 | bool multiply_bloch_phase) { |
| 2558 | (void)size_in_re; |
| 2559 | (void)size_in_im; |
| 2560 | (void)size_out_re; |
| 2561 | |
| 2562 | mpb_real s[3]; /* phase difference per cell in each lattice direction */ |
| 2563 | mpb_real min_out_re = 1e20, max_out_re = -1e20, min_out_im = 1e20, max_out_im = -1e20; |
| 2564 | mpb_real shiftx, shifty, shiftz; |
| 2565 | |
| 2566 | CHECK(d_in_re && d_out_re, "invalid arguments"); |
| 2567 | CHECK((d_out_im && d_in_im) || (!d_out_im && !d_in_im), |
| 2568 | "both input and output must be real or complex"); |
| 2569 | |
| 2570 | coord_map.c0 = vector3_scale(1.0 / n_out[0], coord_map.c0); |
| 2571 | coord_map.c1 = vector3_scale(1.0 / n_out[1], coord_map.c1); |
| 2572 | coord_map.c2 = vector3_scale(1.0 / n_out[2], coord_map.c2); |
| 2573 | |
| 2574 | for (int i = 0; i < 3; ++i) { |
| 2575 | if (kvector) |
| 2576 | s[i] = kvector[i] * TWOPI; |
| 2577 | else |
| 2578 | s[i] = 0; |
| 2579 | } |
| 2580 | |
| 2581 | /* Compute shift so that the origin of the output cell |
| 2582 | is mapped to the origin of the original primitive cell: */ |
| 2583 | shiftx = 0.5 - (coord_map.c0.x * 0.5 * n_out[0] + coord_map.c1.x * 0.5 * n_out[1] + |
| 2584 | coord_map.c2.x * 0.5 * n_out[2]); |
| 2585 | shifty = 0.5 - (coord_map.c0.y * 0.5 * n_out[0] + coord_map.c1.y * 0.5 * n_out[1] + |
| 2586 | coord_map.c2.y * 0.5 * n_out[2]); |
| 2587 | shiftz = 0.5 - (coord_map.c0.z * 0.5 * n_out[0] + coord_map.c1.z * 0.5 * n_out[1] + |
| 2588 | coord_map.c2.z * 0.5 * n_out[2]); |
| 2589 | |
| 2590 | for (int i = 0; i < n_out[0]; ++i) |
| 2591 | for (int j = 0; j < n_out[1]; ++j) |
| 2592 | for (int k = 0; k < n_out[2]; ++k) { |
| 2593 | mpb_real x, y, z; |
| 2594 | double xi, yi, zi, xi2, yi2, zi2; |
| 2595 | double dx, dy, dz, mdx, mdy, mdz; |
| 2596 | int i1, j1, k1, i2, j2, k2; |
| 2597 | int ijk; |
| 2598 | |
| 2599 | ijk = (i * n_out[1] + j) * n_out[2] + k; |
| 2600 | |
| 2601 | /* find the point corresponding to d_out[i,j,k] in |
| 2602 | the input array, and also find the next-nearest |
| 2603 | points. */ |
| 2604 | x = coord_map.c0.x * i + coord_map.c1.x * j + coord_map.c2.x * k + shiftx; |
| 2605 | y = coord_map.c0.y * i + coord_map.c1.y * j + coord_map.c2.y * k + shifty; |
| 2606 | z = coord_map.c0.z * i + coord_map.c1.z * j + coord_map.c2.z * k + shiftz; |
| 2607 | MODF_POSITIVE(x, xi); |
| 2608 | MODF_POSITIVE(y, yi); |
| 2609 | MODF_POSITIVE(z, zi); |
| 2610 | |
| 2611 | if (multiply_bloch_phase) { |
no test coverage detected