# include "data.h"

int minterp (int argc, char **argv) {
  
  int i, j, status, nx, ny, Nx, Ny;
  char temp[1024];
  double scale, scale2, dX, dY;
  float *V00, *V01, *V10, *V11, *Vout, dV1, dV2, dV3;
  Buffer *in, *out;

  if (argc != 4) {
    fprintf (stderr, "USAGE: minterpolate <from> <to> scale \n");
    return (FALSE);
  }

  if ((in  = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
  if ((out = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);

  fits_free_matrix (&out[0].matrix);
  fits_free_header (&out[0].header);

  scale  = atof (argv[3]);
  scale2 = scale*scale;
  Nx = in[0].header.Naxis[0];
  Ny = in[0].header.Naxis[1];
  nx = Nx * scale;
  ny = Ny * scale;

  /* create new matrix */
  out[0].bitpix = in[0].bitpix;
  out[0].unsign = in[0].unsign;
  out[0].bscale = in[0].bscale;
  out[0].bzero  = in[0].bzero;
  fits_copy_header (&in[0].header, &out[0].header);

  fits_modify (&out[0].header, "NAXIS1", "%d", 1, nx);
  fits_modify (&out[0].header, "NAXIS2", "%d", 1, ny);
  out[0].header.Naxis[0] = nx;
  out[0].header.Naxis[1] = ny;
  fits_create_matrix (&out[0].header, &out[0].matrix);

  /* fix astrometric terms */
  status =  fits_scan (&out[0].header, "CDELT1", "%lf", 1, &dX);
  status &= fits_scan (&out[0].header, "CDELT2", "%lf", 1, &dY);
  dX /= scale;
  dY /= scale;
  if (status) {
    fits_modify (&out[0].header, "CDELT1", "%lf", 1, dX);
    fits_modify (&out[0].header, "CDELT2", "%lf", 1, dY);
  }
  status =  fits_scan (&out[0].header, "CRPIX1", "%lf", 1, &dX);
  status &= fits_scan (&out[0].header, "CRPIX2", "%lf", 1, &dY);
  dX *= scale;
  dY *= scale;
  if (status) {
    fits_modify (&out[0].header, "CRPIX1", "%lf", 1, dX);
    fits_modify (&out[0].header, "CRPIX2", "%lf", 1, dY);
  }

  /* adjust filename */
  temp[0] = 0;
  if ((in[0].file[0] != '*') && (in[0].file[0] != '(')) {
    strcpy (temp, "*");
  }
  strcat (temp, in[0].file);
  strcpy (out[0].file, temp);

  dX = dY = scale;

  {

    int dx, dy, Npix;
    float I, J, x, y, xs, xe, ys, ye;

    Npix = 0;
    for (j = 0; j < Ny - 1; j++) {
      for (i = 0; i < Nx - 1; i++) {
	V00 = (float *)(in[0].matrix.buffer) + j*Nx + i;
	V10 = V00 + 1;
	V01 = V00 + Nx;
	V11 = V01 + 1;
	dV1 = (*V11 + *V00 - *V01 - *V10) / scale2;
	dV2 = (*V01 - *V00) / scale;
	dV3 = (*V10 - *V00) / scale;

	x = (i + 0.5) * scale;
	y = (j + 0.5) * scale;

	xs = 0;
	if (i == 0) { xs = -0.5*scale; }
	xe = scale;
	if (i == Nx - 2) { xe = 1.5*scale; }

	ys = 0;
	if (j == 0) { ys = -0.5*scale; }
	ye = scale;
	if (j == Ny - 2) { ye = 1.5*scale; }

	for (J = ys; J < ye; J += 1.0) {
	  dx = (x + xs);
	  dy = (y + J);
	  Vout = (float *)(out[0].matrix.buffer) + dy*nx + dx;
	  for (I = xs; I < xe; I += 1.0, Vout++) {
	    *Vout = dV1 * (I*J) + dV2 * J + dV3 * I + *V00;
	    Npix ++;
	  }
	}
      }
    }
  }

  return (TRUE);

}

