Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile	(revision 37732)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/Makefile	(revision 37733)
@@ -39,4 +39,5 @@
 $(SRC)/densify.$(ARCH).o	\
 $(SRC)/device.$(ARCH).o	\
+$(SRC)/dft2d.$(ARCH).o	\
 $(SRC)/dimendown.$(ARCH).o	\
 $(SRC)/dimenup.$(ARCH).o	\
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/dft2d.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/dft2d.c	(revision 37733)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/dft2d.c	(revision 37733)
@@ -0,0 +1,142 @@
+# include "data.h"
+
+// perform a 2D fourier transform.  
+
+static int NxLast = 0;
+static int NyLast = 0;
+static float *Fterms = NULL;
+
+int dft2d (int argc, char **argv) {
+  
+  int ix, iy, wx, wy;
+  Buffer *src = NULL;
+  Buffer *tgt = NULL;
+
+  if (argc != 4) goto usage;
+  
+  if (strcasecmp(argv[2], "to")) goto usage;
+  if ((src = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) goto usage;
+  if ((tgt = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) goto usage;
+
+  int Nx = src[0].matrix.Naxis[0];
+  int Ny = src[0].matrix.Naxis[1];
+  ResetBuffer (tgt, Nx, Ny, -32, 0.0, 1.0);
+
+  if (Nx % 2) goto usage;
+  if (Ny % 2) goto usage;
+
+  // generate and save the fourier terms, if needed
+  if (Fterms && (NxLast == Nx) && (NyLast == Ny)) goto reuse;
+
+  if (!Fterms) {
+    ALLOCATE (Fterms, float, SQ(Nx*Ny));
+    NxLast = Nx;
+    NyLast = Ny;
+  } else {
+    REALLOCATE (Fterms, float, SQ(Nx*Ny));
+    NxLast = Nx;
+    NyLast = Ny;
+  }
+
+  // x-dir cosine terms
+  float Wx = 2*M_PI/(float)Nx;
+  float Wy = 2*M_PI/(float)Ny;
+  for (wx = 0; wx <= Nx / 2; wx++) {
+    for (wy = 0; wy <= Ny / 2; wy++) {
+      float rx = ((wx == 0) || (wx == Nx/2)) ? (1.0/Nx) : (2.0/Nx);
+      float ry = ((wy == 0) || (wy == Ny/2)) ? (1.0/Ny) : (2.0/Ny);
+      float rn = rx*ry;
+      for (ix = 0; ix < Nx; ix++) {
+	float fxc = cos(wx*Wx*ix);
+	for (iy = 0; iy < Ny; iy++) {
+	  float fyc = cos(wy*Wy*iy);
+	  int Nf = wx + wy*Ny;
+	  int out = ix + iy*Nx + Nf*Nx*Ny;
+	  Fterms[out] = fxc*fyc*rn;
+	  if (wy == 0) continue;
+	  if (wy == Ny/2) continue;
+	  float fys = sin(wy*Wy*iy);
+	  Nf = wx + (Ny - wy)*Ny;
+	  out = ix + iy*Nx + Nf*Nx*Ny;
+	  Fterms[out] = fxc*fys*rn;
+	}
+	if (wx == 0) continue;
+	if (wx == Nx/2) continue;
+	float fxs = sin(wx*Wx*ix);
+	for (iy = 0; iy < Ny; iy++) {
+	  float fyc = cos(wy*Wy*iy);
+	  int Nf = (Nx - wx) + wy*Ny;
+	  int out = ix + iy*Nx + Nf*Nx*Ny;
+	  Fterms[out] = fxs*fyc*rn;
+	  if (wy == 0) continue;
+	  if (wy == Ny/2) continue;
+	  float fys = sin(wy*Wy*iy);
+	  Nf = (Nx - wx) + (Ny - wy)*Ny;
+	  out = ix + iy*Nx + Nf*Nx*Ny;
+	  Fterms[out] = fxs*fys*rn;
+	}
+      }
+    }
+  }
+
+ reuse:
+  {
+    // generate and save the dot products
+    float *sv = (float *)src->matrix.buffer;
+    float *tv = (float *)tgt->matrix.buffer;
+    for (wx = 0; wx <= Nx/2; wx++) {
+      for (wy = 0; wy <= Ny/2; wy++) {
+	float fsum = 0.0;
+	int Nf = (wx + wy*Nx)*Nx*Ny;
+	for (ix = 0; ix < Nx; ix++) {
+	  for (iy = 0; iy < Ny; iy++) {
+	    float Ft = Fterms[Nf + ix + iy*Nx];
+	    float Fi = sv[ix + iy*Nx];
+	    fsum += Ft*Fi;
+	    // fprintf (stderr, "%d,%d : %d,%d : %f %f %f\n", ix, iy, wx, wy, Ft, Fi, fsum);
+	  }
+	}
+	tv[wx + wy*Nx] = fsum;
+      }
+    }
+  }
+  return (TRUE);
+
+ usage:
+  gprint (GP_ERR, "USAGE: dft2d (input) to (output)\n");
+  gprint (GP_ERR, "  NOTE: Nx & Ny must both be even\n");
+  return (FALSE);
+}
+
+
+/*** 
+
+     2D direct fourier transform:
+
+     we have an image src of size [Nx,Ny].  in each direction,
+     we can generate the following fourier components:
+
+     cos(0*2pi*x/Nx), cos(1*2pi*x/Nx),... cos((Nx/2)  *2pi*x/Nx), 
+     sin((Nx/2-1)*2pi*x/Nx), sin((Nx/2-2)*2pi*x/Nx),.. sin(1*2pi*x/Nx)
+
+     Note that there are 2 fewer sin elements than there are cos
+     elements:
+       cos(0) = 1 -> DC term
+       sin(0) = 0
+       sin(pi*x) = 0 for integer values of x
+		     
+     we represent the output with the highest frequency in the middle and the DC element
+     in the 0,0 corner
+
+     to generate the output elements, we need to take the dot product 
+     of the input image with an image for the given x,y frequency
+ 
+     thus, for an Nx*Ny input image we need a cube of fourier terms of dimension
+     Nx*Ny*(Nx*Ny) (only one of which is trivial)
+
+     In the assumption that we will likely re-do the same size image multiple times, I
+     store the fourier factor array once it is generated.
+
+     NOTE: do I need to require Nx & Ny even?
+
+ ***/
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/dimenup.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/dimenup.c	(revision 37732)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/dimenup.c	(revision 37733)
@@ -23,8 +23,5 @@
     return (FALSE);
   }
-
-  gfits_free_matrix (&buf[0].matrix);
-  gfits_free_header (&buf[0].header);
-  CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0);
+  ResetBuffer (buf, Nx, Ny, -32, 0.0, 1.0);
 
   out = (float *) buf[0].matrix.buffer;
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c	(revision 37732)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/init.c	(revision 37733)
@@ -32,4 +32,5 @@
 int densify          PROTO((int, char **));
 int device           PROTO((int, char **));
