/* mesh.c -- mesh handling routines
//
// Written and Copyright (C) 1994, 1996 by Michael J. Gourlay
//
// NO WARRANTEES, EXPRESS OR IMPLIED.
*/

#include <stdio.h>
#include <stdlib.h>

#include "my_malloc.h"
#include "file.h"
#include "mesh.h"




/* ---------------------------------------------------------------------- */

/* meshFree: free memory of internal arrays of a MeshT
//
// Note that the memory of the MeshT instance is not freed here.
//
// proc (in): name of the calling procedure (for debugging.
*/
void
meshFree(MeshT *meshP, char *proc)
{
#ifdef DEBUG_ALLOC
  printf("freeing mesh %p %p from %s\n", meshP->x, meshP->y, proc);
#endif

  if(meshP->x != NULL) {
    my_free(meshP->x, "meshFree");
    meshP->x = NULL;
  }

  if(meshP->y != NULL) {
    my_free(meshP->y, "meshFree");
    meshP->y = NULL;
  }
}




/* ---------------------------------------------------------------------- */

/* meshAlloc: allocate memory for internal arrays of a MeshT
//
// meshP (in/out): mesh.  meshP->nx and meshP->ny are used on input to
//   determine mesh size.  Memory for meshP->x and meshP->y are allocated.
//
// proc (in): name of the calling procedure (for debugging)
*/
int
meshAlloc(MeshT *meshP, const char *proc)
{
  if(meshP->nx <= 0 || meshP->ny <= 0) {
    fprintf(stderr, "meshAlloc: non-positive size\n");
    return 1;
  }

  if(meshP->x != NULL || meshP->y != NULL) {
    fprintf(stderr,
            "meshAlloc: warning: allocating over un-freed mesh from %s\n",
            proc);
  }

  /* Allocate mesh arrays */
  if((meshP->x=MY_CALLOC(meshP->nx * meshP->ny, double, "meshAlloc"))==NULL) {
    fprintf(stderr, "%s: Bad Alloc\n", proc);
    return 1;
  }
  if((meshP->y=MY_CALLOC(meshP->nx * meshP->ny, double, "meshAlloc"))==NULL) {
    my_free(meshP->x, "meshAlloc");
    fprintf(stderr, "%s: Bad Alloc\n", proc);
    return 1;
  }

#ifdef DEBUG_ALLOC
  printf("allocated mesh %p %p from %s\n", meshP->x, meshP->y, proc);
#endif

  return 0;
}




/* ---------------------------------------------------------------------- */

/* meshInitialize: Initialize a mesh
//
// As each MeshT instance is created, it should be initialized
// using this routine, even if the size is zero.
//
// If the size is non-zero then memory for the mesh arrays is allocated here.
*/
void
meshInitialize(MeshT *meshP, int nx, int ny)
{
  /* Set the mesh size */
  meshP->nx = nx;
  meshP->ny = ny;

  if(nx*ny > 0) {
    meshP->x = meshP->y = NULL;
    meshAlloc(meshP, "meshInitialize");
  } else {
    /* Set arrays to NULL to indicate that they are not allocated */
    meshP->x = meshP->y = NULL;
  }
}




/* ---------------------------------------------------------------------- */
/* meshPrint: print some info about a mesh, for debugging
*/
void
meshPrint(const MeshT *meshP)
{
  printf("size=%i,%i    max corner=%lg,%lg\n", meshP->nx, meshP->ny,
         meshP->x[meshP->nx * meshP->ny - 1],
         meshP->y[meshP->nx * meshP->ny - 1]);
}




/* ---------------------------------------------------------------------- */

