/* $Id$ */

/*****************************************************************************
** 
**      Date: 28th July 2000
**
**	Author: Scott A. Belmonte
**
******************************************************************************
******************************************************************************
*     FOUE is                                                                
*           a program that extracts Fourier map information from
*           the binary Fourier map files produced by GSAS:
*           EXPNAM.DELF EXPNAM.FOBS, EXPNAM.FCLC, EXPNAM.NFDF,
*           EXPNAM.PTSN and EXPNAM.DPTS
* 
*           
*
*     Program Web-site at:                                                   
*            http://www.ccp14.ac.uk/ccp/web-mirrors/scott-belmonte-software/ 
*                                                                            
*                            Scott A. Belmonte - July 2000             
*                            s.a.belmonte@dl.ac.uk                           
*                            Department of Physics and Astronomy             
*                            University of Edinburgh                         
*                            UK.                                   
*                                                                            
******************************************************************************
******************************************************************************
*
*    Copyright (C) 2000 Scott Belmonte
*
* 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.
*
* 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.
*
*
***********************************************************************/


/*
 * $Log$
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <ctype.h>
#include <math.h>

#define PI 3.1415926536
#define MAXPHASES 10
#define MAXATOMS 1000

typedef struct fourierData {
  float *rho;        /* The actual Fourier density data for one section. */
  float romx, romn;  /* Max and min density in the section. */
} FourierData;

/* Variables contained in the Fourier map file header. */
typedef struct fourierHeader {
  char title[67], phase[67], maptyp[5];
  float cell[7];
  int nxyzi[3], nxyzo[3], nxyzt[3];
  int msect, nrho;
  float rmax, rmin, srho;
} FourierHeader;

typedef struct phaseData {
  float cell[6];
  int natoms;
  char type[MAXATOMS][3];
  float x[MAXATOMS];
  float y[MAXATOMS];
  float z[MAXATOMS];
  float occ[MAXATOMS];
  char label[MAXATOMS][12];
  float cryst2orth[12];
  float orth2cryst[12];
} PhaseData;

typedef struct wingxmapHeader {
  int nx, ny, nplanes;
  float dx, dy, dz;
  float a, b, c, alpha, beta, gamma;
  int startx, starty, startz;
  int msect, i;
  float rmin, rmax;
} WinGXMapHeader;

void output_ascii(PhaseData *phase, FourierHeader *hdr,
		  FourierData *foudata, int interactive,
		  int precision, double scale, char *expnam);
void output_wingx(PhaseData *phase, FourierHeader *hdr,
		  FourierData *foudata, int interactive,
		  char *expnam);
void output_crystals(PhaseData *phase, FourierHeader *hdr,
		     FourierData *foudata, int interactive,
		     char *expnam);
void output_grd(PhaseData *phase, FourierHeader *hdr,
		FourierData *foudata, int interactive,
		char *expnam);
void parse_expfile(char *expnam, PhaseData *data, int interactive);
void get_gsas_line(char *line, FILE *file);
float interpolate_density(PhaseData *phase, FourierHeader *hdr,
                          FourierData *data,
                          float x, float y, float z);
void get_box_limits(PhaseData *phase, FourierHeader *hdr,
		    float *xstart, float *xstep, float *xend,
		    float *ystart, float *ystep, float *yend,
		    float *zstart, float *zstep, float *zend);
void setup_trans_matrices(PhaseData *phase);
void trans_coords(float *m,
		  float *x, float *y, float *z,
		  float a, float b, float c);
void grid_to_cryst(PhaseData *phase, FourierHeader *hdr,
		   float *a, float *b, float *c,
		   int i, int j, int k);
void cryst_to_grid(PhaseData *phase, FourierHeader *hdr,
		   int *i, int *j, int *k,
 		   float a, float b, float c);

