Index: /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/tiltystreak.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/tiltystreak.c	(revision 24521)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/tiltystreak.c	(revision 24522)
@@ -1,2 +1,3 @@
+/* -*- c-file-style: "Ellemtel" -*- */
 /* tiltystreak.c - fit and subtract tilted streaks from a row of cells */
 /* Assumes that biastool has already been run -- zero streakiness near bias */
@@ -52,4 +53,5 @@
 static int parabola_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]);
 int linearfit(int npt, double *y, int nx, double *x, double *param);
+static void syntax(const char *prog);
 
 static int nmedian_buf=0;
@@ -59,5 +61,5 @@
 int main(int argc, const char* argv[])
 {
-   int i, row, col, ext;
+   int i, j, row, col, ext;
    HeaderUnit ihu = fh_create();
    HeaderUnit ehu = NULL;
@@ -66,10 +68,13 @@
    int nextend;
    CELL otarow[NCELL];
-   int naxis, naxis1, naxis2, naxis3;
+   int naxis, naxis1=0, naxis2=0;
+   int n1, n2;
    int prescan1, ovrscan1, ovrscan2;
    double bzero_d;
    int bzero = 32768, firstime = 1, algo=2;
+   int cellmask[NCELL*NCELL];
 
    if(argc < 2) {
+      syntax(argv[0]);
       fprintf(stderr, "error: syntax: tiltystreak Uncompressed_GPC1_MEF.fits\n");
       exit(EXIT_FAILURE);
@@ -77,17 +82,43 @@
    ifilename = argv[1];
 
-/* Choose different algorithms? */
-   if(argc > 2 && strncmp(argv[2], "alg=", 4) == 0) {
-      if(strcmp(argv[2]+4, "line") == 0) {
-	 algo = 1;
-      } else if(strcmp(argv[2]+4, "parabola") == 0) {
-	 algo = 2;
-      } else if(strcmp(argv[2]+4, "bent") == 0) {
-	 algo = 3;
+/* Defaults */
+   for(i=0; i<NCELL*NCELL; i++) cellmask[i] = 1;	/* Use all cells */
+   algo = 2;			/* Default algorithm is a parabola */
+
+/* Process arguments */
+   for(i=2; i<argc; i++) {
+
+/* Veto cells by a 64 digit mask (0/1 for cells 0-63) */
+      if(strncmp(argv[i], "mask=", 5) == 0) {	/* mask=64_digits */
+	 if(strlen(argv[i]) < 5+NCELL*NCELL) {
+	    fprintf(stderr, "\rerror: cannot get mask from `%s'\n", argv[i]);
+	    exit(EXIT_FAILURE);
+	 } else {
+	    for(j=0; j<NCELL*NCELL; j++) cellmask[j] = argv[i][5+j] != '0';
+	 }
+
+      } else if(strncmp(argv[2], "alg=", 4) == 0) {	/* Algorithm */
+	 if(strcmp(argv[2]+4, "line") == 0) {
+	    algo = 1;
+	 } else if(strcmp(argv[2]+4, "parabola") == 0) {
+	    algo = 2;
+	 } else if(strcmp(argv[2]+4, "bent") == 0) {
+	    algo = 3;
+	 } else {
+	    fprintf(stderr, "error: `%s' is not a valid algorithm\n", argv[2]);
+	    exit(EXIT_FAILURE);
+	 }
+
+      } else if(strncmp(argv[i], "help", 4) == 0) {	/* help output */
+	 syntax(argv[0]);
+	 exit(EXIT_SUCCESS);
+
       } else {
-	 fprintf(stderr, "error: `%s' is not a valid algorithm\n", argv[2]);
+	 fprintf(stderr, "\rerror: unrecognized option `%s'\n", argv[i]);
+	 syntax(argv[0]);
 	 exit(EXIT_FAILURE);
       }
    }
+
 
    if (fh_file(ihu, ifilename, FH_FILE_RDWR) != FH_SUCCESS) {
@@ -139,11 +170,10 @@
 
 /* Positioned at the desired cell now */
-	 naxis3 = 1;
 	 if (fh_get_PRESCAN1(ehu, &prescan1) != FH_SUCCESS ||
 	     fh_get_OVRSCAN1(ehu, &ovrscan1) != FH_SUCCESS ||
 	     fh_get_OVRSCAN2(ehu, &ovrscan2) != FH_SUCCESS ||
 	     fh_get_NAXIS(ehu, &naxis) != FH_SUCCESS ||
-	     fh_get_NAXIS1(ehu, &naxis1) != FH_SUCCESS ||
-	     fh_get_NAXIS2(ehu, &naxis2) != FH_SUCCESS) {
+	     fh_get_NAXIS1(ehu, &n1) != FH_SUCCESS ||
+	     fh_get_NAXIS2(ehu, &n2) != FH_SUCCESS) {
 	    fprintf(stderr, "error: Cannot get NAXIS* in extension %d row %d col %d\n", ext, row, col);
 	    exit(EXIT_FAILURE);
@@ -157,4 +187,6 @@
 /* Make some space now that we finally know the size */
 	 if(firstime) {
+	    naxis1 = n1;
+	    naxis2 = n2;
 	    firstime = 0;
 	    for(i=0; i<NCELL; i++) {
@@ -167,4 +199,11 @@
 	    Bbuf = (double *)calloc(naxis2, sizeof(double));
 	    Xbuf = (double *)calloc(NPOLY*naxis2, sizeof(double));
+	 }
+	 else {
+	    if(naxis1 != n1 || naxis2 != n2) {
+	       fprintf(stderr, "error: All cells must be same size (%dx%d != %dx%d)\n",
+		       n1, n2, naxis1, naxis2);
+	       exit(EXIT_FAILURE);
+	    }
 	 }
 
@@ -214,6 +253,7 @@
 #ifdef DEBUG
       for(col=0; col<NCELL; col++) {
-	 printf("%1d %1d %3d %3d %6d %6d %6d %6d %8.3f\n",
+	 printf("%1d %1d %3d %3d 0x%08x %6d %6d %6d %6d %8.3f\n",
 		col, row, otarow[col].cell, otarow[col].mefxt, 
+		(int)otarow[col].hu,
 		otarow[col].bias, otarow[col].sky, otarow[col].rms, 
 		otarow[col].nsat, otarow[col].ampl);
@@ -258,9 +298,19 @@
 /* Return to all the cells in the MEF and rewrite the data */
       for(col=0; col<NCELL; col++) {
+#ifdef DEBUG
+	 printf("Now writing cell %d %d ext %d (hu 0x%08x)\n", 
+		col, row, otarow[col].mefxt, (int)otarow[col].hu);
+#endif
 	 if(otarow[col].mefxt < 0) continue;
+/* Skip if cell mask not set */
+	 if(cellmask[col+NCELL*row] != 1) continue;
 /* This seeks back to the start of data from this cell again... */
 	 ehu = otarow[col].hu;
 	 fh_ehu(ehu, 0);
 
+#ifdef DEBUG
+	 printf("Now writing cell %d %d ext %d (hu 0x%08x)\n", 
+		col, row, otarow[col].mefxt, (int)otarow[col].hu);
+#endif
 	 if (fh_write_padded_image(ehu, fh_file_desc(ehu), otarow[col].data,
 				   naxis1*naxis2*sizeof(short),
@@ -272,5 +322,5 @@
       }
 
-   }
+   }	/* Loop over rows */
 
    fh_destroy(ihu);
@@ -439,5 +489,5 @@
       A = Abuf[j] - (par[0] + j*par[1] + j*j*par[2]);
 #ifdef DEBUG
-      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f\n", j, Abuf[j], A);
+//      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f\n", j, Abuf[j], A);
 #endif
       for(k=0; k<NCELL; k++) {
@@ -504,6 +554,6 @@
       B = Bbuf[j] - (Bpar[0] + j*Bpar[1] + j*j*Bpar[2]);
 #ifdef DEBUG
-      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f %9.3f %9.3f\n", 
-				      j, Abuf[j], A, Bbuf[j], B);
+//      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f %9.3f %9.3f\n", 
+//				      j, Abuf[j], A, Bbuf[j], B);
 #endif
 /* Parabola passes through 0 at i=sx+nx, A at sx+nx/6, B at sx+nx/2 */
@@ -590,6 +640,6 @@
 
 #ifdef DEBUG
-      if(otarow[3].mefxt == 4) printf("%4d %4d %9.3f %9.3f %9.3f\n", 
-				      j, i0, A, B, C);
+//      if(otarow[3].mefxt == 4) printf("%4d %4d %9.3f %9.3f %9.3f\n", 
+//				      j, i0, A, B, C);
 #endif
 
@@ -609,2 +659,10 @@
    return(0);
 }
+
+static void syntax(const char *prog)
+{
+   printf("Syntax: %s mef_file [options]\n", prog);
+   printf("   where options include:\n\n");
+   printf(" mask=0101...   64 digits to work on cells 0:63.\n");
+   printf(" alg=N          Use fit algorithm 1=linear, 2=parabola (default)\n");   printf("                   3=bent linear (not yet fully implemented)\n");
+}