+int dft2d            PROTO((int, char **));
 int dimendown        PROTO((int, char **));
 int dimenup          PROTO((int, char **));
@@ -196,4 +197,5 @@
   {1, "densify",      densify,          "create an image histogram from a set of vectors"},
   {1, "device",       device,           "set / get current graphics device"},
+  {1, "dft2d",        dft2d,            "2D discrete fourier transform"},
   {1, "dimendown",    dimendown,        "convert image to vector"},
   {1, "dimenup",      dimenup,          "convert vector to image"},
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/read_vectors.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 37732)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 37733)
@@ -5,4 +5,5 @@
 
 void read_vectors_cleanup ();
+int read_table_sizes (Header *header);
 
 int datafile (int argc, char **argv) {
@@ -331,15 +332,14 @@
 
   off_t Nbytes;
-  int i, j, N, Nextend, Ny, Binary, vecType, padIfShort;
-  char type[16], ID[80], *CCDKeyword;
+  int i, j, N, Ny, Binary, vecType;
+  char type[16], ID[80];
   FTable table;
   Header header;
   Vector **vec;
-  int FITS_TRANSPOSE;
 
   table.buffer = NULL;
   header.buffer = NULL;
 
-  FITS_TRANSPOSE = FALSE;
+  int FITS_TRANSPOSE = FALSE;
   if ((N = get_argument (argc, argv, "-transpose"))) {
     remove_argument (N, &argc, argv);
@@ -347,5 +347,5 @@
   }
 
-  CCDKeyword = NULL;
+  char *CCDKeyword = NULL;
   if ((N = get_argument (argc, argv, "-keyword"))) {
     remove_argument (N, &argc, argv);
@@ -354,5 +354,5 @@
   }
 
-  padIfShort = FALSE;
+  int padIfShort = FALSE;
   if ((N = get_argument (argc, argv, "-pad-if-short"))) {
     remove_argument (N, &argc, argv);
@@ -360,8 +360,24 @@
   }
 
-  Nextend = -1;
+  int getSizes = FALSE;
+  if ((N = get_argument (argc, argv, "-sizes"))) {
+    remove_argument (N, &argc, argv);
+    getSizes = TRUE;
+  }
+
+  int Nextend = -1;
   if ((N = get_argument (argc, argv, "-extnum"))) {
     remove_argument (N, &argc, argv);
     Nextend = atoi (extname);
+  }
+
+  int start = 0;
+  int Nrows = -1; // -1 : read entire table
+  if ((N = get_argument (argc, argv, "-range"))) {
+    remove_argument (N, &argc, argv);
+    start = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+    Nrows = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -406,6 +422,24 @@
     }
     if (!gfits_load_header (f, &header)) ESCAPE ("error reading header for extension %d\n", Nextend);
-    if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension %d\n", Nextend);
-
+
+    if (getSizes) {
+      read_table_sizes (&header);
+      if (CCDKeyword != NULL) free (CCDKeyword); 
+      gfits_free_header (&header); 
+      return TRUE;
+    }
+
+    if (Nrows == -1) {
+      Nrows = header.Naxis[1] - start;
+    }
+    if (start < 0) ESCAPE ("invalid range: start < 0\n");
+    if (start >= header.Naxis[1]) ESCAPE ("invalid range: start >= Ny (%d)\n", header.Naxis[1]);
+    if (Nrows < 0) ESCAPE ("invalid range: Nrows < 0\n");
+    if (start + Nrows > header.Naxis[1]) ESCAPE ("invalid range: start + Nrows > Ny (%d)\n", header.Naxis[1]);
+
+    // Ny = 100, start = 0, Nrows = -1 -> Nrows => 100
+    // Ny = 100, start = 10, Nrows = 90
+
+    if (!gfits_fread_ftable_range (f, padIfShort, &table, start, Nrows)) ESCAPE ("error reading table for extension %d\n", Nextend);
   } else {
     if (CCDKeyword == NULL) {
@@ -437,5 +471,17 @@
 	continue;
       }
-      if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension\n");
+
+      if (Nrows == -1) {
+	Nrows = header.Naxis[1] - start;
+      }
+      if (start < 0) ESCAPE ("invalid range: start < 0\n");
+      if (start >= header.Naxis[1]) ESCAPE ("invalid range: start >= Ny (%d)\n", header.Naxis[1]);
+      if (Nrows < 0) ESCAPE ("invalid range: Nrows < 0\n");
+      if (start + Nrows > header.Naxis[1]) ESCAPE ("invalid range: start + Nrows > Ny (%d)\n", header.Naxis[1]);
+
+      if (!gfits_fread_ftable_range (f, padIfShort, &table, start, Nrows)) ESCAPE ("error reading table for extension %d\n", Nextend);
+
+      // if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension\n");
+
       break;
     }