int main(int argv, char *argc[])
{

  FILE *fouin;
  char fouin_name[FILENAME_MAX];
  char expnam[80], *locate;
  char in_type[80], out_type[80];
  unsigned int ui;
  int i, type_ok, dummy_int, nRead;
  int error, verbose, interactive, precision;
  double scale;
  long int record_size, records_per_section;
  FourierHeader hdr;
  FourierData *data;
  PhaseData phase_data;

  /* Check command line to determin whether to run in
   * interactive or batch mode.
   */

  verbose = 1;

  if (argv <= 2) {
    /* Run in interactive mode. */
    interactive = 0;
    /* Check for the verbose flag. */
    if (argv == 2) { 
      if (strcmp(argc[1], "-v") == 0) {
	verbose = 0;
      }
    }

    printf("\n");
    printf("********************************************\n");
    printf("*                FOUE V1.3                 *\n");
    printf("*     GSAS Fourier Map Data Extractor.     *\n");
    printf("*                                          *\n");
    printf("*    Released under the GNU GPL Licence    *\n");
    printf("*       (c) 2000 Scott A. Belmonte         *\n");
    printf("********************************************\n");
    printf("\n");
    
    printf("Command line options: foue [options] expnam\n");
    printf("-ascii # #.# Output expnam.txt ASCII file with integer precision #\n");
    printf("             and real #.# scale factor. E.g. -ascii 2 120.56\n");
    printf("-wingx       Output expnam.map WinGX format. \n");
    printf("-march       Output expnam.fou Marching Cubes format. \n");
    printf("-xd          Output expnam.grd Project XD format.\n");
    printf("-delf        Input Fourier file format is DELF.\n");
    printf("-fobs        Input Fourier file format is FOBS.\n");
    printf("-fclc        Input Fourier file format is FCLC.\n");
    printf("-nfdf        Input Fourier file format is NFDF.\n");
    printf("-ptsn        Input Fourier file format is PTSN.\n");
    printf("-dpts        Input Fourier file format is DPTS.\n");
    printf("Defaults to FCLC input, ASCII precision 0 output if\n");
    printf("is an error while reading command line. \n\n");

    printf("Enter EXPNAM: ");
    scanf("%s", expnam);
    
    /* Get input type. */
    type_ok = 0;
    do {
      printf("Enter type of Fourier map the extract.\n");
      printf("Type (DELF, FOBS, FCLC, NFDF, PTSN or DPTS): ");
      scanf("%s", in_type);
      
      for (ui = 0; ui < strlen(in_type); ui++) 
	in_type[ui] = toupper(in_type[ui]);  
      if (strcmp(in_type, "DELF") == 0 ||
	  strcmp(in_type, "FOBS") == 0 ||
	  strcmp(in_type, "FCLC") == 0 ||
	  strcmp(in_type, "NFDF") == 0 ||
	  strcmp(in_type, "PTSN") == 0 ||
	  strcmp(in_type, "DPTS") == 0) {
	type_ok = 1;
      } else {
	printf("Type must be one of DELF, FOBS, FCLC, NFDF, PTSN or DPTS\n\n");
	type_ok = 0;
      }
    } while (type_ok == 0);
    printf("\n");
    
    /* Get output type. */
    type_ok = 0;
    do {
      printf("Enter output format type.\n");
      printf("Type     Description\n");
      printf("ASCII    Plain ASCII dump\n");
      printf("WINGX    WinGX map format (binary)\n");
      printf("MARCH    Marching Cube format (text)\n");
      printf("XD       Project XD *.grd format (text)\n");
      printf("Type: ");
      scanf("%s", out_type);
      
      for (ui = 0; ui < strlen(out_type); ui++)
	out_type[ui] = toupper(out_type[ui]);  
      if (strcmp(out_type, "ASCII") == 0 ||
	  strcmp(out_type, "WINGX") == 0 ||
	  strcmp(out_type, "MARCH") == 0 ||
	  strcmp(out_type, "XD") == 0) {
	type_ok = 1;
      } else {
	printf("Type must be one of ASCII, WINGX, MARCH or XD\n\n");
	type_ok = 0;
      }
    } while (type_ok == 0);
    printf("\n");
    
  } else {
    /* Batch mode: parse options. 
     * Default input type: FCLC
     * Default output type: ASCII, precision 0, scale 1.0
     */
    interactive = 1;
    strcpy(in_type, "FCLC");
    strcpy(out_type, "ASCII");
    precision = 0;
    scale = 1.0;

    /* The experiment name should always be the last
     * thing on the command line.
     */
    strcpy(expnam, argc[argv-1]);
    for (i = 1; i < argv - 1; i++) {
      if (strcmp(argc[i], "-a") == 0||
	  strcmp(argc[i], "-ascii") == 0) {
	if ((i + 2) > argv) {
	  precision = 0;
	  scale = 1.0;
	} else {
	  precision = atoi(argc[i+1]);
	  scale = atof(argc[i+2]);
	}
	if (precision < 0 || precision > 6) precision = 0;
	strcpy(out_type, "ASCII");
      }

      if (strcmp(argc[i], "-w") == 0 || strcmp(argc[i], "-wingx") == 0) 
	strcpy(out_type, "WINGX");
      
      if (strcmp(argc[i], "-m") == 0 || strcmp(argc[i], "-march") == 0) 
	strcpy(out_type, "MARCH");
      
      if (strcmp(argc[i], "-x") == 0 || strcmp(argc[i], "-xd") == 0) 
	strcpy(out_type, "XD");

      if (strcmp(argc[i], "-delf") == 0) 
	strcpy(in_type, "DELF");

      if (strcmp(argc[i], "-fobs") == 0) 
	strcpy(in_type, "FOBS");

      if (strcmp(argc[i], "-fclc") == 0) 
	strcpy(in_type, "FCLC");

      if (strcmp(argc[i], "-nfdf") == 0) 
	strcpy(in_type, "NFDF");

      if (strcmp(argc[i], "-ptsn") == 0) 
	strcpy(in_type, "PTSN");

      if (strcmp(argc[i], "-dpts") == 0) 
	strcpy(in_type, "DPTS");
    }
  }

  /* Check to see if the filename has an extension such as .EXP.
   * If it does, remove it.
   */
  if ((locate = strstr(expnam, ".")) != NULL) {
    *locate = '\0';
  } 
    
    
  /* First of all try opening with the 4 letter extension
   * used by the older UNIX versions of GSAS.
   */
  error = 0;
  strcpy(fouin_name, expnam);
  strcat(fouin_name, ".");
  strcat(fouin_name, in_type);
  
  /* Try to open the Fourier map file in all lower case. */
  for (ui = 0; ui < strlen(fouin_name); ui++) 
    fouin_name[ui] = tolower(fouin_name[ui]);
  if ((fouin = fopen(fouin_name, "rb")) == NULL) {
    
    /* Can't open with lower case filename so try all upper case. */ 
    for (ui = 0; ui < strlen(fouin_name); ui++)
      fouin_name[ui] = toupper(fouin_name[ui]);  
    if ((fouin = fopen(fouin_name, "rb")) == NULL) {
      error = 1;
    }
  }
  
  /* If there is no file with a four letter extension, try
   * three letter extensions.
   */
  if (error == 1) {
    in_type[3] = '\0';
    strcpy(fouin_name, expnam);
    strcat(fouin_name, ".");
    strcat(fouin_name, in_type);
    
    /* Try to open the Fourier map file in all lower case. */
    for (ui = 0; ui < strlen(fouin_name); ui++) 
      fouin_name[ui] = tolower(fouin_name[ui]);
    if ((fouin = fopen(fouin_name, "rb")) == NULL) {
      
      /* Can't open with lower case filename so try all upper case. */ 
      for (ui = 0; ui < strlen(fouin_name); ui++)
	fouin_name[ui] = toupper(fouin_name[ui]);  
      if ((fouin = fopen(fouin_name, "rb")) == NULL) {
        fprintf(stderr, "Could not find Fourier map file: %s.\n", fouin_name);
	exit(1);
      }
    }
  }
    
  /* The format of a Fourier map file is not exactly as
   * given in the GSAS manual. For one, the first header
   * actually has two 66 character strings at the start
   * (phase name and title), not one. And I'm not sure
   * whether the lattice parameter and volume info. 
   * is actually included---they are always zero in the
   * example files I have seen (I could be wrong).
   * Secondly, the file is not strictly unformatted as it
   * is split up into records. Just to make it a little
   * bit more of a challenge, the record length depends on
   * the number of entries in each section of the map, which
   * you don't know until you've read the next record, the
   * start address of which is unknown unless you know the
   * number of entries in each section! So we have to use 
   * another technique to determine the record length, which
   * relies on the x direction having greater the zero steps.
   * If is doesn't, then FOUE will have problems.
   *
   *************************************************************
   *   As of Sept 2000 there is a new version of the GSAS manual
   *   which give the correct description of the Fourier map
   *   files. There is no cell info. and records are mentioned.
   *   You still need to do some jiggery-pokery to calculate
   *   the record size.
   *
   *************************************************************
   */
    
    
  /* Read in header information. */ 
  fread(hdr.phase, sizeof(char), 66, fouin); hdr.phase[66] = '\0';
  fread(hdr.title, sizeof(char), 66, fouin); hdr.title[66] = '\0';
  fread(hdr.maptyp, sizeof(char), 4, fouin); hdr.maptyp[4] = '\0';
/*
  fread(hdr.cell, sizeof(float), 7, fouin);
*/
  if (strcmp(hdr.maptyp, "DELF") != 0 &&
      strcmp(hdr.maptyp, "FOBS") != 0 &&
      strcmp(hdr.maptyp, "FCLC") != 0 &&
      strcmp(hdr.maptyp, "NFDF") != 0 &&
      strcmp(hdr.maptyp, "PTSN") != 0 &&
      strcmp(hdr.maptyp, "DPTS") != 0) {
    fprintf(stderr, "Unknown map type: %s\n", hdr.maptyp);
    fprintf(stderr, "This may be a new type or this file is \n");
    fprintf(stderr, "is not a GSAS Fourier map file.\n\n");
    exit(1);
  }
    
  /* Determine where the start of the next record is by
   * searching for the next non-zero set of bytes in the file
   * which I assume to be nxyzi[0]. If nxyzi[0] == 0 then
   * the record length will not be calculated correctly.
   */
  do {
    nRead = fread(&dummy_int, sizeof(int), 1, fouin);
  } while (dummy_int == 0 && nRead == 1);
  
  if (nRead != 1) {
    /* We've reached the end-of-file without finding
     * another non-zero number so exit with error.
     */ 
    fprintf(stderr, "Error reading Fourier map file; no data found.\n");
    fclose(fouin);
    exit(1);
  }
  
  /* We've found a number so assume that we are at the start of 
   * the next record and read in the rest of the header information.
   */
  record_size = ftell(fouin) - sizeof(int);
  hdr.nxyzi[0] = dummy_int;
  fread(&hdr.nxyzi[1], sizeof(int), 2, fouin);
  fread(hdr.nxyzo, sizeof(int), 3, fouin);
  fread(hdr.nxyzt, sizeof(int), 3, fouin);
  fread(&hdr.msect, sizeof(int), 1, fouin);
  fread(&hdr.nrho, sizeof(int), 1, fouin);
  fread(&hdr.rmax, sizeof(float), 1, fouin);
  fread(&hdr.rmin, sizeof(float), 1, fouin);
  fread(&hdr.srho, sizeof(float), 1, fouin);
  if ((hdr.nrho + 2)*4 > record_size) {
    records_per_section = ((hdr.nrho + 2)*4)/record_size + 1;
  } else {
    records_per_section = 1;
  }
   
  /* If msect > 3 then is is most likely that the Fourier
   * map file was created on a machine with a different 
   * endianess than the machine that foue is running on.
   */
  if (hdr.msect > 3 || hdr.msect < 1) {
    fprintf(stderr, "Bad msect value: %d\n", hdr.msect);
    fprintf(stderr, "Should be 1, 2 or 3.\n");
    fprintf(stderr, "Are you sure that the map file was\n");
    fprintf(stderr, "created on this machine? See README.\n\n");
    exit(1);
  }
 
  if (verbose == 0) {
    printf("title = %s\n", hdr.title);
    printf("phase = %s\n", hdr.phase);
    printf("maptyp = %s\n", hdr.maptyp);
/*
    printf("cell: ");
    for (i = 0; i < 7; i++) printf("%f ", hdr.cell[i]);
*/
    printf("\nrecord_size = %d\n", record_size);
    printf("records_per_section = %d\n", records_per_section);
    printf("nxyzi: ");
    for (i = 0; i < 3; i++) printf("%d ", hdr.nxyzi[i]);
    printf("\nnxyzo: ");
    for (i = 0; i < 3; i++) printf("%d ", hdr.nxyzo[i]);
    printf("\nnxyzt: ");
    for (i = 0; i < 3; i++) printf("%d ", hdr.nxyzt[i]);
    printf("\nmsect = %d\n", hdr.msect);
    printf("nrho = %d\n", hdr.nrho);
    printf("rmax = %f, rmin = %f, srho = %f\n", hdr.rmax, hdr.rmin, hdr.srho);
  }
  
  /* Allocate the memory for the Fourier map data. */
  data = malloc(hdr.nxyzt[2]*sizeof(FourierData));
  if (data == NULL) {
    fprintf(stderr, "Error allocating memory for Fourier data.\n");
    fclose(fouin);
    exit(1);
  }
    
  for (i = 0; i < hdr.nxyzt[2]; i++) {
    data[i].rho = malloc(hdr.nrho*sizeof(float));
    if (data[i].rho == NULL) {
      fprintf(stderr, "Error allocating memory for Fourier sections.\n");
      free(data);
      fclose(fouin);
      exit(1);
    }
  }

  for (i = 0; i < hdr.nxyzt[2]; i++) {
    /* Jump to start of section i. */
    fseek(fouin, record_size*(2 + records_per_section*i), SEEK_SET);
    fread(data[i].rho, sizeof(float), hdr.nrho, fouin);
    fread(&data[i].romx, sizeof(float), 1, fouin);
    fread(&data[i].romn, sizeof(float), 1, fouin);    
  }
  fclose(fouin);
  
  /* Now try to open the coressponding EXPNAM.EXP file
   * to get unit cell information.
   */
  parse_expfile(expnam, &phase_data, interactive);
  
  if (interactive == 0) printf("\n");
    
    /* Now output in the chosen format. */
    if (strcmp(out_type, "ASCII") == 0) {
      output_ascii(&phase_data, &hdr, data, interactive,
		   precision, scale, expnam);
    } else if (strcmp(out_type, "WINGX") == 0) {
      output_wingx(&phase_data, &hdr, data, interactive, expnam);
    } else if (strcmp(out_type, "MARCH") == 0) {
      output_crystals(&phase_data, &hdr, data, interactive, expnam);
    } else if (strcmp(out_type, "XD") == 0) {
      if (interactive == 0)
	printf("Project XD output type not fully implemented.\n");
      output_grd(&phase_data, &hdr, data, interactive, expnam);
    }
    
    for (i = 0; i < hdr.nxyzt[2]; i++) free(data[i].rho);
    free(data);
    
    return 0;
} /* End of main(). */



