#include "stdio.h"

/* ########################### wrtraw ############################ */
/*
 * Function to write raw output files.  Input array may be a two-
 * or three-dimensional array.
 * 
 * If the "field" array is two-dimensional, it should appear in the
 * function call argument as &field[0][0].  If it is three-dimensional, it
 * should appear in the argument as &field[0][0][0].
 */
void wrtraw(int nframe, int istart, int iend, int ihop,
       int jstart, int jend, int jhop, int koffst, char *base,
       float *field, int jlead, int klead, int ijob, float *xmax) {
    /* ^^^^^                                         ^^^^^
       Change these to double if appropriate. */

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

/*
 *  Local variables:
 *      i,j = dummy loop counters
 *      indx = index to a point within a 2D or 3D array
 *  Input variables:
 *      field = two- or three-dimensional array (array pointer)
 *      jlead,klead = number of columns (y) and planes (z) in array
 *      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"
 *      ijob = job number
 *        ijob = 0 -> field is 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
 */

  /* 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.  */

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