@@ -546,2 +592,19 @@
   if (vec) free (vec);
 }
+
+// read -fits foo -sizes -- Nx, Ny, Nfields -> $table:Nx, $table:Ny, $table:$Nfields
+// read -fits foo -fields 
+
+int read_table_sizes (Header *header) {
+  
+  int Nfields;
+
+  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);
+
+  set_int_variable ("table:Nx", header->Naxis[0]);
+  set_int_variable ("table:Nx", header->Naxis[0]);
+  set_int_variable ("table:Nfields", Nfields);
+
+  return TRUE;
+}
+
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/include/dvomath.h
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/include/dvomath.h	(revision 37732)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/include/dvomath.h	(revision 37733)
@@ -187,4 +187,5 @@
 int           ListBuffersToList     PROTO((char *name));
 int           CreateBuffer          PROTO((Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale));
+int           ResetBuffer           PROTO((Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale));
 Buffer       *SelectBuffer          PROTO((char *name, int mode, int verbose));
 void          dump_buffers          PROTO((int n));  /* deprecated? */
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/lib.shell/BufferOps.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/lib.shell/BufferOps.c	(revision 37732)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/lib.shell/BufferOps.c	(revision 37733)
@@ -88,4 +88,12 @@
 }
 
+int ResetBuffer (Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale) {
+
+  gfits_free_matrix (&buf[0].matrix);
+  gfits_free_header (&buf[0].header);
+  CreateBuffer (buf, Nx, Ny, bitpix, bzero, bscale);
+  return TRUE;
+}
+
 int CreateBuffer (Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale) {