mpb-data.c

麻省理工的计算光子晶体的程序

/* Copyright (C) 1999, 2000, 2001, 2002, Massachusetts Institute of Technology.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <math.h>

#include <ctl.h>

#include "../src/config.h"

#include <check.h>
#include <matrixio.h>

int verbose = 0;

/* A macro to set x = fractional part of x input, xi = integer part,
   with 0 <= x < 1.0.   Note that we need the second test (if x >= 1.0)
   below, because x may start out as -0 or -1e-23 or something so that
   it is < 0 but x + 1.0 == 1.0, thanks to the wonders of floating point.
   (This has actually happened, on an Alpha.) */
#define MODF_POSITIVE(x, xi) { \
     x=modf(x, &xi); \
     if (x < 0) { x += 1.0; if (x >= 1.0) x = 0; else xi -= 1.0; } \
}

#define ADJ_POINT(i1, i2, nx, dx, xi, xi2) { \
     if (dx >= 0.0) { \
	  i2 = i1 + 1; \
	  if (i2 >= nx) { \
	       i2 -= nx; \
	       xi2 = xi + 1.0; \
	  } \
	  else \
	       xi2 = xi; \
     } \
     else { \
	  i2 = i1 - 1; \
	  if (i2 < 0) { \
	       i2 += nx; \
	       xi2 = xi - 1.0; \
	  } \
	  else \
	       xi2 = xi; \
          dx = -dx; \
     } \
}

void add_cmplx_times_phase(real *sum_re, real *sum_im,
			   real d_re, real d_im,
			   double ix, double iy, double iz, real *s,
			   real scale_by)
{
     static real phase = 0.0, p_re = 1.0, p_im = 0.0;
     real new_phase;

     new_phase = ix * s[0] + iy * s[1] + iz * s[2];
     if (new_phase != phase) {
	  phase = new_phase;
	  p_re = cos(phase);
	  p_im = sin(phase);
     }
     *sum_re += (d_re * p_re - d_im * p_im) * scale_by;
     *sum_im += (d_re * p_im + d_im * p_re) * scale_by;
}

#define TWOPI 6.2831853071795864769252867665590057683943388

#define MAX2(a,b) ((a) >= (b) ? (a) : (b))
#define MIN2(a,b) ((a) < (b) ? (a) : (b))