/* meshCheckCompatibility: make sure two meshes are compatible
//
// Two meshes must be compatible in order to create a "tween" mesh.
//
// In order for meshes to be compatible, they must have the same
// number of points in each direction.
//
// If both meshes are the "input" meshes for an interpolation, they
// should also have the same maximum values, i.e., their corners should
// be at the same places.  However, this routine will be used to check
// for compatibility between input and output meshes, in which case the
// output mesh will initially have no mesh point values set.
//
//
// Return nonzero if meshes don't match.
// Return zero if they match.
*/
int
meshCheckCompatibility(const MeshT *mesh1, const MeshT *mesh2)
{
  if(mesh1->nx != mesh2->nx) {
    return 1;
  } else if(mesh1->ny != mesh2->ny) {
    return 2;
  }

#ifdef MESH_CHECK_CORNERS
  if(mesh1->x[0] != mesh2->x[0]) {
    return 3;
  } else if(mesh1->x[mesh1->nx-1] != mesh2->x[mesh1->nx-1]) {
    printf("mCC: %lg %lg\n", mesh1->x[mesh1->nx-1], mesh2->x[mesh1->nx-1]);
    return 4;
  } else if(mesh1->y[0] != mesh2->y[0]) {
    return 5;
  } else if(mesh1->y[mesh1->ny-1] != mesh2->y[mesh1->ny-1]) {
    return 6;
  }
#endif /* MESH_CHECK_CORNERS */
  return 0;
}




/* ---------------------------------------------------------------------- */

/* meshChannelLinInterp: linear interpolation between two meshes, channel
//
// Note that this routine only operates on a single channel of the
// meshes.  (A mesh has two channels: the list of x-values and the list
// of y-values.)
//
// mi1 (in): a channel of input mesh 1
// mi2 (in): a channel of input mesh 2
// nx (in): number of mesh points in x-direction
// my (in): number of mesh points in y-direction
// mo (out): the respective channel of the output mesh
// t (in): tween parameter:
//   when 0<t<1, mo = (1-t)*mi1 + t*mi2.
//   e.g. when t==0, mo = mi1.  when t==1, mo = mi2.
//
// MJG 18jul94:
// The roundoff error here, although small, sometimes triggers the
// bounds check in the spline evaluator.  The effect should be
// harmless, though, if the spline evaluates slightly out of range.
*/
static void
meshChannelLinInterp(const double *mi1, const double *mi2, int nx, int ny, double t, double *mo)
{
  int xi, yi;

  for(yi=0; yi<ny; yi++) {
    for(xi=0; xi<nx; xi++) {
      mo[yi*nx + xi] = (1.0-t)*mi1[yi*nx + xi] + t*mi2[yi*nx + xi];
    }
  }
}




/* ---------------------------------------------------------------------- */

/* meshInterpolate: interpolate meshes
// See meshChannelLinInterp for semantics of tween_param.
//
// m1p (in): "source" mesh pointer
//
// m2p (in): "destination" mesh pointer
//
// tween_param: (in) parameter indicating output mesh configuration.
//   0 < tween_param < 1
//   When tween_param == 0, moP is m1p.
//   When tween_param == 1, moP is m2p.
//
// moP (out): "tween" mesh, somewhere between m1P and m2P.
*/
void
meshInterpolate(MeshT *moP, const MeshT *m1P, const MeshT *m2P, double tween_param)
{
  int m_c_c;

  if(m_c_c = meshCheckCompatibility(m1P, m2P)) {
    fprintf(stderr, "meshInterpolate: input mesh sizes mismatch %i\n", m_c_c);
    return;
  }

  if(m_c_c = meshCheckCompatibility(m1P, moP)) {
    fprintf(stderr, "meshInterpolate: input mesh size mismatches output mesh %i\n", m_c_c);
    return;
  }

  meshChannelLinInterp(m1P->x, m2P->x, m1P->nx, m1P->ny, tween_param, moP->x);
  meshChannelLinInterp(m1P->y, m2P->y, m1P->nx, m1P->ny, tween_param, moP->y);
}




/* ---------------------------------------------------------------------- */

/* meshReset : set image warp mesh to be a regularly spaced mesh
//
// mP: (in/out) mesh pointer
//
// img_width, img_height: width and height of the image that goes with
//   this mesh
//
//
// There ought to be a meshRescale, which simply rescales the mesh to
// fit the image geometry.  In many cases where meshReset is currently called, the
// meshRescale should probably be called instead, because a rescale is
// less destructive than a meshReset. -- 29nov96
*/
void
meshReset(MeshT *mP, const int img_width, const int img_height)
{
  int xi, yi;
  const double mp_dx = (double)(img_width - 1)  / (double)(mP->nx - 1) ;
  const double mp_dy = (double)(img_height - 1) / (double)(mP->ny - 1) ;

  for(yi=0; yi < mP->ny; yi++) {
    for(xi=0; xi < mP->nx; xi++) {
      mP->x[yi * mP->nx + xi] = (double)((int)(mp_dx*(double)xi + 0.5));
      mP->y[yi * mP->nx + xi] = (double)((int)(mp_dy*(double)yi + 0.5));
    }
  }
}