void parse_expfile(char *expnam, PhaseData *data, int interactive) {

  FILE *exp_file;
  char line[81], exp_filename[FILENAME_MAX], dummy_string[80];
  unsigned int ui;
  int i, nphases, phase, natoms[MAXPHASES];
  char phase_name[MAXPHASES][80];
  float cell[MAXPHASES][6];
  float xyzo[MAXPHASES][MAXATOMS][4];
  char type[MAXPHASES][MAXATOMS][3];
  char label[MAXPHASES][MAXATOMS][12];
  strcpy(exp_filename, expnam);
  strcat(exp_filename, ".exp");

  /* Try to open the *.EXP file. */
  if ((exp_file = fopen(exp_filename, "rb")) == NULL) {
    
    /* Can't open with lower case filename so try all upper case. */ 
    for (ui = 0; ui < strlen(exp_filename); ui++) exp_filename[ui] =
      toupper(exp_filename[ui]);  
    if ((exp_file = fopen(exp_filename, "rb")) == NULL) {
      if (interactive == 0) {
	printf("Can't find %s.\n", exp_filename);
	printf("Can't get unit cell information.\n");
	printf("Setting unit cell to zero.\n");
	printf("This may make some output formats unreadable.\n");
      }
      data->cell[0] = data->cell[1] = data->cell[2] = 0.0;
      data->cell[2] = data->cell[3] = data->cell[4] = 0.0;
      data->natoms = 0;
      return;
    }
  }
  
  /* Parse the *.EXP file to get the phase information. */
  nphases = 0;
  get_gsas_line(line, exp_file);
  while (!feof(exp_file) && nphases < MAXPHASES) {
    if ((strncmp(&line[4], "    PNAM", 8)) == 0) {
      /* Read name. */
      sscanf(&line[14], "%s", phase_name[nphases]);

      /* Read number of atoms. */
      get_gsas_line(line, exp_file);
      while (strncmp(&line[6], " NATOM", 6) != 0) {
	      get_gsas_line(line, exp_file);
      }
      sscanf(&line[14], "%d", &natoms[nphases]);
      if (natoms[nphases] > MAXATOMS) {
        fprintf(stderr,
		"Error: Phase %d has %d atoms.\n", nphases, natoms[nphases]);
        fprintf(stderr, "Can't handle more than %d atoms.\n", MAXATOMS);
        fprintf(stderr, "Recompile after increasing MAXATOMS.\n");
        exit(1);
      }

      /* Read a, b, c. */
      get_gsas_line(line, exp_file);
      while (strncmp(&line[6], "ABC", 2) != 0) {
	      get_gsas_line(line, exp_file);
      }
      sscanf(&line[12], "%f%f%f%s",
	           &cell[nphases][0], &cell[nphases][1], &cell[nphases][2],
	           dummy_string);

      /* Read alpha, beta, gamma. */
      get_gsas_line(line, exp_file);
      while (strncmp(&line[6], "ANGLES", 6) != 0) {
	      get_gsas_line(line, exp_file);
      }
      sscanf(&line[12], "%f%f%f",
	           &cell[nphases][3], &cell[nphases][4], &cell[nphases][5]);

      /* Read atom info. */
      get_gsas_line(line, exp_file);
      while (strncmp(&line[6], "AT", 2) != 0) {
	      get_gsas_line(line, exp_file);
      }
      for (i = 0; i < natoms[nphases]; i++) {
	strncpy(type[nphases][i], &line[14], 2);
	type[nphases][i][2] = '\0';
        strncpy(label[nphases][i], &line[62], 11);
        label[nphases][i][11] = '\0';
        line[62] = '\0';
        sscanf(&line[22], "%f%f%f%f", 
               &xyzo[nphases][i][0], &xyzo[nphases][i][1],
               &xyzo[nphases][i][2], &xyzo[nphases][i][3]);
        get_gsas_line(line, exp_file);
	      get_gsas_line(line, exp_file);
      }
      nphases++;
    }
    get_gsas_line(line, exp_file);
  }
  
  if (interactive == 0) 
    printf("Number of phases found in EXP file: %d\n", nphases);

  if (nphases == 0) {
    if (interactive == 0) {
      printf("No phase information. Setting unit cell to zero.\n");
      printf("This may make sone output formats unreadable.\n");
    }
    data->cell[0] = data->cell[1] = data->cell[2] = 0.0;
    data->cell[2] = data->cell[3] = data->cell[4] = 0.0;
    data->natoms = 0;
    return;
  }

  phase = 0;
  if (interactive == 0) {
    for (i = 0; i < nphases; i++) {
      printf("Phase %d: Phase Name = %s\n", (i+1), phase_name[i]);
      printf("Cell: %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f\n",
	     cell[i][0], cell[i][1], cell[i][2], 
	     cell[i][3], cell[i][4], cell[i][5]);
      printf("Number of atoms in asymmetric unit: %d\n", natoms[i]);
    }
  
    if (nphases > 1) {
      do {
	printf("Which phase is the fourier map from (1 - %d): ", nphases);
	scanf("%d", &phase);
	phase--;
      } while(phase <= 0 || phase > nphases);
    } else {
      phase = 0;
    }
  }

  data->cell[0] = cell[phase][0];
  data->cell[1] = cell[phase][1];
  data->cell[2] = cell[phase][2];
  data->cell[3] = cell[phase][3];
  data->cell[4] = cell[phase][4];
  data->cell[5] = cell[phase][5];
  data->natoms = natoms[phase];
  for (i = 0; i < natoms[phase]; i++) {
    data->x[i] = xyzo[phase][i][0];
    data->y[i] = xyzo[phase][i][1];
    data->z[i] = xyzo[phase][i][2];
    data->occ[i] = xyzo[phase][i][3];
    strcpy(data->type[i], type[phase][i]);
    strcpy(data->label[i], label[phase][i]);
  }

  fclose(exp_file);
  return;
} /* End of parse_expfile(). */



