mpb-data.c

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

					    rank, out_dims2);
	  matrixio_write_real_data(data_id, out_dims2, start, 1, d_out_im);
	  matrixio_close_dataset(data_id);
     }

     if (verbose)
	  printf("Successfully wrote out data.\n");

 done:
     free(d_in_re);
     free(d_in_im);
     free(d_out_re);
     free(d_out_im);
}

void handle_file(const char *fname, const char *out_fname,
		 const char *data_name,
		 int rectify,  int have_ve, vector3 ve,
		 int resolution, real multiply_size[3],
		 int pick_nearest, int transpose)
{
     matrixio_id in_file, out_file;
     real *R, *kvector, *copies;
     int dims[2], rank;
     matrix3x3 Rin = {{1,0,0},{0,1,0},{0,0,1}}, Rout, coord_map;
#define NUM_DATANAMES 16
     char datanames[NUM_DATANAMES][30] = {
	  "data", "x", "y", "z",
	  "epsilon.xx", 
	  "epsilon.xy", 
	  "epsilon.xz", 
	  "epsilon.yy", 
	  "epsilon.yz", 
	  "epsilon.zz", 
	  "epsilon_inverse.xx", 
	  "epsilon_inverse.xy", 
	  "epsilon_inverse.xz", 
	  "epsilon_inverse.yy", 
	  "epsilon_inverse.yz", 
	  "epsilon_inverse.zz"
     };
     int i;

     if (verbose)
	  printf("Reading file %s...\n", fname);
     in_file = matrixio_open(fname, out_fname != NULL);

     if (data_name && !data_name[0])
	  data_name = NULL;

     R = matrixio_read_data_attr(in_file, "lattice vectors",
				 &rank, 2, dims);
     if (R && rank == 2 && dims[0] == 3 && dims[1] == 3) {
	  Rin.c0.x = R[0*3+0]; Rin.c0.y = R[0*3+1]; Rin.c0.z = R[0*3+2];
	  Rin.c1.x = R[1*3+0]; Rin.c1.y = R[1*3+1]; Rin.c1.z = R[1*3+2];
	  Rin.c2.x = R[2*3+0]; Rin.c2.y = R[2*3+1]; Rin.c2.z = R[2*3+2];
	  if (verbose)
	       printf("Read lattice vectors.\n");
     }
     free(R);

     kvector = matrixio_read_data_attr(in_file, "Bloch wavevector",
				       &rank, 1, dims);
     if (rank != 1 || dims[0] != 3) {
	  free(kvector);
	  kvector = NULL;
     }
     else if (verbose)
	  printf("Read Bloch wavevector (%g, %g, %g)\n",
		 kvector[0], kvector[1], kvector[2]);
     
     copies = matrixio_read_data_attr(in_file, "lattice copies",
				      &rank, 1, dims);
     if (copies && rank == 1 && dims[0] == 3) {
	  Rin.c0 = vector3_scale(copies[0], Rin.c0);
	  Rin.c1 = vector3_scale(copies[1], Rin.c1);
	  Rin.c2 = vector3_scale(copies[2], Rin.c2);
	  if (kvector) {
	       kvector[0] *= copies[0];
	       kvector[1] *= copies[1];
	       kvector[2] *= copies[2];
	  }
	  if (verbose)
	       printf("Read lattice copies (%g, %g, %g)\n",
		      copies[0], copies[1], copies[2]);
     }
     free(copies);

     if (verbose)
	  printf("Input lattice = (%g,%g,%g), (%g,%g,%g), (%g,%g,%g)\n",
		 Rin.c0.x, Rin.c0.y, Rin.c0.z,
		 Rin.c1.x, Rin.c1.y, Rin.c1.z,
		 Rin.c2.x, Rin.c2.y, Rin.c2.z);

     Rout = Rin;

     if (rectify) {
	  double V;

	  /* Orthogonalize the output lattice vectors.  If have_ve
	     is true, then the first new lattice vector should be in
	     the direction of the ve unit vector; otherwise, the first
	     new lattice vector is the first original lattice vector.
	     Note that we do this in such a way as to preserve the
	     volume of the unit cell, and so that our first vector
	     (in the direction of ve) smoothly interpolates between
	     the original lattice vectors. */

	  if (have_ve)
	       ve = unit_vector3(ve);
	  else
	       ve = unit_vector3(Rout.c0);

	  /* First, compute c0 in the direction of ve by smoothly
	     interpolating the old c0/c1/c2 (formula is slightly tricky): */
	  V = vector3_dot(vector3_cross(Rout.c0, Rout.c1), Rout.c2);
	  Rout.c1 = vector3_minus(Rout.c1,Rout.c0);
	  Rout.c2 = vector3_minus(Rout.c2,Rout.c0);
	  Rout.c0 = vector3_scale(V / vector3_dot(vector3_cross(Rout.c1,
								Rout.c2),
						  ve),
				  ve);
	  
	  /* Now, orthogonalize c1 and c2: */
	  Rout.c1 = vector3_minus(Rout.c1, 
				  vector3_scale(vector3_dot(ve, Rout.c1), ve));
	  Rout.c2 = vector3_minus(Rout.c2, 
				  vector3_scale(vector3_dot(ve, Rout.c2), ve));
	  Rout.c2 = vector3_minus(Rout.c2, 
				vector3_scale(vector3_dot(Rout.c1, Rout.c2) /
					      vector3_dot(Rout.c1, Rout.c1),
					      Rout.c1));
     }

     Rout.c0 = vector3_scale(multiply_size[0], Rout.c0);
     Rout.c1 = vector3_scale(multiply_size[1], Rout.c1);
     Rout.c2 = vector3_scale(multiply_size[2], Rout.c2);

     if (verbose)
	  printf("Output lattice = (%g,%g,%g), (%g,%g,%g), (%g,%g,%g)\n",
		 Rout.c0.x, Rout.c0.y, Rout.c0.z,
		 Rout.c1.x, Rout.c1.y, Rout.c1.z,
		 Rout.c2.x, Rout.c2.y, Rout.c2.z);

     coord_map = matrix3x3_mult(matrix3x3_inverse(Rin), Rout);

     if (out_fname) {
	  if (verbose)
	       printf("Creating output file %s...\n", out_fname);
	  out_file = matrixio_create(out_fname);
     }
     else {
	  if (verbose)
	       printf("Writing output datasets to input file %s...\n", fname);
	  out_file = in_file;
     }

     for (i = 0; i < NUM_DATANAMES; ++i) {
	  const char *dname = datanames[i];
	  char name_re[300], name_im[300];

	  if (data_name)
	       dname = data_name;

	  strcpy(name_re, dname);
	  handle_dataset(in_file, out_file, name_re, NULL,
			 Rout, coord_map, kvector, resolution,
			 multiply_size, pick_nearest, transpose);

	  sprintf(name_re, "%s.r", dname);
	  sprintf(name_im, "%s.i", dname);
	  handle_dataset(in_file, out_file, name_re, name_im,
			 Rout, coord_map, kvector, resolution,
			 multiply_size, pick_nearest, transpose);

	  if (data_name)
	       break;
     }

     free(kvector);

     matrixio_close(in_file);
     if (out_file != in_file)
	  matrixio_close(out_file);
}