void map_data(real *d_in_re, real *d_in_im, int n_in[3], 
	      real *d_out_re, real *d_out_im, int n_out[3], 
	      matrix3x3 coord_map,
	      real *kvector,
	      short pick_nearest, short transpose)
{
     int i, j, k;
     real s[3]; /* phase difference per cell in each lattice direction */
     real min_out_re = 1e20, max_out_re = -1e20, 
	  min_out_im = 1e20, max_out_im = -1e20;
     real shiftx, shifty, shiftz;

     CHECK(d_in_re && d_out_re, "invalid arguments");
     CHECK((d_out_im && d_in_im) || (!d_out_im && !d_in_im),
	   "both input and output must be real or complex");

     coord_map.c0 = vector3_scale(1.0 / n_out[0], coord_map.c0);
     coord_map.c1 = vector3_scale(1.0 / n_out[1], coord_map.c1);
     coord_map.c2 = vector3_scale(1.0 / n_out[2], coord_map.c2);

     for (i = 0; i < 3; ++i) {
	  if (kvector)
	       s[i] = kvector[i] * TWOPI;
	  else
	       s[i] = 0;
     }

     /* Compute shift so that the origin of the output cell
	is mapped to the origin of the original primitive cell: */
     shiftx = 0.5 - (coord_map.c0.x*0.5*n_out[0] +
		     coord_map.c1.x*0.5*n_out[1] +
		     coord_map.c2.x*0.5*n_out[2]);
     shifty = 0.5 - (coord_map.c0.y*0.5*n_out[0] +
		     coord_map.c1.y*0.5*n_out[1] +
		     coord_map.c2.y*0.5*n_out[2]);
     shiftz = 0.5 - (coord_map.c0.z*0.5*n_out[0] +
		     coord_map.c1.z*0.5*n_out[1] +
		     coord_map.c2.z*0.5*n_out[2]);

     for (i = 0; i < n_out[0]; ++i)
	  for (j = 0; j < n_out[1]; ++j)
	       for (k = 0; k < n_out[2]; ++k) {
		    real x, y, z;
		    double xi, yi, zi, xi2, yi2, zi2;
		    double dx, dy, dz, mdx, mdy, mdz;
		    int i1, j1, k1, i2, j2, k2;
		    int ijk;

		    if (transpose)
			 ijk = (j * n_out[0] + i) * n_out[2] + k;
		    else
			 ijk = (i * n_out[1] + j) * n_out[2] + k;
		    
		    /* find the point corresponding to d_out[i,j,k] in
		       the input array, and also find the next-nearest
		       points. */
		    x = coord_map.c0.x*i + coord_map.c1.x*j + coord_map.c2.x*k
			 + shiftx;
		    y = coord_map.c0.y*i + coord_map.c1.y*j + coord_map.c2.y*k
			 + shifty;
		    z = coord_map.c0.z*i + coord_map.c1.z*j + coord_map.c2.z*k
			 + shiftz;
		    MODF_POSITIVE(x, xi);
		    MODF_POSITIVE(y, yi);
		    MODF_POSITIVE(z, zi);
		    i1 = x * n_in[0]; j1 = y * n_in[1]; k1 = z * n_in[2];
		    dx = x * n_in[0] - i1;
		    dy = y * n_in[1] - j1;
		    dz = z * n_in[2] - k1;
		    ADJ_POINT(i1, i2, n_in[0], dx, xi, xi2);
		    ADJ_POINT(j1, j2, n_in[1], dy, yi, yi2);
		    ADJ_POINT(k1, k2, n_in[2], dz, zi, zi2);

		    /* dx, mdx, etcetera, are the weights for the various
		       points in the input data, which we use for linearly
		       interpolating to get the output point. */
		    if (pick_nearest) {
			 /* don't interpolate */
			 dx = dx <= 0.5 ? 0.0 : 1.0;
			 dy = dy <= 0.5 ? 0.0 : 1.0;
			 dz = dz <= 0.5 ? 0.0 : 1.0;
		    }
		    mdx = 1.0 - dx;
		    mdy = 1.0 - dy;
		    mdz = 1.0 - dz;
		    
		    /* Now, linearly interpolate the input to get the
		       output.  If the input/output are complex, we
		       also need to multiply by the appropriate phase
		       factor, depending upon which unit cell we are in. */

#define IN_INDEX(i,j,k) ((i * n_in[1] + j) * n_in[2] + k)
		    if (d_out_im) {
			 d_out_re[ijk] = 0.0;
			 d_out_im[ijk] = 0.0;
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i1,j1,k1)],
					       d_in_im[IN_INDEX(i1,j1,k1)],
					       xi, yi, zi, s,
					       mdx * mdy * mdz);
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i1,j1,k2)],
					       d_in_im[IN_INDEX(i1,j1,k2)],
					       xi, yi, zi2, s,
					       mdx * mdy * dz);
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i1,j2,k1)],
					       d_in_im[IN_INDEX(i1,j2,k1)],
					       xi, yi2, zi, s,
					       mdx * dy * mdz);
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i1,j2,k2)],
					       d_in_im[IN_INDEX(i1,j2,k2)],
					       xi, yi2, zi2, s,
					       mdx * dy * dz);
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i2,j1,k1)],
					       d_in_im[IN_INDEX(i2,j1,k1)],
					       xi2, yi, zi, s,
					       dx * mdy * mdz);
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i2,j1,k2)],
					       d_in_im[IN_INDEX(i2,j1,k2)],
					       xi2, yi, zi2, s,
					       dx * mdy * dz);
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i2,j2,k1)],
					       d_in_im[IN_INDEX(i2,j2,k1)],
					       xi2, yi2, zi, s,
					       dx * dy * mdz);
			 add_cmplx_times_phase(d_out_re + ijk, d_out_im + ijk,
					       d_in_re[IN_INDEX(i2,j2,k2)],
					       d_in_im[IN_INDEX(i2,j2,k2)],
					       xi2, yi2, zi2, s,
					       dx * dy * dz);
			 min_out_im = MIN2(min_out_im, d_out_im[ijk]);
			 max_out_im = MAX2(max_out_im, d_out_im[ijk]);
		    }
		    else {
			 d_out_re[ijk] =
			      d_in_re[IN_INDEX(i1,j1,k1)] * mdx * mdy * mdz +
			      d_in_re[IN_INDEX(i1,j1,k2)] * mdx * mdy * dz +
			      d_in_re[IN_INDEX(i1,j2,k1)] * mdx * dy * mdz +
			      d_in_re[IN_INDEX(i1,j2,k2)] * mdx * dy * dz +
			      d_in_re[IN_INDEX(i2,j1,k1)] * dx * mdy * mdz +
			      d_in_re[IN_INDEX(i2,j1,k2)] * dx * mdy * dz +
			      d_in_re[IN_INDEX(i2,j2,k1)] * dx * dy * mdz +
			      d_in_re[IN_INDEX(i2,j2,k2)] * dx * dy * dz;
		    }
		    min_out_re = MIN2(min_out_re, d_out_re[ijk]);
		    max_out_re = MAX2(max_out_re, d_out_re[ijk]);