void get_gsas_line(char *line, FILE *file) {

  char ch;
  int nRead;
    
  /* The line parameter must point to a character array
   * of at least 81 bytes, enough for 80 characters 
   * followed by a NULL.
   *
   *
   * Routine that reads lines from all types of GSAS 
   * ASCII files (UNIX and PC). 
   */
  nRead = fread(line, 1, 80, file);
  if (nRead == 80) {
    fread(&ch, 1, 1, file);
    switch (ch) {
    case '\n':
      /* Don't do anything. */
      break;
    case '\r':
      /* Jump forward one char because this is 
       * probably a PC format text file with \r\n
       * pairs terminating the lines. 
       */
      fseek(file, 1, SEEK_CUR);
      break;
    default:
      /* If there is no \n or \r after 80 chars then
       * this is probably a UNIX format GSAS file.
       */
      fseek(file, -1, SEEK_CUR);
      break;
    }
  }

  line[nRead] = '\0';
} /* End of get_gsas_line(). */



void output_ascii(PhaseData *phase, FourierHeader *hdr,
		  FourierData *foudata, int interactive,
		  int precision, double scale, char *expnam) {

  FILE *fouout;
  char fouout_name[FILENAME_MAX];
  int i, j, k;
  double density;

  if (interactive == 0) {
    printf("Enter output file name: ");
    scanf("%s", fouout_name);
  } else {
    strcpy(fouout_name, expnam);
    strcat(fouout_name, ".txt");
  }

  if ((fouout = fopen(fouout_name, "w")) == NULL) {
    fprintf(stderr, "Can't open %s for output.\n", fouout_name);
    exit(1);
  }

  if (interactive == 0) {
    do {
     printf("Enter precision with which to output ascii density data (0-6): ");
      scanf("%d", &precision);
    } while(precision < 0 || precision > 6);
    
    printf("Enter scale factor to multiply density data by: ");
    scanf("%lf", &scale);
  }

  /* The Fourier map consists of nxyzt[2] sections of
   * dimensions nxyzt[0] x nxyzt[1].
   */
  
  fprintf(fouout, "Unit cell: %f %f %f %f %f %f\n", 
	  phase->cell[0], phase->cell[1], phase->cell[2], 
	  phase->cell[3], phase->cell[4], phase->cell[5]);
  fprintf(fouout, "Number of atoms: %d\n", phase->natoms);
  for (i = 0; i < phase->natoms; i++) {
    fprintf(fouout, "%s %f %f %f %f\n", 
            phase->label[i],
            phase->x[i], phase->y[i], phase->z[i], phase->occ[i]);
  }
  fprintf(fouout, "Number of sections: %d\n", hdr->nxyzt[2]);
  fprintf(fouout, "Each section: %d grid points across by %d grid points down.\n",
	  hdr->nxyzt[0], hdr->nxyzt[1]);
  switch (hdr->msect) {
  case 1:
    fprintf(fouout, "Sections increase up x in steps of %f Angstroms.\n",
	    phase->cell[0]/hdr->nxyzi[2]);
    fprintf(fouout, "y increasing left to right in %f Angstrom steps.\n",
	    phase->cell[1]/hdr->nxyzi[0]);
    fprintf(fouout, "z increasing top to bottom in %f Angstrom steps.\n",
	    phase->cell[2]/hdr->nxyzi[1]);
    break;
  case 2:
    fprintf(fouout, "Sections increase up y in steps of %f Angstroms.\n",
	    phase->cell[1]/hdr->nxyzi[2]);
    fprintf(fouout, "z increasing left to right in %f Angstrom steps.\n",
	    phase->cell[2]/hdr->nxyzi[0]);
    fprintf(fouout, "x increasing top to bottom in %f Angstrom steps.\n",
	    phase->cell[0]/hdr->nxyzi[1]);
    break;
  case 3:
    fprintf(fouout, "Sections increase up z in steps of %f Angstroms.\n",
	    phase->cell[2]/hdr->nxyzi[2]);
    fprintf(fouout, "x increasing left to right in %f Angstrom steps.\n",
	    phase->cell[0]/hdr->nxyzi[0]);
    fprintf(fouout, "y increasing top to bottom in %f Angstrom steps.\n",
	    phase->cell[1]/hdr->nxyzi[1]);
    break;
  default:
    fprintf(fouout, "Unknown section direction.\n");
    break;
  }

  /* Now output the density data. */
  for (k = 0; k < hdr->nxyzt[2]; k++) {
    fprintf(fouout, "\nSection %d\n", k);
    for (j = 0; j < hdr->nxyzt[1]; j++) {
      for (i = 0; i < hdr->nxyzt[0]; i++) {
	density = scale*(foudata[k].rho[i + j*hdr->nxyzt[0]]); 
	fprintf(fouout, "%*.*f ", (precision+6), precision, density);
      }
      fprintf(fouout, "\n");
    }
  }
  fclose(fouout);
} /* End of output_ascii(). */



