/* 
 * NOTE: Here is a more intelligent C wrtraw routine from Riaz
 * Siushansian (riaz@uda.dgcd.doc.ca).  This code, unlike the other,
 * works well when using dynamically allocated arrays.  See the bottom
 * of this file for the dynamic array-allocation functions that are
 * known to work well with this routine.  Riaz asked to make sure that
 * I remind folks that this routine comes with no guarantees (but
 * hopefully it is bug free).
 */

/*############################ wrtraw ##################################*/
/*									*/
/* Function to write raw output files.  Input array MUST be a three	*/
/* -dimensional array.							*/
/* 									*/
/* Ref. for functions used. IEEE AP. Mag. V.38, N. 6, Dec. 1996. pp. 7	*/
/* changes made by Riaz Sishansian (riaz@gauss.engga.uwo.ca), Feb. 1997.*/
/* NOTE: The code has not been tested for Full 3D arrays		*/
/* 									*/
/* In all cases "field" array is a pointer to a three-dimensional array.*/
/* Hence it must apprear as "field" in the function call argument.	*/
/*									*/
/*  Local variables:							*/
/*      i,j = dummy loop counters					*/
/*  Input variables:							*/
/*      field = pointer to a three-dimensional array (array pointer)	*/
/*      base   = base name of output files				*/
/*      nframe = frame number						*/
/*      istart,jstart = starting points in first and second "directions"*/
/*      iend,jend = ending points in first and second "directions"	*/
/*      ihop,jhop = strides in first and second "directions"		*/
/*      koffst = offset in third "direction" (i.e. constant plane)	*/
/*      ijob = job number						*/
/*        ijob = 0 -> field is a 3D array with z index of "0" 		*/
/*					  (i.e. a 2D array)		*/
/*          output = f(istart <= i <= iend, jstart <= j <= jend)	*/
/*          "x" coordinate corresponds to i				*/
/*          "y" coordinate corresponds to j				*/
/*          koffst and klead have no affect				*/
/*        ijob = 1 -> field is a 3D array, output over constant x plane	*/
/*          output = f(koffst, istart <= i <= iend, jstart <= j <= jend)*/
/*          "x" = koffst						*/
/*          "y" coordinate corresponds to i				*/
/*          "z" coordinate corresponds to j				*/
/*        ijob = 2 -> field is a 3D array, output over constant y plane	*/
/*          output = f(istart <= i <= iend, koffst, jstart <= j <= jend)*/
/*          "x" coordinate corresponds to i				*/
/*          "y" = koffst						*/
/*          "z" coordinate corresponds to j				*/
/*        ijob = 3 -> field is a 3D array, output over constant z plane	*/
/*          output = f(istart <= i <= iend, jstart <= j <= jend, koffst)*/
/*          "x" coordinate corresponds to i				*/
/*          "y" coordinate corresponds to j				*/
/*          "z" = koffst						*/
/*  Output variable:							*/
/*      xmax   = largest absolute value of output			*/
/*									*/
/*######################################################################*/

void wrtraw(int nframe, int istart, int iend, int ihop,			 
       int jstart, int jend, int jhop, int koffst, char base[127],
       float ***field, int ijob, float *xmax) {
/*     ^^^^^                     ^^^^^                                  */
/*     Change these to double if appropriate                            */

  char filnam[1024];
  float dim1, dim2, tmp;
  int i, j, k ;
  FILE *raw_out;

  /* Construct output file name. */
  sprintf(filnam,"%s.%d",base,nframe);

  /* Open the output file and write header info. */
  if((raw_out=fopen(filnam,"wb"))==NULL) {
    printf("Couldn't open %s\n", filnam);
    exit(1);
  }  

  dim1 = (float)((iend-istart)/ihop+1);
  dim2 = (float)((jend-jstart)/jhop+1);

  fwrite(&dim1,sizeof(float),1,raw_out);
  fwrite(&dim2,sizeof(float),1,raw_out);

/* Loop through writing of data.  The first entry is for the upper	*/
/* left corner, the next is for the next column over, etc.  		*/

  if (ijob == 0) {
    for (j=jend;j>=jstart;j-=jhop)
      for (i=istart;i<=iend;ihop++) {
	if(fabs(field[i][j][0]) > *xmax)
	  *xmax=fabs(field[i][j][0]);
	tmp = (float)field[i][j][0];
	fwrite(&tmp,sizeof(float),1,raw_out);
      }
  } else if (ijob == 1) {
    /* Constant x plane. */
    for (k=jend;k>=jstart;k-=jhop)
      for (j=istart;j<=iend;j+=ihop) {
	if(fabs(field[koffst][j][k]) > *xmax)
	  *xmax=fabs(field[koffst][j][k]);
	tmp = (float)field[koffst][j][k];
	fwrite(&tmp,sizeof(float),1,raw_out);
      }
  } else if (ijob == 2) {
    /* Constant y plane. */
    for (k=jend;k>=jstart;k-=jhop)
      for (i=istart;i<=iend;i+=ihop) {
	if(fabs(field[i][koffst][k]) > *xmax)
	  *xmax=fabs(field[i][koffst][k]);
	tmp = (float)field[i][koffst][k];
	fwrite(&tmp,sizeof(float),1,raw_out);
      }
  } else {
    /* Constant z plane. */
    for (j=jend;j>=jstart;j-=jhop) {
      for (i=istart;i<=iend;i+=ihop) {
	if(fabs(field[i][j][koffst]) > *xmax)
	  *xmax=fabs(field[i][j][koffst]);
	tmp = (float)field[i][j][koffst];
	fwrite(&tmp,sizeof(float),1,raw_out);
      }
    }
  }
  fclose(raw_out);
}
/* ------------------------ end of wrtraw ----------------------- */

/*=======================================================================*/
float ***matrix3d(int size_x, int size_y, int size_z) {
  int i,j;
  float ***m;
  m=(float ***) calloc((unsigned) (size_x), sizeof(float**));
  if (!m) fderror("allocation failure 1 in matrix3d()");
  for(i = 0; i < size_x; i++) {
    m[i]=(float **) calloc((unsigned) (size_y), sizeof(float*));
    if (!m[i]) fderror("allocation failure 2 in matrix3d()");
  }
  for(i = 0; i < size_x; i++) {
    for(j = 0; j < size_y; j++) {
      m[i][j]=(float *) calloc((unsigned) (size_z), sizeof(float));
      if (!m[i][j]) fderror("allocation failure 3 in matrix3d()");
    }
  }
  return m;
}
/*=======================================================================*/


/*=======================================================================*/
void free_matrix3d(float ***m, int size_x, int size_y, int size_z) {
  int i, j;
  for(i = size_x-1; i >= 0; i--) {
    for (j = size_y-1; j >= 0; j--)
      free((char*) (m[i][j]));

    free((char*) (m[i]));
  }
  free((char*) m);
}
/*=======================================================================*/

/*=======================================================================*/
void fderror(char *error_text) {
  fprintf(stderr,"\nRun-time error...\n");
  fprintf(stderr,"%s\n",error_text);
  fprintf(stderr,"...now exiting to system...\n\n");
  exit(1);
}
/*=======================================================================*/