/* ---------------------------------------------------------------------- */

/* meshFunctionalize : set image warp mesh to be functional and bounded
//
// mP (in/out): mesh pointer
// mw, mh: width and height of the mesh's image
//
// Return number of changes.
//
// This routine only enforces vertical and horizontal functional lines.
// (I.e. lines can cross diagonally.)
//
// The problem with this routine is that if a point is out of its box,
// it's mathematically ambiguous whether that point should be moved,
// or whether the adjacent point should be moved.  Moving either fixes
// the "functionality", but usually there is an intuitive choice which
// this algorithm does not see.  This algorithm moves both points.
// To fix this problem, a heuristic could be employed to place a point
// within some weighted average of its neighbors.  Another posibility
// would be to generate the spline and use the values the spline is
// forced to use.  The problem with this is that the spline already
// expects functional data, so forcing a spline might break the
// spline.  Yet another possibility is to make a first-pass through
// the data to see which points need fixing, by the criteria used in
// this routine, along with an additional backwards criterion to ensure
// symmetry..  Then, the second pass would weight the changes according to
// which points require changes.  This would, at least, keep major
// crossovers effects localized (such as when a single point crosses
// over many points).
//
// This could be looked at as a bare-bones functionalizer-- one which
// simply guarentees that meshes are functional.  It is probably the
// job of another algorithm to make the mesh look more like what the
// user intended, but the user should have done what it intended in the
// first place...
*/
int
meshFunctionalize(MeshT *mP, int mw, int mh)
{
  register int xi, yi;
  double mxv, myv;
  int loop_change;
  int mesh_change=0;

  /* Repeat the mesh changes until the mesh stops changing */
  /* (but stop trying after a while to avoid long or infinite loops) */
  do {
    loop_change = 0;

    /* Force top and bottom edges to be at borders */
    for(xi=0; xi < mP->nx; xi++) {
      if(mP->y[xi]!=0) {
        mP->y[xi] = 0;
        loop_change++;
      }
      if(mP->y[(mP->ny-1) * mP->nx + xi] != (mh-1)) {
        mP->y[(mP->ny-1) * mP->nx + xi] = mh-1;
        loop_change++;
      }
    }

    mP->y[0] = 0;
    for(yi=1; yi < mP->ny; yi++) {
      /* Force left and right edges to be at borders */
      if(mP->x[yi * mP->nx + 0] != 0) {
        mP->x[yi* mP->nx + 0] = 0;
        loop_change++;
      }
      if(mP->x[yi* mP->nx + (mP->nx-1)] != (mw-1)) {
        mP->x[yi* mP->nx + (mP->nx-1)] = mw-1;
        loop_change++;
      }

      /* Enforce functionality */
      for(xi=1; xi < mP->nx; xi++) {
        /* make current point right of previous point */
        if(mP->x[yi* mP->nx + xi] <= mP->x[yi* mP->nx + (xi-1)]) {
          mxv = (mP->x[yi* mP->nx + xi] + mP->x[yi* mP->nx + (xi-1)]) / 2;
          mP->x[yi* mP->nx + xi]     = mxv + 1;
          mP->x[yi* mP->nx + (xi-1)] = mxv - 1;
          loop_change++;
        }
        /* make current point below point in previous row */
        if(mP->y[yi* mP->nx + xi] <= mP->y[(yi-1)* mP->nx + xi]) {
          myv = (mP->y[yi* mP->nx + xi] + mP->y[(yi-1)* mP->nx + xi]) / 2;
          mP->y[yi* mP->nx + xi]     = myv + 1;
          mP->y[(yi-1)* mP->nx + xi] = myv - 1 ;
          loop_change++;
        }

        /* make current point inside image boundary */
        if(mP->x[yi* mP->nx + xi] > (mw- mP->nx+xi)) {
          mP->x[yi* mP->nx + xi]  =  mw- mP->nx+xi;
          loop_change++;
        }
        /* make current point inside image boundary */
        if(mP->y[yi* mP->nx + xi] > (mh- mP->ny+yi)) {
          mP->y[yi* mP->nx + xi]  =  mh- mP->ny+yi;
          loop_change++;
        }
      }
    }
    if(loop_change) mesh_change++;
  } while ((mesh_change < (mP->nx+ mP->ny)) && loop_change);

  return mesh_change;
}