void output_wingx(PhaseData *phase, FourierHeader *hdr, 
		  FourierData *foudata, int interactive,
		  char *expnam) {

  FILE *fouout;
  char ch, fouout_name[FILENAME_MAX];
  WinGXMapHeader wingx_hdr;
  int i, j, k;

  if (interactive == 0) {
    printf("Enter output file name: ");
    scanf("%s", fouout_name);
  } else {
    strcpy(fouout_name, expnam);
    strcat(fouout_name, ".map");
  }

  if ((fouout = fopen(fouout_name, "wb")) == NULL) {
    fprintf(stderr, "Can't open %s for output.\n", fouout_name);
    exit(1);
  }

  wingx_hdr.nx = hdr->nxyzt[0];
  wingx_hdr.ny = hdr->nxyzt[1];
  wingx_hdr.nplanes = hdr->nxyzt[2];
  wingx_hdr.dx = (float) 1.0/(hdr->nxyzi[0]);
  wingx_hdr.dy = (float) 1.0/(hdr->nxyzi[1]);
  wingx_hdr.dz = (float) 1.0/(hdr->nxyzi[2]);
  wingx_hdr.a = phase->cell[0];
  wingx_hdr.b = phase->cell[1];
  wingx_hdr.c = phase->cell[2];
  wingx_hdr.alpha = phase->cell[3];
  wingx_hdr.beta = phase->cell[4];
  wingx_hdr.gamma = phase->cell[5];
  wingx_hdr.startx = hdr->nxyzo[0];
  wingx_hdr.starty = hdr->nxyzo[1];
  wingx_hdr.startz = hdr->nxyzo[2];
  wingx_hdr.msect = hdr->msect;
  wingx_hdr.i = 1;
  wingx_hdr.rmin = hdr->rmin;
  wingx_hdr.rmax = hdr->rmax;

  fwrite(&wingx_hdr, sizeof(WinGXMapHeader), 1, fouout);

  for (k = 0; k < hdr->nxyzt[2]; k++) {
    for (i = 0; i < hdr->nxyzt[0]; i++) {
      for (j = 0; j < hdr->nxyzt[1]; j++) {
        fwrite(&foudata[k].rho[i + j*hdr->nxyzt[0]], sizeof(float), 1, fouout);
      }
    }
  }

  ch = '\0';
  fwrite(&phase->natoms, sizeof(int), 1, fouout);
  for (i = 0; i < phase->natoms; i++) {
    /* 7 bytes of label followed by NULL char. */
    fwrite(phase->label[i], sizeof(char), 7, fouout);
    fwrite(&ch, sizeof(char), 1, fouout);
    fwrite(&phase->x[i], sizeof(float), 1, fouout);
    fwrite(&phase->y[i], sizeof(float), 1, fouout);
    fwrite(&phase->z[i], sizeof(float), 1, fouout);
  }
  fclose(fouout);
} /* End of output_wingx. */