void usage(FILE *f)
{
     fprintf(f, "Usage: mpb-data [options] [<filenames>]\n"
	     "Options:\n"
             "         -h : this help message\n"
             "         -V : print version number and copyright\n"
             "         -v : verbose output\n"
	     "  -o <file> : output to <file> (first input file only)\n"
	     "         -r : output rectangular cell\n"
	     " -e <x,y,z> : as -r, but first axis of cell is along <x,y,z>\n"
	     "     -n <n> : output resolution of n grid points per a\n"
	     "    -x <mx>\n"
	     "    -y <my>\n"
	     "    -z <mx> : output mx/my/mz periods in the x/y/z directions\n"
	     "     -m <s> : same as -x <s> -y <s> -z <s>\n"
	     "         -T : transpose first two dimensions (x & y) of data\n"
	     "         -p : pixellized output (no grid interpolation)\n"
	     "  -d <name> : use dataset <name> in the input files (default: all mpb datasets)\n"
	     "           -- you can also specify a dataset via <filename>:<name>\n"
	  );
}

/* given an fname of the form <filename>:<data_name>, return a pointer
   to a newly-allocated string containing <filename>, and point data_name
   to the position of <data_name> in fname.  The user must free() the
<filename> string. */
static char *split_fname(char *fname, char **data_name)
{
     int fname_len;
     char *colon, *filename;

     fname_len = strlen(fname);
     colon = strchr(fname, ':');
     if (colon) {
          int colon_len = strlen(colon);
          filename = (char*) malloc(sizeof(char) * (fname_len-colon_len+1));
          CHECK(filename, "out of memory");
          strncpy(filename, fname, fname_len-colon_len+1);
          filename[fname_len-colon_len] = 0;
          *data_name = colon + 1;
     }
else { /* treat as if ":" were at the end of fname */
          filename = (char*) malloc(sizeof(char) * (fname_len + 1));
          CHECK(filename, "out of memory");
          strcpy(filename, fname);
          *data_name = fname + fname_len;
     }
     return filename;
}