/* ---------------------------------------------------------------------- */

/* meshNearestPoint: find the nearest meshpoint and return square of distance
//
// Fill in the indices of the meshpoint and the x,y distances.
// Distances are dx=(px - mP->x[]) , dy=(py - mP->y[])
//
// mP (in): mesh pointer
// px, py: (in) the pointer coordinates (and can be out of range)
// mi, mj: (out) indices of closest meshpoint
// dx, dy: (out) x, y distance of pointer from nearest meshpoint
//
// Returns:  square distance between pointer and meshpoint
*/
int
meshPointNearest(const MeshT *mP, int px, int py, int *mi, int *mj, int *dx, int *dy)
{
  register int      xi, yi;
  register int      m_dx, m_dy;
  register long int m_d;
  long int          m_d_min;

  m_d_min = 15000;

  /* Guarentee p[xy] is in range */
  if(px < mP->x[0]) {
    px = mP->x[0];
  }
  if(py < mP->y[0]) {
    py = mP->y[0];
  }
  if(px > mP->x[mP->ny* mP->nx-1]) {
    px = mP->x[mP->nx* mP->ny-1];
  }
  if(py > mP->y[mP->ny* mP->nx-1]) {
    py = mP->y[mP->nx* mP->ny-1];
  }

  for(yi=0; yi < mP->ny; yi++) {
    for(xi=0; xi < mP->nx; xi++) {
      m_dx = px - mP->x[yi * mP->nx + xi];
      m_dy = py - mP->y[yi * mP->nx + xi];
      m_d = m_dx*m_dx + m_dy*m_dy;

      if(m_d < m_d_min) {
        m_d_min = m_d;
        *mi = xi;
        *mj = yi;
        if(dx!=NULL) *dx = m_dx;
        if(dy!=NULL) *dy = m_dy;
      }
    }
  }

  return m_d_min ;
}




/* ---------------------------------------------------------------------- */