void output_crystals(PhaseData *phase, FourierHeader *hdr,
		     FourierData *foudata, int interactive,
		     char *expnam) {

  FILE *fouout;
  char fouout_name[FILENAME_MAX];
  int i, j, k, nx, ny, nz, nxy;
  float x, y, z;
  float xstart, ystart, zstart;
  float xstep, ystep, zstep;
  float xend, yend, zend;

  if (interactive == 0) {
    printf("Enter output file name: ");
    scanf("%s", fouout_name);
  } else {
    strcpy(fouout_name, expnam);
    strcat(fouout_name, ".fou");
  }

  if ((fouout = fopen(fouout_name, "w")) == NULL) {
    fprintf(stderr, "Can't open %s for output.\n", fouout_name);
    exit(1);
  }

  setup_trans_matrices(phase);

  get_box_limits(phase, hdr,
		 &xstart, &xstep, &xend,
		 &ystart, &ystep, &yend,
		 &zstart, &zstep, &zend);

  fprintf(fouout, "INFO  DOWN, ACROSS AND SECTION \n");
  fprintf(fouout, "TRAN \n");
  for (i = 0; i < 12; i++) fprintf(fouout, "%15.8f \n", phase->cryst2orth[i]);

  fprintf(fouout, "CELL \n");
  fprintf(fouout, "%15.8f \n", phase->cell[0]);
  fprintf(fouout, "%15.8f \n", phase->cell[1]);
  fprintf(fouout, "%15.8f \n", phase->cell[2]);
  fprintf(fouout, "%15.8f \n", phase->cell[3]*PI/180.0);
  fprintf(fouout, "%15.8f \n", phase->cell[4]*PI/180.0);
  fprintf(fouout, "%15.8f \n", phase->cell[5]*PI/180.0);

  fprintf(fouout, "L14 \n");
  fprintf(fouout, "%15.8f \n", xstart);
  fprintf(fouout, "%15.8f \n", xstep);
  fprintf(fouout, "%15.8f \n", xend);
  fprintf(fouout, "%15.8f \n", 1.0);

  fprintf(fouout, "%15.8f \n", ystart);
  fprintf(fouout, "%15.8f \n", ystep);
  fprintf(fouout, "%15.8f \n", yend);
  fprintf(fouout, "%15.8f \n", 1.0);

  fprintf(fouout, "%15.8f \n", zstart);
  fprintf(fouout, "%15.8f \n", zstep);
  fprintf(fouout, "%15.8f \n", zend);
  fprintf(fouout, "%15.8f \n", 1.0);

  nx = (int) ceil((xend - xstart)/xstep);
  ny = (int) ceil((yend - ystart)/ystep);
  nz = (int) ceil((zend - zstart)/zstep);
  fprintf(fouout, "SIZE \n");
  fprintf(fouout, "%8d \n", nx);
  fprintf(fouout, "%8d \n", ny);
  fprintf(fouout, "%8d \n", nz);

  nxy = nx*ny;
  for (k = 0, z = zstart; k < nz; k++, z += zstep) {
    fprintf(fouout, "BLOCK \n");
    fprintf(fouout, "%8d \n", nxy);    
    for (i = 0, x = xstart; i < nx; i++, x += xstep) {
      for (j = 0, y = ystart; j < ny; j++, y += ystep) {
        fprintf(fouout, "%15.8f\n", 
                interpolate_density(phase, hdr, foudata, x, y, z));
      }
    }
  }

  fprintf(fouout, "LIST5 \n");
  fprintf(fouout, "%8d \n", phase->natoms);
  fprintf(fouout, "%8d \n", 14);

  for (i = 0; i < phase->natoms; i++) {
    fprintf(fouout, "%s \n", phase->type[i]);
    fprintf(fouout, "%15.8f\n", 1.0);
    fprintf(fouout, "%15.8f\n", phase->occ[i]);
    fprintf(fouout, "%15.8f\n", 0.0);
    fprintf(fouout, "%15.8f\n", phase->x[i]);
    fprintf(fouout, "%15.8f\n", phase->y[i]);
    fprintf(fouout, "%15.8f\n", phase->z[i]);
    fprintf(fouout, "%15.8f\n", 0.0);
    fprintf(fouout, "%15.8f\n", 0.0);
    fprintf(fouout, "%15.8f\n", 0.0);
    fprintf(fouout, "%15.8f\n", 0.0);
    fprintf(fouout, "%15.8f\n", 0.0);
    fprintf(fouout, "%15.8f\n", 0.0);
    fprintf(fouout, "%15.8f\n", 0.0);
  }

  fclose(fouout);
} /* End of output_crystals(). */



void output_grd(PhaseData *phase, FourierHeader *hdr, 
		FourierData *foudata, int interactive,
		char *expnam) {

  FILE *fouout;
  char fouout_name[FILENAME_MAX];
  int i, j, k, counter, nx, ny, nz;
  float x, y, z;
  float xstart, ystart, zstart;
  float xstep, ystep, zstep;
  float xend, yend, zend;

  if (interactive == 0) {
    printf("Enter output file name: ");
    scanf("%s", fouout_name);
  } else {
    strcpy(fouout_name, expnam);
    strcat(fouout_name, ".grd");
  }

  if ((fouout = fopen(fouout_name, "w")) == NULL) {
    fprintf(stderr, "Can't open %s for output.\n", fouout_name);
    exit(1);
  }

  setup_trans_matrices(phase);

  get_box_limits(phase, hdr,
		 &xstart, &xstep, &xend,
		 &ystart, &ystep, &yend,
		 &zstart, &zstep, &zend);

  nx = (int) ceil((xend - xstart)/xstep);
  ny = (int) ceil((yend - ystart)/ystep);
  nz = (int) ceil((zend - zstart)/zstep);

  fprintf(fouout, "3DGRDFIL  0\n");
  fprintf(fouout, "         ESP\n\n");
  fprintf(fouout, "! Gridpoints, Origin, Physical Dimensions\n");
  fprintf(fouout, "%15d%15d%15d\n", nx, ny, nz);
  fprintf(fouout, "%15.5E%15.5E%15.5E\n", xstart, ystart, zstart);
  fprintf(fouout, "%11.4f    %11.4f    %11.4f    \n", 
          (xend - xstart), (yend - ystart), (zend - zstart));
  fprintf(fouout, "! Objects\n");
  fprintf(fouout, "%10d\n", phase->natoms);
  for (i = 0; i < phase->natoms; i++) {
    trans_coords(phase->cryst2orth, &x, &y, &z,
		 phase->x[i], phase->y[i], phase->z[i]);		 
    fprintf(fouout, "%s%8.5f  %8.5f  %8.5f ATOM\n",
	    phase->label[i], x, y, z);
  }
  fprintf(fouout, "! Connections\n");
  fprintf(fouout, "%10d\n", 0);
  fprintf(fouout, "! Values");

  for (k = 0, z = zstart; k < nz; k++, z += zstep) {
    for (j = 0, y = ystart; j < ny; j++, y += ystep) {
      counter = 0;
      for (i = 0, x = xstart; i < nx; i++, x += xstep) {
        if (counter % 6 == 0) fprintf(fouout, "\n");
        fprintf(fouout, "%13.4E", 
                interpolate_density(phase, hdr, foudata, x, y, z));
	counter++;
      }
    }
  }

  fclose(fouout);
} /* End of output_grd(). */



