Index: /trunk/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 16106)
+++ /trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 16107)
@@ -44,4 +44,5 @@
 $(SRC)/fft1d.old.$(ARCH).o		\
 $(SRC)/fft1d.$(ARCH).o		\
+$(SRC)/fft2d.old.$(ARCH).o		\
 $(SRC)/fft2d.$(ARCH).o		\
 $(SRC)/fit2d.$(ARCH).o		\
Index: /trunk/Ohana/src/opihi/cmd.data/fft1d.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/fft1d.c	(revision 16106)
+++ /trunk/Ohana/src/opihi/cmd.data/fft1d.c	(revision 16107)
@@ -3,6 +3,12 @@
 int fft1d (int argc, char **argv) {
   
-  int Npix, Nbit, ZeroImaginary;
+  int N, Npix, Nbit, ZeroImaginary, forward;
   Vector *Ire, *Iim, *Ore, *Oim;
+
+  forward = TRUE;
+  if ((N = get_argument (argc, argv, "-inverse"))) {
+    remove_argument (N, &argc, argv);
+    forward = FALSE;
+  }
 
   if ((argc != 6) || (strcmp (argv[3], "to"))) {
@@ -22,5 +28,5 @@
   Npix = Ire[0].Nelements;
   if (!ZeroImaginary && (Npix != Iim[0].Nelements)) {
-    gprint (GP_ERR, "vector mismatch in size\n");
+    gprint (GP_ERR, "vector size mismatch\n");
     return (FALSE);
   }
@@ -50,5 +56,5 @@
   }    
 
-  fft (Ore[0].elements, Oim[0].elements, Npix, Nbit, TRUE); 
+  fft1D (Ore[0].elements, Oim[0].elements, Npix, Nbit, forward); 
   
   return (TRUE);
Index: /trunk/Ohana/src/opihi/cmd.data/fft2d.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/fft2d.c	(revision 16106)
+++ /trunk/Ohana/src/opihi/cmd.data/fft2d.c	(revision 16107)
@@ -3,13 +3,12 @@
 int fft2d (int argc, char **argv) {
   
-  int i, N, Nx, Ny, Naxis[2];
-  int Npix, ZeroImaginary, isign;
-  float *t1, *t2, *out, *temp;
+  int ZeroImaginary, forward;
+  int N, Nx, Ny;
   Buffer *Ire, *Iim, *Ore, *Oim;;
 
-  isign = 1;
+  forward = TRUE;
   if ((N = get_argument (argc, argv, "-inverse"))) {
     remove_argument (N, &argc, argv);
-    isign = -1;
+    forward = FALSE;
   }
 
@@ -28,4 +27,25 @@
   }    
   if ((Ire = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+
+  Nx = Ire[0].header.Naxis[0];
+  Ny = Ire[0].header.Naxis[1];
+
+  // check input image dimensions (match lengths? binary lengths?)
+  if (!ZeroImaginary) {
+    if ((Nx != Iim[0].header.Naxis[0]) || 
+	(Ny != Iim[0].header.Naxis[1])) {
+      gprint (GP_ERR, "image size mismatch\n");
+      return (FALSE);
+    }
+  }
+  if (!IsBinary (Nx, NULL)) {
+    gprint (GP_ERR, "Nx is not a binary number!\n");
+    return (FALSE);
+  }
+  if (!IsBinary (Ny, NULL)) {
+    gprint (GP_ERR, "Ny is not a binary number!\n");
+    return (FALSE);
+  }
+
   if ((Ore = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
   if ((Oim = SelectBuffer (argv[5], ANYBUFFER, TRUE)) == NULL) return (FALSE);
@@ -36,77 +56,23 @@
   gfits_free_matrix (&Oim[0].matrix);
   gfits_free_header (&Oim[0].header);
+  
+  /* fix up output headers (real) & allocate data buffer */
+  CreateBuffer (Ore, Nx, Ny, -32, 0.0, 1.0);
+  CreateBuffer (Oim, Nx, Ny, -32, 0.0, 1.0);
 
-  /* get image dimensions, check value */
-  Npix = Ire[0].header.Naxis[0]*Ire[0].header.Naxis[1];
-  Nx = Ire[0].header.Naxis[0];
-  Ny = Ire[0].header.Naxis[1];
-  Naxis[0] = Ny; Naxis[1] = Nx;
-  if (!IsBinaryOld (Npix)) {
-    gprint (GP_ERR, "dimensions are not binary!\n");
-    return (FALSE);
+  gfits_copy_header (&Ire[0].header, &Ore[0].header);
+  gfits_copy_header (&Ire[0].header, &Oim[0].header);
+
+  // copy data to output buffers (fft is done in place)
+  memcpy (Ore[0].matrix.buffer, Ire[0].matrix.buffer, Nx*Ny*sizeof(float));
+
+  if (ZeroImaginary) {
+    memset (Oim[0].matrix.buffer, 0, Nx*Ny*sizeof(float));
+  } else {
+    memcpy (Oim[0].matrix.buffer, Iim[0].matrix.buffer, Nx*Ny*sizeof(float));
   }
-  
-  /* create working space */
-  t1 = (float *) Ire[0].matrix.buffer;
-  ALLOCATE (temp, float, 2*Npix);
-  out = temp;
 
-  /* copy input to working space */
-  if (ZeroImaginary) {
-    for (i = 0; i < Npix; i++, t1++) {
-      *out = *t1;
-      out++;
-      *out = 0;
-      out++;
-    }
-  } else {
-    t2 = (float *) Iim[0].matrix.buffer;
-    for (i = 0; i < Npix; i++, t1++, t2++) {
-      *out = *t1;
-      out++;
-      *out = *t2;
-      out++;
-    }
-  } 
-    
   /* run the fft */
-  fftN (temp, Naxis, 2, isign);
-
-  /* fix up output headers (real) */
-  gfits_copy_header (&Ire[0].header, &Ore[0].header);
-  gfits_modify (&Ore[0].header, "NAXIS1", "%d", 1, Nx);
-  gfits_modify (&Ore[0].header, "NAXIS2", "%d", 1, Ny);
-  Ore[0].header.Naxis[0] = Nx;
-  Ore[0].header.Naxis[1] = Ny;
-  Ore[0].bitpix = Ire[0].bitpix;
-  Ore[0].unsign = Ire[0].unsign;
-  Ore[0].bscale = Ire[0].bscale;
-  Ore[0].bzero  = Ire[0].bzero;
-  gfits_create_matrix (&Ore[0].header, &Ore[0].matrix);
-
-  /* fix up output headers (imaginary) */
-  gfits_copy_header (&Ire[0].header, &Oim[0].header);
-  gfits_modify (&Oim[0].header, "NAXIS1", "%d", 1, Nx);
-  gfits_modify (&Oim[0].header, "NAXIS2", "%d", 1, Ny);
-  Oim[0].header.Naxis[0] = Nx;
-  Oim[0].header.Naxis[1] = Ny;
-  Oim[0].bitpix = Ire[0].bitpix;
-  Oim[0].unsign = Ire[0].unsign;
-  Oim[0].bscale = Ire[0].bscale;
-  Oim[0].bzero  = Ire[0].bzero;
-  gfits_create_matrix (&Oim[0].header, &Oim[0].matrix);
-
-  /* move data from working space to output buffers */
-  out = temp;
-  t1 = (float *) Ore[0].matrix.buffer;
-  t2 = (float *) Oim[0].matrix.buffer;
-  for (i = 0; i < Npix; i++, t1++, t2++) {
-    *t1 = *out / Npix;
-    out ++;
-    *t2 = *out / Npix;
-    out ++;
-  }    
-
-  free (temp);
+  fftND ((float *)Ore[0].matrix.buffer, (float *)Oim[0].matrix.buffer, 2, Ire[0].header.Naxis, forward);
 
   return (TRUE);
Index: /trunk/Ohana/src/opihi/cmd.data/fft2d.old.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/fft2d.old.c	(revision 16107)
+++ /trunk/Ohana/src/opihi/cmd.data/fft2d.old.c	(revision 16107)
@@ -0,0 +1,113 @@
+# include "data.h"
+
+int fft2dold (int argc, char **argv) {
+  
+  int i, N, Nx, Ny, Naxis[2];
+  int Npix, ZeroImaginary, isign;
+  float *t1, *t2, *out, *temp;
+  Buffer *Ire, *Iim, *Ore, *Oim;;
+
+  isign = 1;
+  if ((N = get_argument (argc, argv, "-inverse"))) {
+    remove_argument (N, &argc, argv);
+    isign = -1;
+  }
+
+  if ((argc != 6) || (strcmp (argv[3], "to"))) {
+    gprint (GP_ERR, "USAGE: fft2d (real) (imag) to (real) (imag)\n");
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  Iim = NULL;
+  ZeroImaginary = TRUE; /* Input(imaginary) may be 0, in which case we create a 0 filled image */
+  if (!strcmp (argv[2], "0")) { 
+  } else {
+    ZeroImaginary = FALSE;
+    if ((Iim = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  }    
+  if ((Ire = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((Ore = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((Oim = SelectBuffer (argv[5], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
+  /* free up output space */
+  gfits_free_matrix (&Ore[0].matrix);
+  gfits_free_header (&Ore[0].header);
+  gfits_free_matrix (&Oim[0].matrix);
+  gfits_free_header (&Oim[0].header);
+
+  /* get image dimensions, check value */
+  Npix = Ire[0].header.Naxis[0]*Ire[0].header.Naxis[1];
+  Nx = Ire[0].header.Naxis[0];
+  Ny = Ire[0].header.Naxis[1];
+  Naxis[0] = Ny; Naxis[1] = Nx;
+  if (!IsBinaryOld (Npix)) {
+    gprint (GP_ERR, "dimensions are not binary!\n");
+    return (FALSE);
+  }
+  
+  /* create working space */
+  t1 = (float *) Ire[0].matrix.buffer;
+  ALLOCATE (temp, float, 2*Npix);
+  out = temp;
+
+  /* copy input to working space */
+  if (ZeroImaginary) {
+    for (i = 0; i < Npix; i++, t1++) {
+      *out = *t1;
+      out++;
+      *out = 0;
+      out++;
+    }
+  } else {
+    t2 = (float *) Iim[0].matrix.buffer;
+    for (i = 0; i < Npix; i++, t1++, t2++) {
+      *out = *t1;
+      out++;
+      *out = *t2;
+      out++;
+    }
+  } 
+    
+  /* run the fft */
+  fftNold (temp, Naxis, 2, isign);
+
+  /* fix up output headers (real) */
+  gfits_copy_header (&Ire[0].header, &Ore[0].header);
+  gfits_modify (&Ore[0].header, "NAXIS1", "%d", 1, Nx);
+  gfits_modify (&Ore[0].header, "NAXIS2", "%d", 1, Ny);
+  Ore[0].header.Naxis[0] = Nx;
+  Ore[0].header.Naxis[1] = Ny;
+  Ore[0].bitpix = Ire[0].bitpix;
+  Ore[0].unsign = Ire[0].unsign;
+  Ore[0].bscale = Ire[0].bscale;
+  Ore[0].bzero  = Ire[0].bzero;
+  gfits_create_matrix (&Ore[0].header, &Ore[0].matrix);
+
+  /* fix up output headers (imaginary) */
+  gfits_copy_header (&Ire[0].header, &Oim[0].header);
+  gfits_modify (&Oim[0].header, "NAXIS1", "%d", 1, Nx);
+  gfits_modify (&Oim[0].header, "NAXIS2", "%d", 1, Ny);
+  Oim[0].header.Naxis[0] = Nx;
+  Oim[0].header.Naxis[1] = Ny;
+  Oim[0].bitpix = Ire[0].bitpix;
+  Oim[0].unsign = Ire[0].unsign;
+  Oim[0].bscale = Ire[0].bscale;
+  Oim[0].bzero  = Ire[0].bzero;
+  gfits_create_matrix (&Oim[0].header, &Oim[0].matrix);
+
+  /* move data from working space to output buffers */
+  out = temp;
+  t1 = (float *) Ore[0].matrix.buffer;
+  t2 = (float *) Oim[0].matrix.buffer;
+  for (i = 0; i < Npix; i++, t1++, t2++) {
+    *t1 = *out / Npix;
+    out ++;
+    *t2 = *out / Npix;
+    out ++;
+  }    
+
+  free (temp);
+
+  return (TRUE);
+}
Index: /trunk/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/init.c	(revision 16106)
+++ /trunk/Ohana/src/opihi/cmd.data/init.c	(revision 16107)
@@ -29,4 +29,5 @@
 int fft1d            PROTO((int, char **));
 int fft2d            PROTO((int, char **));
+int fft2dold         PROTO((int, char **));
 int fit2d            PROTO((int, char **));
 int fit              PROTO((int, char **));
@@ -149,4 +150,5 @@
   {"fft1d",        fft1d,	     "fft on the pixel-stream in an image"},
   {"fft2d",        fft2d,	     "fft on an image"},
+  {"fft2dold",     fft2dold,	     "fft on an image"},
   {"fit",     	   fit,		     "fit polynomial to vector pair"},
   {"fit2d",        fit2d,	     "fit 2-d polynomial to vector triplet"},
Index: /trunk/Ohana/src/opihi/cmd.data/test/fft2d.sh
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/test/fft2d.sh	(revision 16107)
+++ /trunk/Ohana/src/opihi/cmd.data/test/fft2d.sh	(revision 16107)
@@ -0,0 +1,24 @@
+
+list tests
+ test1
+ test2
+ test3
+end
+
+macro test1
+ $PASS = 1
+ break -auto off
+
+ mcreate t 2048 2048
+ set x = xramp(t)
+ set y = yramp(t)
+
+ # set f = dsin(3*x)*dcos(5*y)
+ set f = exp(-0.5*((x-1024)^2 + (y-1024)^2)*0.01)
+
+ # date; fft2d f 0 to Frn Fin; date
+
+ # date; fft2dold f 0 to Fro Fio; date
+
+ tv f -1 1
+end
Index: /trunk/Ohana/src/opihi/include/data.h
===================================================================
--- /trunk/Ohana/src/opihi/include/data.h	(revision 16106)
+++ /trunk/Ohana/src/opihi/include/data.h	(revision 16107)
@@ -84,6 +84,6 @@
 
 /* in fft.c */
-void fft (float *dataRe, float *dataIm, int N, int Nbit, int forward);
-void fftN (float *data, int *nn, int ndim, int isign);
+void fft1D (float *dataRe, float *dataIm, int N, int Nbit, int forward);
+int fftND (float *dataRe, float *dataIm, int Ndim, int *Nsize, int forward);
 int IsBinary (int N, int *Nbit);
 
Index: /trunk/Ohana/src/opihi/lib.data/fft.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.data/fft.c	(revision 16106)
+++ /trunk/Ohana/src/opihi/lib.data/fft.c	(revision 16107)
@@ -2,11 +2,5 @@
 
 // fft based on code by Douglas L. Jones (see note at EOF). modified for Ohana C style
-void fftN (float *data, int *nn, int ndim, int isign) {
-  
-  fprintf (stderr, "broken\n");
-  return;
-}
-
-void fft (float *x, float *y, int n, int Nbit, int forward) {
+void fft1D (float *x, float *y, int n, int Nbit, int forward) {
 
   int i,j,k,n1,n2;
@@ -75,4 +69,67 @@
 }                          
 
+// This implementation uses the 1-D fft above for each of the vectors in each dimension.
+// This requires 2(Nx*Ny*...) mem copies, but the fft operations are likely to happen in
+// cache.
+int fftND (float *x, float *y, int Ndim, int *Nsize, int forward) {
+
+  int i, nIndex, minor, major, iDim;
+  int step, Nmajor, Nminor, Nmax, Ntotal;
+  int *Nbit;
+  float *tmpX, *tmpY;
+
+  ALLOCATE (Nbit, int, Ndim);
+
+  // find the longest axis and allocate storage for that length
+  Nmax = 0;
+  Ntotal = 1;
+  for (i = 0; i < Ndim; i++) {
+    Nmax = MAX(Nmax, Nsize[i]);
+    Ntotal *= Nsize[i];
+    if (!IsBinary (Nsize[i], &Nbit[i])) {
+      free (Nbit);
+      return (FALSE);
+    }
+  }
+  ALLOCATE (tmpX, float, Nmax);
+  ALLOCATE (tmpY, float, Nmax);
+  
+  step = 1;
+  Nminor = 1;
+  Nmajor = Ntotal;
+  for (iDim = 0; iDim < Ndim; iDim++) {
+    step *= Nsize[iDim];
+    Nmajor /= Nsize[iDim];
+
+    // we perform the FFT along all other dimensions 
+    for (major = 0; major < Nmajor; major++) {
+      for (minor = 0; minor < Nminor; minor++) {
+	// nIndex = minor + i*Nminor + major*step;
+	// extract the data values to the temp vector
+	nIndex = minor + major*step;
+	for (i = 0; i < Nsize[iDim]; i++) {
+	  tmpX[i] = x[nIndex];
+	  tmpY[i] = y[nIndex];
+	  nIndex += Nminor;
+	}
+
+	fft1D (tmpX, tmpY, Nsize[iDim], Nbit[iDim], forward);
+
+	// replace the result vectors
+	nIndex = minor + major*step;
+	for (i = 0; i < Nsize[iDim]; i++) {
+	  x[nIndex] = tmpX[i];
+	  y[nIndex] = tmpY[i];
+	  nIndex += Nminor;
+	}
+      }
+    }
+    Nminor *= Nsize[iDim];
+  }
+  free (tmpX);
+  free (tmpY);
+  return (TRUE);
+}
+
 // check that a number is binary (2^Nbit).  returns int(log_2(N)) in Nbit
 int IsBinary (int N, int *Nbit) {
@@ -80,14 +137,15 @@
   int i, Nset;
 
-  *Nbit = 0;
+  if (Nbit != NULL) *Nbit = 0;
   Nset = 0;
   for (i = 0; i < 8*sizeof(int); i++) {
     if (N & 0x01) {
       Nset ++;
-      *Nbit = i;
+      if (Nbit != NULL) *Nbit = i;
     }
     N >>= 1;
   }
   if (Nset > 1) return (FALSE);
+  return (TRUE);  
 }
 
Index: /trunk/Ohana/src/opihi/lib.data/fftold.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.data/fftold.c	(revision 16106)
+++ /trunk/Ohana/src/opihi/lib.data/fftold.c	(revision 16107)
@@ -49,58 +49,4 @@
     mmax = istep;
   }
-}
-
-void fftold_Nindices (float *Data, int N, int isign) {
-
-  int n,ifp1,m,j,ipf2,i;
-  double wtemp,wr,wpr,wpi,wi,theta;
-  float tempr,tempi, *data;
-
-  data = Data;
-  n = N << 1;
-
-  ip1 = 2;
-  ip2 = n;
-  i2rev = 0;
-  for (i2 = 0; i2 < ip2; i2 += ip1) {
-    if (i2 < i2ref) {
-      FSWAP (data[i2rev], data[i2]);
-      FSWAP (data[i2rev+1], data[i2+1]);
-    }
-    ibit = ip2 >> 1;
-    while (ibit >= ip1 && i2rev >= ibit) {
-      i2rev -= ibit;
-      ibit >>= 1;
-    }
-    i2rev += ibit;
-  }
-
-  ifp1 = 2;
-  while (ifp1 < ip2) {
-    ifp2 = ifp1 << 1;
-    theta = isign*6.28318530717959 / (2*ifp1/ip1);
-    wtemp = sin(0.5*theta);
-    wpr = -2.0*wtemp*wtemp;
-    wpi = sin(theta);
-    wr = 1.0;
-    wi = 0.0;
-    for (i3 = 0; i3 < ifp1; i3+=ip1) {
-      // in the N-d version, this loop is broken into two loops
-      for (i1 = i3; i1 < n; i1+=ipf2) {
-	k1 = i1;
-	k2 = k1 + ifp1;
-	tempr = wr*data[k2] - wi*data[k2+1];
-	tempi = wr*data[k2+1] + wi*data[k2];
-	data[k2] = data[k1] - tempr;
-	data[k2+1] = data[k1+1] - tempi;
-	data[k1] += tempr;
-	data[k1+1] += tempi;
-      }
-      wr = (wtemp = wr)*wpr - wi*wpi+wr;
-      wi = wi*wpr + wtemp*wpi + wi;
-    }
-    ifp1 = ifp2;
-  }
-
 }