/* meshAddLine: add a mesh line
//
// mP (in/out): mesh pointer
//
// mi (in): upper left index of the quadrangle enclosing the new line
//   (i.e. mi is less than the index of the new line in the new mesh.)
//   For adding vertical lines, mi is the index of the nearby left column.
//   For adding horizontal lines, mi is the index of the nearby upper row.
//
//
// mt (in): relative distance between the surrounding mesh lines
//
// type (in): 1 for vertical lines or 2 for horizontal lines
//
// Note that adding a mesh line has a subtlety concerning
// "location" that is not an issue with deleting mesh lines or
// picking mesh points.  When adding a mesh line, it's not quite so
// simple to figure out where the user is indicating to add the
// line because in general the line could be quite curvy and
// twisted.  This algorithm tries its best, but be careful when
// reading this routine and providing its input arguments.
//
// Allocates memory for the new mesh arrays.
// Sets the incoming mesh array pointer to the newly allocated array.
// Frees the old mesh arrays.
//
// Return zero if okay.
//
// Returns nonzero if fails.
//   Failure can happen if memory runs out or if a bad value for "type"
//   is provided, or if mi is out of bounds.
*/
int
meshAddLine(MeshT *mP, const int mi, const double mt, const int type)
{
  int xi, yi;
  MeshT new;

  meshInitialize(&new, 0, 0);

  /* Set up the new mesh size */
  switch(type) {
    /* Add vertical */
    case 1:
      /* Add a column */
      new.nx = mP->nx + 1;
      new.ny = mP->ny;
      if((mi<0) || (mi > mP->nx)) {
        fprintf(stderr, "meshAddLine: bad value: 0>mi=%i>nx=%i\n", mi, mP->nx);
        return -2;
      }
      break;

    /* Add horizontal */
    case 2:
      /* Add a row */
      new.nx = mP->nx;
      new.ny = mP->ny + 1;
      if((mi<0) || (mi > mP->ny)) {
        fprintf(stderr, "meshAddLine: bad value: 0>mi=%i>ny=%i\n", mi, mP->ny);
        return -3;
      }
      break;

    /* Invalid type */
    default:
      fprintf(stderr, "meshAddLine: Bad Value: type: %i\n", type);
      return -1 ;
  }

  /* Allocate the new mesh */
  if(meshAlloc(&new, "meshAddLine")) return 1 ;

  /* Make the change */
  switch(type) {
    /* --- Add vertical line --- */
    case 1:
      /* Copy the left columns from old into new */
      for(yi=0; yi < mP->ny; yi++) {
        for(xi=0; xi <= mi; xi++) {
          new.x[yi*new.nx + xi] = mP->x[yi * mP->nx + xi];
          new.y[yi*new.nx + xi] = mP->y[yi * mP->nx + xi];
        }
      }

      /* Copy the right columns from old into new */
      for(yi=0; yi < mP->ny; yi++) {
        for(xi=mi+1; xi < mP->nx; xi++) {
          new.x[yi*new.nx + (xi+1)] = mP->x[yi * mP->nx + xi];
          new.y[yi*new.nx + (xi+1)] = mP->y[yi * mP->nx + xi];
        }
      }

      /* Add the new column */
      {
        double mx1, mx2, mxv;
        double my1, my2, myv;
        for(yi=0; yi < mP->ny; yi++) {
          /* Place new line between two horizontally adjacent lines */
          mx1 = mP->x[yi* mP->nx + mi];
          mx2 = mP->x[yi* mP->nx + (mi+1)];
          mxv = (1.0-mt)*mx1 + mt*mx2;
          new.x[yi*new.nx + (mi+1)] = mxv;

          my1 = mP->y[yi* mP->nx + mi];
          my2 = mP->y[yi* mP->nx + (mi+1)];
          myv = (1.0-mt)*my1 + mt*my2;
          new.y[yi*new.nx + (mi+1)] = myv;
        }
      }
      break;

    /* --- Add horizontal line --- */
    case 2:
      /* Copy the top rows from old to new */
      for(yi=0; yi <= mi; yi++) {
        for(xi=0; xi< mP->nx; xi++) {
          new.x[yi*new.nx + xi] = mP->x[yi * mP->nx + xi];
          new.y[yi*new.nx + xi] = mP->y[yi * mP->nx + xi];
        }
      }

      /* Copy the bottom rows from old to new */
      for(yi=mi+1; yi < mP->ny; yi++) {
        for(xi=0; xi < mP->nx; xi++) {
          new.x[(yi+1)*new.nx + xi] = mP->x[yi * mP->nx + xi];
          new.y[(yi+1)*new.nx + xi] = mP->y[yi * mP->nx + xi];
        }
      }

      /* Add the new row */
      {
        double mx1, mx2, mxv;
        double my1, my2, myv;
        for(xi=0; xi < mP->nx; xi++) {
          /* Place new line between two vertically adjacent lines */
          mx1 = mP->x[(mi)* mP->nx + xi];
          mx2 = mP->x[(mi+1)* mP->nx + xi];
          mxv = (1.0-mt)*mx1 + mt*mx2;
          new.x[(mi+1)*new.nx + xi] = mxv;

          my1 = mP->y[(mi)* mP->nx + xi];
          my2 = mP->y[(mi+1)* mP->nx + xi];
          myv = (1.0-mt)*my1 + mt*my2;
          new.y[(mi+1)*new.nx + xi] = myv;
        }
      }
      break;

    /* --- Invalid type --- */
    default:
      fprintf(stderr, "meshAddLine: Bad Value: type: %i\n", type);
      return -1 ;
  }

  /* Free the old mesh arrays */
  meshFree(mP, "meshAddLine");

  /* Point to the new mesh arrays */
  mP->x  = new.x;
  mP->y  = new.y;
  mP->nx = new.nx;
  mP->ny = new.ny;

  return 0 ;
}




/* ---------------------------------------------------------------------- */