float interpolate_density(PhaseData *phase, FourierHeader *hdr,
                        FourierData *data,
                        float x, float y, float z) {
  /* 3D linear interpolation of Fourier density in unit
   * cell space into cartesian space.
   */

  float a, b, c, a0, b0, c0, a1, b1, c1, da, db, dc, I[14];
  int i, j, k;

  trans_coords(phase->orth2cryst, &a, &b, &c, x, y, z);
  cryst_to_grid(phase, hdr, &i, &j, &k, a, b, c);

  /* If outside grid then return zero density
   * else interpolate from surrounding grid points.
   */
  if (i < 0 || j < 0 || k < 0) return 0.0;

  I[0] = data[k].rho[i + j*hdr->nxyzt[0]];
  I[1] = data[k+1].rho[i + j*hdr->nxyzt[0]];
  I[2] = data[k].rho[(i+1) + j*hdr->nxyzt[0]];
  I[3] = data[k+1].rho[(i+1) + j*hdr->nxyzt[0]];
  I[4] = data[k].rho[i + (j+1)*hdr->nxyzt[0]];
  I[5] = data[k+1].rho[i + (j+1)*hdr->nxyzt[0]];
  I[6] = data[k].rho[(i+1) + (j+1)*hdr->nxyzt[0]];
  I[7] = data[k+1].rho[(i+1) + (j+1)*hdr->nxyzt[0]];

  grid_to_cryst(phase, hdr, &a0, &b0, &c0, i, j, k);
  grid_to_cryst(phase, hdr, &a1, &b1, &c1, i+1, j+1, k+1);

  da = (a - a0)/(a1 - a0);
  db = (b - b0)/(b1 - b0);
  dc = (c - c0)/(c1 - c0);

  I[8] = (I[1] - I[0])*dc + I[0];
  I[9] = (I[3] - I[2])*dc + I[2];
  I[10] = (I[5] - I[4])*dc + I[4];
  I[11] = (I[7] - I[6])*dc + I[6];

  I[12] = (I[10] - I[8])*db + I[8];
  I[13] = (I[11] - I[9])*db + I[9];

  return ((I[13] - I[12])*da + I[12]);
} /* End of interpolate_density(). */



void get_box_limits(PhaseData *phase, FourierHeader *hdr,
		                float *xstart, float *xstep, float *xend,
		                float *ystart, float *ystep, float *yend,
		                float *zstart, float *zstep, float *zend) {

  /* Get the minimum and maximum extent in cartesian coordinates
   * of the GSAS grid, which is in unit cell coordinates.
   */
   
  int ci, cj, ck, i, j, k;
  float a, b, c, x, y, z;
  float xmin, ymin, zmin, xmax, ymax, zmax;

  xmin = 10000000.0;
  ymin = 10000000.0;
  zmin = 10000000.0;
  xmax = -10000000.0;
  ymax = -10000000.0;
  zmax = -10000000.0;
  
  for (ck = 0; ck <= 1; ck++) {
    for (cj = 0; cj <= 1; cj++) {
      for (ci = 0; ci <= 1; ci++) {
	      i = ci*(hdr->nxyzt[0] - 1);
 	      j = cj*(hdr->nxyzt[1] - 1);
	      k = ck*(hdr->nxyzt[2] - 1);
	      grid_to_cryst(phase, hdr, &a, &b, &c, i, j, k);
	      trans_coords(phase->cryst2orth, &x, &y, &z, a, b, c);
	      xmin = (x < xmin) ? x : xmin;
	      xmax = (x > xmax) ? x : xmax;
	      ymin = (y < ymin) ? y : ymin;
	      ymax = (y > ymax) ? y : ymax;
	      zmin = (z < zmin) ? z : zmin;
	      zmax = (z > zmax) ? z : zmax;
      }
    }
  }

  *xstart = xmin;
  *xstep = phase->cell[0]/hdr->nxyzi[0];
  *xend = xmax;
  *ystart = ymin;
  *ystep = phase->cell[1]/hdr->nxyzi[1];
  *yend = ymax;
  *zstart = zmin;
  *zstep = phase->cell[2]/hdr->nxyzi[2];
  *zend = zmax;
} /* End of get_box_limits(). */