#undef IN_INDEX
	       }

     if (verbose) {
	  printf("real part range: %g .. %g\n", min_out_re, max_out_re);
	  if (d_out_im)
	       printf("imag part range: %g .. %g\n", min_out_im, max_out_im);
     }
}

void handle_dataset(matrixio_id in_file, matrixio_id out_file, 
		    const char *name_re, const char *name_im,
		    matrix3x3 Rout, matrix3x3 coord_map,
		    real *kvector, int resolution, real multiply_size[3],
		    int pick_nearest, int transpose)
{
     real *d_in_re = NULL, *d_in_im = NULL, *d_out_re = NULL, *d_out_im = NULL;
     int in_dims[3] = {1,1,1}, out_dims[3] = {1,1,1}, out_dims2[3], rank = 3;
     int i, N;
     int start[3] = {0,0,0};
     matrixio_id data_id;
     char out_name[1000];

     d_in_re = matrixio_read_real_data(in_file, name_re, &rank, in_dims,
				       0, 0, 0, NULL);
     if (!d_in_re)
	  goto done;

     if (verbose)
	  printf("Found dataset %s...\n", name_re);

     if (name_im) {
	  d_in_im = matrixio_read_real_data(in_file, name_im, &rank, out_dims,
					    0, 0, 0, NULL);
	  if (!d_in_im) {
	       fprintf(stderr, "mpb-data: found %s dataset but not %s\n",
		       name_re, name_im);
	       goto done;
	  }
	  
	  for (i = 0; i < 3; ++i) {
	       CHECK(out_dims[i] == in_dims[i],
		     "re/im datasets must have same size!");
	  }

	  if (verbose)
	       printf("   and imaginary part dataset %s...\n", name_im);

     }

     if (verbose)
	  printf("Input data is rank %d, size %dx%dx%d.\n",
		 rank, in_dims[0], in_dims[1], in_dims[2]);

     if (resolution > 0) {
	  out_dims[0] = vector3_norm(Rout.c0) * resolution + 0.5;
	  out_dims[1] = vector3_norm(Rout.c1) * resolution + 0.5;
	  out_dims[2] = vector3_norm(Rout.c2) * resolution + 0.5;
     }
     else {
	  for (i = 0; i < 3; ++i)
	       out_dims[i] = in_dims[i] * multiply_size[i];
     }
     for (i = rank; i < 3; ++i)
	  out_dims[i] = 1;
     for (N = 1, i = 0; i < 3; ++i)
	  N *= (out_dims[i] = MAX2(out_dims[i], 1));

     if (transpose) {
	  out_dims2[0] = out_dims[1];
	  out_dims2[1] = out_dims[0];
	  out_dims2[2] = out_dims[2];
     }
     else {
	  out_dims2[0] = out_dims[0];
	  out_dims2[1] = out_dims[1];
	  out_dims2[2] = out_dims[2];
     }

     if (verbose)
	  printf("Output data %dx%dx%d.\n",
		 out_dims2[0], out_dims2[1], out_dims2[2]);

     CHK_MALLOC(d_out_re, real, N);
     if (d_in_im) {
	  CHK_MALLOC(d_out_im, real, N);
     }

     map_data(d_in_re, d_in_im, in_dims, d_out_re, d_out_im, out_dims,
	      coord_map, kvector, pick_nearest, transpose);

     strcpy(out_name, name_re);
     if (out_file == in_file)
	  strcat(out_name, "-new");
     if (verbose)
	  printf("Writing dataset to %s...\n", out_name);
     data_id = matrixio_create_dataset(out_file, out_name,"", rank, out_dims2);
     matrixio_write_real_data(data_id, out_dims2, start, 1, d_out_re);
     matrixio_close_dataset(data_id);

     if (d_out_im) {
	  strcpy(out_name, name_im);
	  if (out_file == in_file)
	       strcat(out_name, "-new");
	  if (verbose)
	       printf("Writing dataset to %s...\n", out_name);
	  data_id = matrixio_create_dataset(out_file, out_name, "",