/* meshDeleteLine: delete a mesh line
//
// mP: mesh pointer
//
// mi is index the of a mesh point on the to-be-deleted line
//  for deleting vertical lines, mi is the index of the column
//  for deleting horizontal lines, mi is the index of the row
//
// type is 1 for vertical lines or 2 for horizontal lines
//
// Allocate memory for the mesh and set the incoming mesh pointers to
//  the newly allocated array
//
// Caller must decrement the appropriate local analogy to nx or ny
//
// Return nonzero if meshDeleteLine fails
*/
int
meshDeleteLine(MeshT *mP, int mi, int type)
{
  int xi, yi;
  MeshT new;

  meshInitialize(&new, 0, 0);

  switch(type) {
    /* Delete vertical */
    case 1:
      /* Delete a column */
      new.nx = mP->nx - 1;
      new.ny = mP->ny;
      break;

    /* Delete horizontal */
    case 2:
      /* Delete a row */
      new.nx = mP->nx;
      new.ny = mP->ny - 1;
      break;

    /* Invalid type */
    default:
      fprintf(stderr, "meshDeleteLine: Bad Value: type: %i\n", type);
      return -1;
  }

  if(meshAlloc(&new, "meshDeleteLine")) return 1;

  switch(type) {
    /* --- Delete vertical line --- */
    case 1:
      /* Copy the left columns */
      for(yi=0; yi < mP->ny; yi++) {
        for(xi=0; xi<mi; xi++) {
          new.x[yi*new.nx + xi] = mP->x[yi* mP->nx + xi];
          new.y[yi*new.nx + xi] = mP->y[yi* mP->nx + xi];
        }
      }

      /* Copy the right columns */
      for(yi=0; yi< mP->ny; yi++) {
        for(xi=mi+1; xi< mP->nx; xi++) {
          new.x[yi*new.nx + (xi-1)] = mP->x[yi* mP->nx + xi];
          new.y[yi*new.nx + (xi-1)] = mP->y[yi* mP->nx + xi];
        }
      }

      break;

    /* --- Delete horizontal line --- */
    case 2:
      /* Copy the top rows */
      for(yi=0; yi<mi; yi++) {
        for(xi=0; xi< mP->nx; xi++) {
          new.x[yi*new.nx + xi] = mP->x[yi* mP->nx + xi];
          new.y[yi*new.nx + xi] = mP->y[yi* mP->nx + xi];
        }
      }

      /* Copy the bottom rows */
      for(yi=mi+1; yi< mP->ny; yi++) {
        for(xi=0; xi< mP->nx; xi++) {
          new.x[(yi-1)*new.nx + xi] = mP->x[yi* mP->nx + xi];
          new.y[(yi-1)*new.nx + xi] = mP->y[yi* mP->nx + xi];
        }
      }

      break;

    /* --- --- --- Invalid type --- --- --- */
    default:
      fprintf(stderr, "meshDeleteLine: Bad Value: type: %i\n", type);
      return -1;
  }

  /* Free the old mesh arrays */
  meshFree(mP, "meshDeleteLine");

  /* Point to the new mesh arrays */
  mP->x = new.x;
  mP->y = new.y;
  mP->nx = new.nx;
  mP->ny = new.ny;

  return 0;
}




/* ---------------------------------------------------------------------- */