void setup_trans_matrices(PhaseData *phase) {

  /* m[0] == m11, m[1] = m12, m[2] = m13
   * m[3] == m21, m[4] = m22, m[5] = m23
   * m[6] == m31, m[7] = m32, m[8] = m33
   * m[9] == x trans, m[10] = y trans, m[11] = z trans
   * phase->cryst2orth and phase->orth2cryst have the
   * layout above.
   * phase->cryst2orth is the matrix with which to multiply 
   * fractional unit cell coordinates to get coordinates in
   * an orthonormal coordinate system with x along a,
   * z along c* and y normal to the x-z plane.
   * phase->orth2cryst is the inverse matrix to do the
   * inverse operation.
   */

  double a, b, c, alpha, beta, gamma;
  double a_star, b_star, c_star;
  double cos_alpha, sin_alpha, cos_beta, sin_beta, cos_gamma, sin_gamma;
  double cos_alpha_star, cos_beta_star;
  double vol;

  a = phase->cell[0];
  b = phase->cell[1];
  c = phase->cell[2];
  alpha = phase->cell[3];
  beta = phase->cell[4];
  gamma = phase->cell[5];

  cos_alpha = cos(alpha*PI/180.0);
  sin_alpha = sin(alpha*PI/180.0);
  cos_beta = cos(beta*PI/180.0);
  sin_beta = sin(beta*PI/180.0);
  cos_gamma = cos(gamma*PI/180.0);
  sin_gamma = sin(gamma*PI/180.0);

  vol = a*a*b*b*c*c*
        (1 - cos_alpha*cos_alpha - cos_beta*cos_beta - cos_gamma*cos_gamma
	 + 2*cos_alpha*cos_beta*cos_gamma);
  if (vol <= 0.0) {
    printf("Error: Bad unit cell volume: vol^2 = %f\n", vol);
    printf("Unit cell params have not been read correctly.\n");
    printf("You can still output in ASCII format.\n");
    exit(1);
  }

  vol = sqrt(vol);

  a_star = b*c*sin_alpha/vol;
  b_star = a*c*sin_beta/vol;
  c_star = a*b*sin_gamma/vol;
  cos_alpha_star = (cos_beta*cos_gamma - cos_alpha)/(sin_beta*sin_gamma);
  cos_beta_star = (cos_alpha*cos_gamma - cos_beta)/(sin_alpha*sin_gamma);

  phase->cryst2orth[0] = (float) a;
  phase->cryst2orth[1] = (float) b*cos_gamma;
  phase->cryst2orth[2] = (float) c*cos_beta;
  phase->cryst2orth[3] = (float) 0.0;
  phase->cryst2orth[4] = (float) b*sin_gamma;
  phase->cryst2orth[5] = (float) -c*sin_beta*cos_alpha_star;
  phase->cryst2orth[6] = (float) 0.0;
  phase->cryst2orth[7] = (float) 0.0;
  phase->cryst2orth[8] = (float) 1.0/c_star;
  phase->cryst2orth[9] = 0.0;
  phase->cryst2orth[10] = 0.0;
  phase->cryst2orth[11] = 0.0;

  phase->orth2cryst[0] = (float) 1.0/a;
  phase->orth2cryst[1] = (float) -cos_gamma/(a*sin_gamma);
  phase->orth2cryst[2] = (float) a_star*cos_beta_star;
  phase->orth2cryst[3] = (float) 0.0;
  phase->orth2cryst[4] = (float) 1.0/(b*sin_gamma);
  phase->orth2cryst[5] = (float) b_star*cos_alpha_star;
  phase->orth2cryst[6] = (float) 0.0;
  phase->orth2cryst[7] = (float) 0.0;
  phase->orth2cryst[8] = (float) c_star;
  phase->orth2cryst[9] = (float) 0.0;
  phase->orth2cryst[10] = (float) 0.0;
  phase->orth2cryst[11] = (float) 0.0;
} /* End of setup_trans_matrices(). */



void trans_coords(float *m,
		              float *x, float *y, float *z,
		              float a, float b, float c) {
  /* Multiply a vector (a, b, c) by matrix m to
   * give a transformed vector (x, y, z).
   */

  *x = m[0]*a + m[1]*b + m[2]*c;
  *y = m[3]*a + m[4]*b + m[5]*c;
  *z = m[6]*a + m[7]*b + m[8]*c;
} /* End of trans_coords(). */



void grid_to_cryst(PhaseData *phase, FourierHeader *hdr,
		               float *a, float *b, float *c,
		               int i, int j, int k) {

  /* Convert grid point coordinates (i, j, k) into  
   * fractional coordinates (a, b, c) w.r.t the unit cell
   * axes.
   */

  switch (hdr->msect) {
  case 1:
    *a = (float) (k + hdr->nxyzo[2])/hdr->nxyzi[2];
    *b = (float) (i + hdr->nxyzo[0])/hdr->nxyzi[0];
    *c = (float) (j + hdr->nxyzo[1])/hdr->nxyzi[1];
    break;
  case 2:
    *a = (float) (j + hdr->nxyzo[1])/hdr->nxyzi[1];
    *b = (float) (k + hdr->nxyzo[2])/hdr->nxyzi[2];
    *c = (float) (i + hdr->nxyzo[0])/hdr->nxyzi[0];
    break;
  case 3:
    *a = (float) (i + hdr->nxyzo[0])/hdr->nxyzi[0];
    *b = (float) (j + hdr->nxyzo[1])/hdr->nxyzi[1];
    *c = (float) (k + hdr->nxyzo[2])/hdr->nxyzi[2];
    break;
  default:
    printf("Error: unknown section direction.\n");
    exit(1);
    break;
  }
} /* End of grid_to_cryst(). */



void cryst_to_grid(PhaseData *phase, FourierHeader *hdr,
		               int *i, int *j, int *k,
 		               float a, float b, float c) {

  /* Converts fractional unit cell coordinates (a, b, c) 
   * to the grid point nearest the origin that has
   * fractional coordinates closest to (a, b, c). 
   * If the point lies outside grid in any direction then
   * i, j, and k are all set to -1.
   */

  switch (hdr->msect) {
  case 1:
    *i = (int) floor(b*hdr->nxyzi[0] - hdr->nxyzo[0]);
    *j = (int) floor(c*hdr->nxyzi[1] - hdr->nxyzo[1]);
    *k = (int) floor(a*hdr->nxyzi[2] - hdr->nxyzo[2]);
    break;
  case 2:
    *i = (int) floor(c*hdr->nxyzi[0] - hdr->nxyzo[0]);
    *j = (int) floor(a*hdr->nxyzi[1] - hdr->nxyzo[1]);
    *k = (int) floor(b*hdr->nxyzi[2] - hdr->nxyzo[2]);
    break;
  case 3:
    *i = (int) floor(a*hdr->nxyzi[0] - hdr->nxyzo[0]);
    *j = (int) floor(b*hdr->nxyzi[1] - hdr->nxyzo[1]);
    *k = (int) floor(c*hdr->nxyzi[2] - hdr->nxyzo[2]);
    break;
  default:
    printf("Error: unknown section direction.\n");
    exit(1);
    break;
  }

  if (*i < 0 || *j < 0 || *k < 0 ||
      *i >= (hdr->nxyzt[0] - 1) ||
      *j >= (hdr->nxyzt[1] - 1) || 
      *k >= (hdr->nxyzt[2] - 1)) {
    *i = -1;
    *j = -1;
    *k = -1;
  }
} /* End of cryst_to_grid(). */