int main(int argc, char **argv)
{
     char *out_fname = NULL, *data_name = NULL;
     int rectify = 0, resolution = 0, have_ve = 0;
     vector3 ve = {1,0,0};
     real multiply_size[3] = {1,1,1};
     int pick_nearest = 0, transpose = 0;
     int ifile, c;
     extern char *optarg;
     extern int optind;

     while ((c = getopt(argc, argv, "hVvo:x:y:z:m:d:n:prTe:")) != -1)
          switch (c) {
              case 'h':
                   usage(stdout);
                   return EXIT_SUCCESS;
              case 'V':
                   printf("mpb-data " MPB_VERSION " by Steven G. Johnson.\n"
"Copyright (C) 1999, 2000, 2001, 2002, Massachusetts Institute of Technology.\n"
"This is free software, and you are welcome to redistribute it under the\n"
"terms of the GNU General Public License (GPL).  mpb-data comes with\n"
"ABSOLUTELY NO WARRANTY; see the GPL for more details.\n");
                   return EXIT_SUCCESS;
              case 'v':
                   verbose = 1;
                   break;
              case 'o':
		   free(out_fname);
                   out_fname = (char*) malloc(sizeof(char) *
                                              (strlen(optarg) + 1));
                   CHECK(out_fname, "out of memory");
                   strcpy(out_fname, optarg);
                   break;
              case 'd':
		   free(data_name);
                   data_name = (char*) malloc(sizeof(char) *
                                              (strlen(optarg) + 1));
                   CHECK(data_name, "out of memory");
                   strcpy(data_name, optarg);
                   break;
              case 'x':
                   multiply_size[0] = atof(optarg);
                   break;
              case 'y':
                   multiply_size[1] = atof(optarg);
                   break;
              case 'z':
                   multiply_size[2] = atof(optarg);
                   break;
              case 'm':
                   multiply_size[0] = atof(optarg);
                   multiply_size[1] = atof(optarg);
                   multiply_size[2] = atof(optarg);
                   break;
              case 'n':
                   resolution = atoi(optarg);
		   CHECK(resolution > 0,
			 "invalid resolution for -n (must be positive)");
                   break;
              case 'p':
                   pick_nearest = 1;
                   break;
              case 'T':
                   transpose = 1;
                   break;
	      case 'e':
		   have_ve = 1;
		   if (3 != sscanf(optarg, "%lf,%lf,%lf", 
				   &ve.x, &ve.y, &ve.z)) {
			fprintf(stderr,
				"Invalid -e argument \"%s\"\n", optarg);
			usage(stderr);
			return EXIT_FAILURE;
		   }
                   rectify = 1;
                   break;
              case 'r':
                   rectify = 1;
                   break;
              default:
                   fprintf(stderr, "Invalid argument -%c\n", c);
                   usage(stderr);
                   return EXIT_FAILURE;
          }
     if (optind == argc) {  /* no parameters left */
          usage(stderr);
          return EXIT_FAILURE;
     }
     
     for (ifile = optind; ifile < argc; ++ifile) {
	  char *dname, *h5_fname;
          h5_fname = split_fname(argv[ifile], &dname);
	  if (!dname[0])
               dname = data_name;

	  handle_file(h5_fname, out_fname, dname, 
		      rectify, have_ve, ve, resolution, 
		      multiply_size, pick_nearest, transpose);
	  
	  if (out_fname)
               free(out_fname);
          out_fname = NULL;
          free(h5_fname);
     }
     free(data_name);

     return EXIT_SUCCESS;
}