/* meshRead: Read a mesh from a file
//
// meshP: (in/out) mesh pointer
// filename: (in) mesh file name
//
// Frees memory of previous mesh.
// Allocates memory for the meshes and sets .nx and .ny members.
//
// Returns zero if load succeeds.
// Returns nonzero if load fails.
*/
int
meshRead(MeshT *meshP, const char *filename)
{
  int xi, yi;        /* loop mesh indices */
  double mesh_point; /* input buffer for mesh point value */
  char magic[2];     /* magic number, for file identification */
  FILE *fP;          /* mesh file pointer */

  /* Open mesh file for reading */
  if((fP=fopen(filename, "r"))==NULL) {
    fprintf(stderr, "meshRead: could not read file '%s'\n", filename);
    return 1;
  }

  /* Read first two characters of mesh file */
  if(get_block(fP, magic, 2)) {
    fprintf(stderr, "meshRead: premature EOF in file '%s'\n", filename);
    fclose(fP);
    return EOF;
  }

  /* "M2" as the first two characters indicates an ASCII mesh file */
  if(magic[0]=='M' && magic[1]=='2') {
    /* Read the mesh geometry */
    if(fscanf(fP, "%li", &(meshP->nx))!=1 || (meshP->nx < 0)) {
      fprintf(stderr, "meshRead: missing nx\n");
      fclose(fP);
      return 2;
    }
    if(fscanf(fP, "%li", &(meshP->ny))!=1 || (meshP->ny < 0)) {
      fprintf(stderr, "meshRead: missing ny\n");
      fclose(fP);
      return 3;
    }

    /* Free the old mesh and allocate memory for the new mesh */
    meshFree(meshP, "meshRead");
    if(meshAlloc(meshP, "meshRead")) {
      fclose(fP);
      return 6;
    }

    /* Read the mesh point values */
#ifdef TRANSPOSE_MESH
    for(xi=0; xi < meshP->nx; xi++) { /* brace for matching } */
      for(yi=0; yi < meshP->ny; yi++) { /* brace for matching } */
#else
    for(yi=0; yi < meshP->ny; yi++) {
      for(xi=0; xi < meshP->nx; xi++) {
#endif
        if(fscanf(fP, "%lf", &mesh_point)!=1) {
          fprintf(stderr, "meshRead: missing float at %i %i\n", xi, yi);
          fclose(fP);
          meshFree(meshP, "meshRead");
          return 4;
        }
        meshP->x[yi* meshP->nx + xi] = mesh_point;

        if(fscanf(fP, "%lf", &mesh_point)!=1) {
          fprintf(stderr, "meshRead: missing float at %i %i\n", xi, yi);
          fclose(fP);
          meshFree(meshP, "meshRead");
          return 4;
        }
        meshP->y[yi* meshP->nx + xi] = mesh_point;
      }
    }

  } else {
    fprintf(stderr, "meshRead: file was not a valid mesh file\n");
    fclose(fP);
    return 5;
  }

  fclose(fP);
  return 0;
}




/* ---------------------------------------------------------------------- */

/* meshWrite: save a mesh to file named filename
//
// meshP: (in) mesh pointer
// filename: (in) mesh file name
//
// Returns zero if load succeeds.
// Returns nonzero if load fails.
*/
int
meshWrite(const MeshT *meshP, char *filename)
{
  int xi, yi;
  FILE *fP;

  if((fP=fopen(filename, "w"))==NULL) {
    fprintf(stderr, "meshWrite: could not write file '%s \n", filename);
    return 1;
  }

  /* M2 indicates an ASCII mesh file */
  fprintf(fP, "M2\n");

  /* Write the mesh geometry */
  fprintf(fP, "%i %i\n", meshP->nx, meshP->ny);

  /* Write the mesh point values */
  for(yi=0; yi < meshP->ny; yi++) {
    for(xi=0; xi < meshP->nx; xi++) {
      fprintf(fP, "%lg ", meshP->x[yi * meshP->nx + xi]);
      fprintf(fP, "%lg\n", meshP->y[yi * meshP->nx + xi]);
    }
  }
  fclose(fP);
  return 0;
}




/* ---------------------------------------------------------------------- */

/* meshMatch: Make mesh dimensions match another mesh
//
// If this mesh needs to be resized, then it is also reset.
// This should probably be rewritten to simply rescale the other mesh,
// instead of destroying it.
//
// this: (in/out) mesh to be made to match the other mesh
// other: (in) mesh to be matched to
*/
void
meshMatch(MeshT *this, const MeshT *other)
{
  if((this->nx != other->nx) || (this->ny != other->ny)) {
    this->nx = other->nx;
    this->ny = other->ny;

    meshFree(this, "meshMatch");
    meshAlloc(this, "meshMatch");

    /* The +1 is because meshReset wants width and height , but the
    // mesh stores pixel coordinates, which range from 0 to size-1 in
    // each direction.  The +0.5 is to avoid roundoff truncation.
    */
    meshReset(this,
      other->x[other->nx * other->ny - 1] + 1.5,
      other->y[other->nx * other->ny - 1] + 1.5);
  }
}