Index: /tags/imclean-1-4-0/Ohana/src/imclean/Makefile
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/Makefile	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/Makefile	(revision 10354)
@@ -0,0 +1,29 @@
+default: imclean
+help:
+	@echo "make options: imclean (default)"
+
+include ../../Configure
+HOME    =       $(ROOT)/src/imclean
+BIN	=	$(HOME)/bin
+LIB	=	$(HOME)/lib
+SRC	=	$(HOME)/src
+MAN	=	$(HOME)/doc
+INC	= 	$(HOME)/include
+include ../../Makefile.Common
+
+LDFLAGS := 	-ldvo -lFITS -lohana $(LDFLAGS)
+
+imclean: $(BIN)/imclean.$(ARCH)
+install: $(DESTBIN)/imclean
+
+IMCLEAN = \
+$(SRC)/imclean.$(ARCH).o	$(SRC)/ConfigInit.$(ARCH).o \
+$(SRC)/args.$(ARCH).o		$(SRC)/sort_stars.$(ARCH).o \
+$(SRC)/AdjustHeader.$(ARCH).o	$(SRC)/find_group.$(ARCH).o \
+$(SRC)/find_trails.$(ARCH).o	$(SRC)/find_line.$(ARCH).o \
+$(SRC)/wstars.$(ARCH).o	        $(SRC)/LoadStarsDophot.$(ARCH).o \
+$(SRC)/LoadStarsSex.$(ARCH).o	$(SRC)/LoadStarsChad.$(ARCH).o	\
+$(SRC)/wfits.$(ARCH).o
+
+$(IMCLEAN): $(INC)/imclean.h
+$(BIN)/imclean.$(ARCH): $(IMCLEAN)
Index: /tags/imclean-1-4-0/Ohana/src/imclean/bin/.cvsignore
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/bin/.cvsignore	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/bin/.cvsignore	(revision 10354)
@@ -0,0 +1,2 @@
+*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix
+*.linrh
Index: /tags/imclean-1-4-0/Ohana/src/imclean/doc/ChangeLog.txt
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/doc/ChangeLog.txt	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/doc/ChangeLog.txt	(revision 10354)
@@ -0,0 +1,12 @@
+
+- imclean-1-3 : 2006.08.23
+  * converted to gfits APIs (forces libfits 1.6)
+
+2005.10.20: imclean-1-2
+	minor changes to fix -Wall errors
+	fixed exit status
+	synched with libohana (v1.5) / libfits (v1.4) changes
+
+2005.07.06: imclean-1-1
+	some changes to use the autocoded version of the SMPData structure (was Stars)
+	support for both fits and text output formats
Index: /tags/imclean-1-4-0/Ohana/src/imclean/doc/VERSIONS
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/doc/VERSIONS	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/doc/VERSIONS	(revision 10354)
@@ -0,0 +1,10 @@
+
+tag names used by imclean:
+
+AG         : Comment
+imclean-1-0 : first version under CVS
+
+USES
+libfits-1-0
+libohana-1-0
+
Index: /tags/imclean-1-4-0/Ohana/src/imclean/include/imclean.h
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/include/imclean.h	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/include/imclean.h	(revision 10354)
@@ -0,0 +1,50 @@
+# include <ohana.h>
+# include <dvo.h>
+
+enum {DOPHOT, CHAD, SEXTRACT};
+
+int    MODE;
+int    FITS_OUTPUT;
+int    VERBOSE;
+int    RESET;
+int    FORCE_RUN;
+int    PROVIDE_ASTROM;
+int    NEWPHOTCODE;
+
+/* global variables set in parameter file */
+char   *PHOTCODE;
+char   PhotCodeFile[256];
+char   AstromFile[256];
+
+double DEFAULT_ERROR;
+double RADIUS;
+double TRAIL_WIDTH;
+int    NBINS;
+int    NPTSINLINE;
+double MIN_DENSITY;
+double NSIGMA;
+
+double RA, DEC, ZERO_POINT, MIN_SN_FSTAT;
+int CHAR_LINE, TYPE_FIELD, AP_FIELD, PSF_FIELD, HEADER_COORDS;
+
+int FIX_KEYWORD;
+char **KEYWORD, **KEYVALU, **KEYFMT;
+
+SMPData *LoadStarsDophot (char *filename, int *nstars, Header *header);
+SMPData *LoadStarsChad (char *filename, int *nstars, Header *header);
+SMPData *LoadStarsSex (char *filename, int *nstars, Header *header);
+
+void ConfigInit (int *argc, char **argv);
+void AdjustHeader (Header *header);
+void sort_stars (SMPData *X, int N);
+void find_trails (SMPData *stars, int Nstars);
+int find_group (SMPData *stars, char *mark, int Npts, int i, double *ANGLE);
+void wstars (char *filename, SMPData *stars, int Nstars, Header *header);
+
+void fix_total (SMPData *stars, int Nstars, Header *header);
+int find_line (SMPData *stars, char *mark, int Npts, int i, double *M, double *B, double Angle);
+void find_better_line (SMPData *stars, char *mark, int Npts, int i, double *M, double *B, int axis);
+
+void help ();
+void args (int argc, char **argv);
+void wfits (char *filename, SMPData *stars, int Nstars, Header *header);
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/AdjustHeader.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/AdjustHeader.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/AdjustHeader.c	(revision 10354)
@@ -0,0 +1,110 @@
+# include "imclean.h"
+
+void AdjustHeader (Header *header) {
+
+  int i;
+  double value;
+  char line[256];
+
+  if (FIX_KEYWORD) {
+    for (i = 0; i < FIX_KEYWORD; i++) {
+      if (!strcmp (KEYFMT[i], "%f")) {
+	value = atof (KEYVALU[i]);
+	gfits_modify (header, KEYWORD[i], "%le", 1, value);
+      } else {
+	gfits_modify (header, KEYWORD[i], "%s", 1, KEYVALU[i]);
+      }
+    }
+  }
+
+  if (NEWPHOTCODE) {
+    /* we are going to write the photcode into the header
+       here we are just checking that the photcode provided
+       is a valid code */
+    if (!LoadPhotcodes (PhotCodeFile)) {
+      fprintf (stderr, "ERROR: can't load photcodes\n");
+      exit (1);
+    }
+    if (!GetPhotcodebyName (PHOTCODE)) {
+      fprintf (stderr, "ERROR: photcode not found in photcode table\n");
+      exit (1);
+    }
+    gfits_modify (header, "PHOTCODE", "%s", 1, PHOTCODE);
+  }    
+
+  if (!HEADER_COORDS) {
+    int rh, rm, dd, dm;
+    float rs, ds;
+    RA = RA / 15.0;
+    rh = RA; /* rh is int */
+    rm = 60.0 * (RA - rh);
+    rs = 3600 * (RA - rh - rm / 60.0);
+    sprintf (line, "%02d:%02d:%05.2f", rh, rm, rs);
+    gfits_modify (header, "RA", "%s", 1, line);
+    dd = DEC;
+    dm = 60.0 * (DEC - dd);
+    ds = 3600 * (DEC - dd - dm / 60.0);
+    sprintf (line, "%02d:%02d:%05.2f", dd, dm, ds);
+    gfits_modify (header, "DEC", "%s", 1, line);
+  }
+ 
+  if (PROVIDE_ASTROM) {
+
+    Header astrom_header;
+    Coords coords;
+
+    if (!gfits_read_header (AstromFile, &astrom_header)) {
+      fprintf (stderr, "ERROR: can't get astrometry from %s\n", AstromFile);
+      exit (1);
+    }
+    if (!GetCoords (&coords, &astrom_header)) {
+      fprintf (stderr, "ERROR: no astrometric solution in header\n");
+      exit (1);
+    }
+    /*** use PutCoords to update header ***/
+    if (coords.Npolyterms > 1) {
+      gfits_modify (header, "CTYPE1",   "%s",  1, "RA---PLY");
+      gfits_modify (header, "CTYPE2",   "%s",  1, "DEC--PLY");
+    } else {
+      gfits_modify (header, "CTYPE1",   "%s",  1, "RA---TAN");
+      gfits_modify (header, "CTYPE2",   "%s",  1, "DEC--TAN");
+    }    
+    gfits_modify (header, "NASTRO",   "%d", 1, 1); 
+
+    gfits_modify (header, "CDELT1",   "%le", 1, coords.cdelt1); 
+    gfits_modify (header, "CDELT2",   "%le", 1, coords.cdelt2);
+    gfits_modify (header, "CRVAL1",   "%lf", 1, coords.crval1);
+    gfits_modify (header, "CRVAL2",   "%lf", 1, coords.crval2);  
+    gfits_modify (header, "CRPIX1",   "%lf", 1, coords.crpix1);
+    gfits_modify (header, "CRPIX2",   "%lf", 1, coords.crpix2);
+    gfits_modify (header, "PC001001", "%le", 1, coords.pc1_1);
+    gfits_modify (header, "PC001002", "%le", 1, coords.pc1_2);
+    gfits_modify (header, "PC002001", "%le", 1, coords.pc2_1);
+    gfits_modify (header, "PC002002", "%le", 1, coords.pc2_2);
+    gfits_modify (header, "NPLYTERM", "%d", 1, coords.Npolyterms);
+    if (coords.Npolyterms > 1) {
+      /* RA Terms */
+      gfits_modify (header, "PCA1X2Y0", "%le", 1, coords.polyterms[0][0]);   /* polyterms[0]); */
+      gfits_modify (header, "PCA1X1Y1", "%le", 1, coords.polyterms[1][0]);   /* polyterms[1]); */
+      gfits_modify (header, "PCA1X0Y2", "%le", 1, coords.polyterms[2][0]);   /* polyterms[2]); */
+      
+      if (coords.Npolyterms > 2) {
+	gfits_modify (header, "PCA1X3Y0", "%le", 1, coords.polyterms[3][0]);   /* polyterms[3]); */
+	gfits_modify (header, "PCA1X2Y1", "%le", 1, coords.polyterms[4][0]);   /* polyterms[4]); */
+	gfits_modify (header, "PCA1X1Y2", "%le", 1, coords.polyterms[5][0]);   /* polyterms[5]); */
+	gfits_modify (header, "PCA1X0Y3", "%le", 1, coords.polyterms[6][0]);   /* polyterms[6]); */
+      }
+      /* Dec Terms */
+      gfits_modify (header, "PCA2X2Y0", "%le", 1, coords.polyterms[0][1]);   /* polyterms[7]); */
+      gfits_modify (header, "PCA2X1Y1", "%le", 1, coords.polyterms[1][1]);   /* polyterms[8]); */
+      gfits_modify (header, "PCA2X0Y2", "%le", 1, coords.polyterms[2][1]);   /* polyterms[9]); */
+      
+      if (coords.Npolyterms > 2) {
+	gfits_modify (header, "PCA2X3Y0", "%le", 1, coords.polyterms[3][1]);   /* polyterms[10]); */
+	gfits_modify (header, "PCA2X2Y1", "%le", 1, coords.polyterms[4][1]);   /* polyterms[11]); */
+	gfits_modify (header, "PCA2X1Y2", "%le", 1, coords.polyterms[5][1]);   /* polyterms[12]); */
+	gfits_modify (header, "PCA2X0Y3", "%le", 1, coords.polyterms[6][1]);   /* polyterms[13]); */
+      }
+    }
+  }
+}
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/ConfigInit.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/ConfigInit.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/ConfigInit.c	(revision 10354)
@@ -0,0 +1,36 @@
+# include "imclean.h"
+
+void ConfigInit (int *argc, char **argv) {
+  
+  char *config, *file;
+
+  /*** load configuration info ***/
+  file = SelectConfigFile (argc, argv, "ptolemy");
+  config = LoadConfigFile (file);
+  if (config == (char *) NULL) {
+    fprintf (stderr, "ERROR: can't find configuration file %s\n", file);
+    if (file != (char *) NULL) free (file);
+    exit (1);
+  }
+  if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
+
+  ScanConfig (config, "ZERO_PT",           "%lf", 0, &ZERO_POINT);
+  ScanConfig (config, "MIN_SN_FSTAT",      "%lf", 0, &MIN_SN_FSTAT);
+  ScanConfig (config, "DEFAULT_ERROR_FSTAT", "%lf", 0, &DEFAULT_ERROR);
+  ScanConfig (config, "DOPHOT_CHAR_LINE",  "%d", 0, &CHAR_LINE);
+  ScanConfig (config, "DOPHOT_TYPE_FIELD", "%d", 0, &TYPE_FIELD);
+  ScanConfig (config, "DOPHOT_PSF_FIELD",  "%d", 0, &PSF_FIELD);
+  ScanConfig (config, "DOPHOT_AP_FIELD",   "%d", 0, &AP_FIELD);
+  ScanConfig (config, "PHOTCODE_FILE",     "%s", 0, PhotCodeFile);
+
+  /* unique to markstar */
+  ScanConfig (config, "SEARCH_RADIUS",   "%lf", 0, &RADIUS);
+  ScanConfig (config, "TRAIL_WIDTH",     "%lf", 0, &TRAIL_WIDTH);
+  ScanConfig (config, "NANGLE_BINS",     "%d",  0, &NBINS);
+  ScanConfig (config, "NPTSINLINE",      "%d",  0, &NPTSINLINE);
+  ScanConfig (config, "MIN_DENSITY",     "%lf", 0, &MIN_DENSITY);
+  ScanConfig (config, "SPACE_SIGMA",     "%lf", 0, &NSIGMA); 
+
+  free (config);
+  free (file);
+}
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsChad.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsChad.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsChad.c	(revision 10354)
@@ -0,0 +1,133 @@
+# include "imclean.h"
+# define NBLOCK 100
+# define HIST_BINS 150
+
+SMPData *LoadStarsChad (char *filename, int *nstars, Header *header) {
+
+  FILE *f;
+  SMPData *stars;
+  int NSTARS, Nstars, i, N;
+  int status;
+  double x, y, m, sky, lsky;
+  char *buffer;
+  int Mhist[HIST_BINS], Shist[HIST_BINS], bin, sum;
+  double FWHMx, FWHMy, angle, flux;
+  int satfound, done;
+  double saturate, complete;
+  char line[256];
+
+  ALLOCATE (buffer, char, CHAR_LINE*NBLOCK);
+
+  f = fopen (filename, "r");
+  if (f == (FILE *) NULL) {
+    fprintf (stderr, "ERROR: can't find object file %s\n", filename);
+    exit (1);
+  }
+
+  /* zero things that will sum */
+  for (i = 0; i < HIST_BINS; i++) { Mhist[i] = Shist[i] = 0; }
+  
+  N = 0;
+  Nstars = 0;
+  NSTARS = 500;
+  ALLOCATE (stars, SMPData, NSTARS);
+
+  /* for now assume file 'header' is fixed-format */
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  if (strncasecmp (line, "#seeing", 7)) {
+    fprintf (stderr, "error in header, skipping\n");
+    exit (1);
+  }
+  sscanf (line, "%*s %lf", &FWHMx);
+  FWHMy = FWHMx;
+  angle = 0;
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+
+  /* read in data from obj file */
+  /* data is not fixed format for lines, read each line one-at-a-time */
+  for (i = 0; (status = fscanf (f, "%lf %lf %*f %*f %lf %lf %*f", &x, &y, &sky, &flux)) != EOF; i++) {
+
+    if (status != 4) {
+      fprintf (stderr, "format error in file %s, line %d\n", filename, i);
+      continue;
+    }
+
+    if (flux <= 0) continue;
+    m = -2.5*log10 (flux);
+
+    if (sky < 1.0) {
+      lsky = 0.0;
+    } else {
+      lsky = log10(sky);
+    }
+
+    bin = MAX (0.0, MIN (HIST_BINS, 10.0 * (m + 15.0)));  /* stick in 0.1 mag bins */
+    Mhist[bin] ++;
+
+    m = MIN (50.0, m);
+    m = MAX (-24.0, m);
+
+    stars[Nstars].fx = 0;
+    stars[Nstars].fy = 0;
+    stars[Nstars].df = 0;
+    stars[Nstars].Mgal = 50.0;;
+    stars[Nstars].Map = 50.0;
+    stars[Nstars].X = x;
+    stars[Nstars].Y = y;
+    stars[Nstars].M = m;
+    stars[Nstars].dM = 0.01;
+    stars[Nstars].dophot = 1;
+    stars[Nstars].sky = lsky;
+    Nstars++;
+    if (Nstars == NSTARS - 1) {
+      NSTARS += 500;
+      REALLOCATE (stars, SMPData, NSTARS);
+    }
+  }
+
+  /* look at histogram, find saturation and completion limits */
+  sum = 0;
+  for (i = 0; i < HIST_BINS; i++) {
+    sum += Mhist[i];
+    Shist[i] = sum;
+  }
+  satfound = done = FALSE;
+  saturate = complete = 0.0;
+  for (i = 0; (i < HIST_BINS) && !done; i++) {
+    if ((!satfound) && (Mhist[i] > 0)) {
+      saturate = 0.1*(i-1) - 15.0;
+      satfound = TRUE;
+    }
+    if (Shist[i] > 0.9*Shist[HIST_BINS - 1]) {
+      complete = 0.1*i - 15.0;
+      done = TRUE;
+    }
+  }
+  
+  gfits_modify (header, "ZERO_PT", "%lf", 1, ZERO_POINT);
+  gfits_modify (header, "FWHM_X", "%lf", 1, FWHMx);
+  gfits_modify (header, "FWHM_Y", "%lf", 1, FWHMy);
+  gfits_modify (header, "ANGLE", "%lf", 1, angle);
+  gfits_modify (header, "FSATUR", "%lf", 1, (saturate + ZERO_POINT));
+  gfits_modify (header, "FLIMIT", "%lf", 1, (complete + ZERO_POINT));
+  gfits_modify (header, "NSTARS", "%d", 1, N);
+  for (i = 1; i <= 9; i++) {
+    sprintf (line, "TDOPHOT%1d", i);
+    gfits_modify (header, line, "%d", 1, 0);
+  }
+
+  *nstars = Nstars;
+  return (stars);
+
+}
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsDophot.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsDophot.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsDophot.c	(revision 10354)
@@ -0,0 +1,176 @@
+# include "imclean.h"
+# define NBLOCK 100
+# define HIST_BINS 300 
+/* the mag histogram has range 0.0 to 30.0 */
+/* m = 0.1*Mhist[i] */
+
+SMPData *LoadStarsDophot (char *filename, int *nstars, Header *header) {
+
+  FILE *f;
+  SMPData *stars;
+  int NSTARS, Nstars, i, Nline, N;
+  int type, status;
+  double x, y, m, dm, sky, lsky, fx, fy, df, Mgal, Map;
+  char *buffer;
+  int Mhist[HIST_BINS], Shist[HIST_BINS], n[20], bin, sum;
+  double FWHMx, FWHMy, angle;
+  int satfound, done;
+  double saturate, complete;
+  char line[256];
+  int MedHist[2002], NMedHist;
+  double SMedHist, dMed;
+
+  N = NMedHist = 0;
+  bzero (MedHist, 2002*sizeof(int));
+
+  ALLOCATE (buffer, char, CHAR_LINE*NBLOCK);
+
+  f = fopen (filename, "r");
+  if (f == (FILE *) NULL) {
+    fprintf (stderr, "ERROR: can't find object file %s\n", filename);
+    exit (1);
+  }
+
+  /* zero things that will sum */
+  for (i = 0; i < HIST_BINS; i++) { Mhist[i] = Shist[i] = 0; }
+  for (i = 1; i <= 9; i++) { n[i] = 0; }
+
+  Nstars = 0;
+  NSTARS = 500;
+  ALLOCATE (stars, SMPData, NSTARS);
+
+  /* read average values from first line */
+  scan_line (f, line);
+  sscanf (line, "%*s %*s %*s %lf %lf %lf", &FWHMx, &FWHMy, &angle);
+
+  /* read in data from obj file */
+  while ((Nline = fread (buffer, CHAR_LINE, NBLOCK, f)) > 0) {
+    for (i = 0; i < Nline; i++) {
+      /* we are now using all entries on the *.obj line */
+      status = sscanf (&buffer[i*CHAR_LINE], "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", 
+		       &type, &x, &y, &m, &dm, &sky, &fx, &fy, &df, &Mgal, &Map);
+      if (status != 11) {
+	fprintf (stderr, "format error in file %s, line %d\n", filename, i);
+	continue;
+      }
+
+      /* dophot magnitudes can range from 99.999 to -99.999 
+	 realistic numbers are between -20 and 0
+	 outside, we should set the value to 50.0 to force saturation */
+      
+      n[type] ++;
+      if (type == 6) continue;
+      if (type == 9) continue;
+      if (type == 8) continue;
+      if (type == 16) continue;
+      
+      if (m > 0) continue;                /* skip stars which totally fail on fit */
+      if (dm == 0.0) dm = DEFAULT_ERROR;  /* stars with poor errors, get 25.5% errors */
+      dm = MIN (0.999, MAX (0.0, dm));    /* truncate dm to fit in range 0 - 999 on output */
+
+      /* need to accumulate the median histogram thingy */
+      dMed = Map - m;
+      if ((fabs(m) < 90) && (fabs(Map) < 90) && (fabs(dMed) < 1)) {
+	bin = 1000 * (dMed + 1);
+	MedHist[bin] ++;
+	NMedHist ++;
+      }
+      
+      /* dophot provides values which are -2.5*log(counts) */
+      m    = ((m    > -25) && (m    < 0)) ? m + ZERO_POINT    : 50.0;
+      Mgal = ((Mgal > -25) && (Mgal < 0)) ? Mgal + ZERO_POINT : 50.0;
+      Map  = ((Map  > -25) && (Map  < 0)) ? Map + ZERO_POINT  : 50.0;
+
+      if (sky < 1.0) {
+	lsky = 0.0;
+      } else {
+	lsky = log10(sky);
+      }
+
+      if (MIN_SN_FSTAT*dm > 1.0) continue;  /* skip stars with errors too large */
+          
+      switch (type) {
+      case 1:
+      case 4:
+      case 7:
+	bin = MAX (0, MIN (HIST_BINS - 1, 10.0 * m));  /* stick in 0.1 mag bins */
+	Mhist[bin] ++;
+      default:
+	if (df < 0.0) df += 360.0;
+	stars[Nstars].X      = x;
+	stars[Nstars].Y      = y;
+	stars[Nstars].M      = m;
+	stars[Nstars].Mgal   = Mgal;
+	stars[Nstars].Map    = Map;
+	stars[Nstars].dM     = dm;
+	stars[Nstars].dophot = type;
+	stars[Nstars].sky    = lsky;
+	stars[Nstars].fx     = fx;
+	stars[Nstars].fy     = fy;
+	stars[Nstars].df     = df;
+	Nstars++;
+	if (Nstars == NSTARS - 1) {
+	  NSTARS += 500;
+	  REALLOCATE (stars, SMPData, NSTARS);
+	}
+      }
+    }
+
+  }    
+
+  /* look at histogram, find saturation and completion limits */
+  sum = 0;
+  for (i = 0; i < HIST_BINS; i++) {
+    sum += Mhist[i];
+    Shist[i] = sum;
+  }
+  satfound = done = FALSE;
+  saturate = complete = 0.0;
+  for (i = 0; (i < HIST_BINS) && !done; i++) {
+    if ((!satfound) && (Mhist[i] > 0)) {
+      saturate = 0.1*(i-1);
+      satfound = TRUE;
+    }
+    if (Shist[i] > 0.9*Shist[HIST_BINS - 1]) {
+      complete = 0.1*i;
+      done = TRUE;
+    }
+  }
+  
+  SMedHist = 0;
+  for (i = 0; (i < 2002) && (SMedHist < NMedHist / 2); i++) {
+    SMedHist += MedHist[i];
+  }
+  if (i == 2002) {
+    fprintf (stderr, "error finding (Ap - Fit) median\n");
+    SMedHist = 0;
+  } else {
+    SMedHist = 0.001*i - 1;
+    fprintf (stderr, "(Ap - Fit) median = %f\n", SMedHist);
+  }
+
+  for (i = 0; i < Nstars; i++) {
+    stars[i].Mgal += SMedHist;
+    stars[i].M += SMedHist;
+  }    
+
+  gfits_modify (header, "ZERO_PT", "%lf", 1, ZERO_POINT);
+  gfits_modify (header, "FWHM_X", "%lf", 1, FWHMx);
+  gfits_modify (header, "FWHM_Y", "%lf", 1, FWHMy);
+  gfits_modify (header, "ANGLE", "%lf", 1, angle);
+  gfits_modify (header, "FSATUR", "%lf", 1, saturate);
+  gfits_modify (header, "FLIMIT", "%lf", 1, complete);
+  gfits_modify (header, "APMIFIT", "%lf", 1, SMedHist);
+  gfits_modify (header, "NSTARS", "%d", 1, N);
+  for (i = 1; i <= 9; i++) {
+    sprintf (line, "TDOPHOT%1d", i);
+    gfits_modify (header, line, "%d", 1, n[i]);
+  }
+
+  *nstars = Nstars;
+  return (stars);
+
+}
+
+/* this function should load the stars and immediately convert them to
+   have the ZERO_PT zero point */
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsSex.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsSex.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/LoadStarsSex.c	(revision 10354)
@@ -0,0 +1,142 @@
+# include "imclean.h"
+# define NBLOCK 100
+# define HIST_BINS 300 
+
+/* good for sextractor */
+
+SMPData *LoadStarsSex (char *filename, int *nstars, Header *header) {
+
+  FILE *f;
+  SMPData *stars;
+  int NSTARS, Nstars, i, Nline, N;
+  int type, status;
+  double x, y, m, dm, sky, lsky, ftype;
+  char *buffer;
+  int Mhist[HIST_BINS], Shist[HIST_BINS], n[10], bin, sum, flags;
+  double A, A2, S2, FWHMx, FWHMy, angle, Mgal, Map;
+  int gotFWHM, satfound, done;
+  double saturate, complete;
+
+  CHAR_LINE = 105;
+  TYPE_FIELD = 0;
+  
+  ALLOCATE (buffer, char, CHAR_LINE*NBLOCK);
+
+  f = fopen (filename, "r");
+  if (f == (FILE *) NULL) {
+    fprintf (stderr, "ERROR: can't find object file %s\n", filename);
+    exit (1);
+  }
+
+  /* zero things that will sum */
+  N = A = A2 = S2 = 0;
+  gotFWHM = FALSE;
+  for (i = 0; i < HIST_BINS; i++) { Mhist[i] = Shist[i] = 0; }
+  for (i = 1; i <= 9; i++) { n[i] = 0; }
+
+  Nstars = 0;
+  NSTARS = 500;
+  ALLOCATE (stars, SMPData, NSTARS);
+
+  /* read in data from obj file */
+  while ((Nline = fread (buffer, CHAR_LINE, NBLOCK, f)) > 0) {
+    for (i = 0; i < Nline; i++) {
+      status = sscanf (&buffer[i*CHAR_LINE + TYPE_FIELD], "%lf %lf %lf %lf %lf   %lf %lf %lf %lf %lf  %lf %d", 
+		       &ftype, &x, &y, &m, &dm, &sky, &FWHMx, &FWHMy, &angle, &Mgal, &Map, &flags);
+      if (status != 12) {
+	fprintf (stderr, "ERROR: format error in file %s, line %d\n", filename, i);
+	continue;
+      } 
+
+      if (flags > 7) continue;
+      /* if (m > 0) continue; skip stars which totally fail on fit */
+      if (dm == 0.0) dm = DEFAULT_ERROR;  
+
+      /* sextract can provide values which are -2.5*log(counts) */
+      m    += ZERO_POINT;
+      Mgal += ZERO_POINT;
+      Map  += ZERO_POINT;
+      m    = MIN (32.767, MAX (-32.767, m));
+      Mgal = MIN (32.767, MAX (-32.767, Mgal));
+      Map  = MIN (32.767, MAX (-32.767, Map));
+
+      if (sky < 1) {
+	lsky = 0.0;
+      } else {
+	lsky = log10(sky);
+      }
+
+      /* type = MAX (0, MIN (9, 5*log10(ftype) + 10)); */
+      switch (flags) {
+      case 4:
+      case 5:
+      case 6:
+      case 7:
+	type = 10;
+	break;
+      case 1:
+      case 2:
+      case 3:
+	type = 3;
+	break;
+      default:
+	type = 1;
+      }
+
+      if (MIN_SN_FSTAT*dm > 1.0) continue;
+      
+      bin = MAX (0, MIN (HIST_BINS - 1, 10.0 * m));  /* stick in 0.1 mag bins */
+      Mhist[bin] ++;
+      
+      stars[Nstars].X = x;
+      stars[Nstars].Y = y;
+      stars[Nstars].M = m;
+      stars[Nstars].dM = dm;
+      stars[Nstars].dophot = type;
+      stars[Nstars].sky = lsky;
+      
+      stars[Nstars].fx = FWHMx;
+      stars[Nstars].fy = FWHMx * (FWHMy/FWHMx);
+      stars[Nstars].df = angle;
+      stars[Nstars].Mgal = Mgal;
+      stars[Nstars].Map = Map;
+
+      Nstars++;
+      if (Nstars == NSTARS - 1) {
+	NSTARS += 500;
+	REALLOCATE (stars, SMPData, NSTARS);
+      }
+    }
+
+  }    
+
+  sum = 0;
+  for (i = 0; i < HIST_BINS; i++) {
+    sum += Mhist[i];
+    Shist[i] = sum;
+  }
+  satfound = done = FALSE;
+  saturate = complete = 0.0;
+  for (i = 0; (i < HIST_BINS) && !done; i++) {
+    if ((!satfound) && (Mhist[i] > 0)) {
+      saturate = 0.1*(i-1);
+      satfound = TRUE;
+    }
+    if (Shist[i] > 0.9*Shist[HIST_BINS - 1]) {
+      complete = 0.1*i;
+      done = TRUE;
+    }
+  }
+  
+  gfits_modify (header, "ZERO_PT", "%lf", 1, ZERO_POINT);
+  gfits_modify (header, "FWHM_X", "%lf", 1, FWHMx);
+  gfits_modify (header, "FWHM_Y", "%lf", 1, FWHMy);
+  gfits_modify (header, "ANGLE", "%lf", 1, angle);
+  gfits_modify (header, "FSATUR", "%lf", 1, saturate);
+  gfits_modify (header, "FLIMIT", "%lf", 1, complete);
+  gfits_modify (header, "NSTARS", "%d", 1, N);
+
+  *nstars = Nstars;
+  return (stars);
+
+}
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/args.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/args.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/args.c	(revision 10354)
@@ -0,0 +1,87 @@
+# include "imclean.h"
+# define NARGS 2  /* minimum is: addstar (filename) */
+
+void help () {
+
+  fprintf (stderr, "USAGE: imclean (file.fits) (file.obj) (file.cmp)\n");
+  exit (2);
+
+}
+
+void args (int argc, char **argv) {
+  
+  int i, N;
+
+  if (get_argument (argc, argv, "-help") ||
+      get_argument (argc, argv, "-h")) {
+    help ();
+  }
+
+  VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  FITS_OUTPUT = FALSE;
+  if ((N = get_argument (argc, argv, "-fits"))) {
+    FITS_OUTPUT = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  HEADER_COORDS = TRUE;
+  if ((N = get_argument (argc, argv, "-coords"))) {
+    remove_argument (N, &argc, argv);
+    RA = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    DEC = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    HEADER_COORDS = FALSE;
+  }
+
+  NEWPHOTCODE = FALSE;
+  if ((N = get_argument (argc, argv, "-p"))) {
+    NEWPHOTCODE = TRUE;
+    remove_argument (N, &argc, argv);
+    PHOTCODE = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  PROVIDE_ASTROM = FALSE;
+  if ((N = get_argument (argc, argv, "-astrom"))) {
+    remove_argument (N, &argc, argv);
+    strcpy (AstromFile, argv[N]);
+    remove_argument (N, &argc, argv);
+    PROVIDE_ASTROM = TRUE;
+  }
+
+  MODE = DOPHOT;
+  if ((N = get_argument (argc, argv, "-chad"))) {
+    MODE = CHAD;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-sex"))) {
+    MODE = SEXTRACT;
+    remove_argument (N, &argc, argv);
+  }
+
+  ALLOCATE (KEYWORD, char *, 64);
+  ALLOCATE (KEYVALU, char *, 64);
+  ALLOCATE (KEYFMT, char *, 64);
+  FIX_KEYWORD = 0;
+  while ((N = get_argument (argc, argv, "-key"))) {
+    i = FIX_KEYWORD;
+    remove_argument (N, &argc, argv);
+    KEYWORD[i] = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+    KEYFMT[i] = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+    KEYVALU[i] = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+    FIX_KEYWORD ++;
+    if (FIX_KEYWORD == 64) break;
+  }
+
+
+  if (argc != 4) help ();
+
+}
+
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/find_group.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/find_group.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/find_group.c	(revision 10354)
@@ -0,0 +1,61 @@
+# include "imclean.h"
+
+int find_group (SMPData *stars, char *mark, int Npts, int i, double *ANGLE) {
+ 
+  int j, N, bin, Ndegbin;
+  unsigned short int *A;
+  double Xo, Yo, dX, dY, angle, limit;
+
+   /* assign some parameter values */
+  if (mark[i]) return (FALSE);
+
+  Ndegbin = NBINS / 180.0;
+  ALLOCATE (A, unsigned short int, NBINS + 1);
+  bzero (A, (NBINS+1)*sizeof(short));
+
+  /* look for points concentrated in an angle bin */
+  Xo = stars[i].X;
+  Yo = stars[i].Y;
+  N = 0;
+  /* points east */
+  for (j = i + 1; (j < Npts) && ((dX = stars[j].X - Xo) < RADIUS); j++) {
+    dY = stars[j].Y - Yo;
+    if (fabs(dY) < RADIUS) {
+      angle = atan2 (dY,dX);
+      if (!finite(angle)) continue;  /* only NaN if dD = dR = 0 */
+      if (angle < 0) angle += M_PI;
+      bin = angle*DEG_RAD*Ndegbin;
+      A[bin] ++;
+      N ++;
+    }
+  }
+  /* points west */
+  for (j = i - 1; (j >= 0) && ((dX = Xo - stars[j].X) < RADIUS); j--) {
+    if (mark[j]) continue;
+    dY = stars[j].Y - Yo;
+    if (fabs(dY) < RADIUS) {
+      angle = atan2 (dY,dX);
+      if (!finite(angle)) continue;
+      if (angle < 0) angle += M_PI;
+      bin = angle*DEG_RAD*Ndegbin;
+      A[bin] ++;	
+      N ++;
+    }
+  }
+  if (N < 5) {
+    free (A);
+    return (FALSE);
+  }
+  limit = MAX (5, NSIGMA*sqrt((double)(N)/(double)(NBINS)));
+  for (j = 0; j < NBINS; j++) {
+    if (A[j] > limit) {
+      *ANGLE = j / (DEG_RAD*Ndegbin);
+      fprintf (stderr, "group: %f (%f %f)\n", *ANGLE, Xo, Yo);
+      free (A);
+      return (TRUE);
+    }
+  }
+  free (A);
+  return (FALSE);
+}
+  
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/find_line.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/find_line.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/find_line.c	(revision 10354)
@@ -0,0 +1,164 @@
+# include "imclean.h"
+
+int find_line (SMPData *stars, char *mark, int Npts, int i, double *M, double *B, double Angle) {
+  
+  int j, N;
+  double X, Y, X2, Y2, XY, m, b, det;
+  double dX, dY, Xo, Yo, angle;
+  char Flipped;
+  
+  Xo = stars[i].X;
+  Yo = stars[i].Y;
+  X = Xo;
+  Y = Yo;
+  X2 = Xo*Xo;
+  Y2 = Yo*Yo;
+  XY = Xo*Yo;
+  N = 1;
+  /* points to the right */
+  for (j = i + 1; (j < Npts) && ((dX = stars[j].X - Xo) < RADIUS); j++) {
+    if (mark[j]) continue;
+    dY = stars[j].Y - Yo;
+    if (fabs(dY) < RADIUS) {
+      angle = atan2 (dY,dX);
+      if (!finite(angle)) continue;
+      if (angle < 0) angle += M_PI;
+      if (fabs(angle - Angle) < 2*RAD_DEG) {
+	X += stars[j].X;
+	Y += stars[j].Y;
+	X2 += stars[j].X*stars[j].X;
+	Y2 += stars[j].Y*stars[j].Y;
+	N ++;
+	XY += stars[j].X*stars[j].Y;
+      }
+    }
+  }
+  /* points to the left */
+  for (j = i - 1; (j >= 0) && ((dX = Xo - stars[j].X) < RADIUS); j--) {
+    if (mark[j]) continue;
+    dY = stars[j].Y - Yo;
+    if (fabs(dY) < RADIUS) {
+      angle = atan2 (dY,dX);
+      if (!finite(angle)) continue;
+      if (angle < 0) angle += M_PI;
+      if (fabs(angle - Angle) < 2*RAD_DEG) {
+	X += stars[j].X;
+	Y += stars[j].Y;
+	X2 += stars[j].X*stars[j].X;
+	Y2 += stars[j].Y*stars[j].Y;
+	N ++;
+	XY += stars[j].X*stars[j].Y;
+      }
+    }
+  }
+  /* determine coeffs */
+  Flipped = 0;
+  det = 1.0 / (X2*N - X*X);
+  m = det * (XY*N - X*Y);
+  b = det * (X2*Y - XY*X);
+  if (fabs(m) > 1.1) { /* use a line of R = m*D + b instead */
+    /* fprintf (stderr, "high slope object: %f %f  -> ", m, b); */
+    det = 1.0 / (Y2*N - Y*Y);
+    m = det * (XY*N - X*Y);
+    b = det * (Y2*X - XY*Y);
+    Flipped = 1;
+    /* fprintf (stderr, "%f %f\n", m, b); */
+  }
+
+  *M = m;
+  *B = b;
+  return (Flipped);
+
+}
+
+
+void find_better_line (SMPData *stars, char *mark, int Npts, int i, double *M, double *B, int axis) {
+  
+  int j, N;
+  double X, Y, X2, Y2, XY, m, b, det, dist;
+  double delta, *path;
+  int NLINE, *line;
+
+  NLINE = 50;
+  ALLOCATE (line, int, NLINE);
+  ALLOCATE (path, double, NLINE);
+  
+  /* fit a line to points near line */
+  X = Y = X2 = Y2 = XY = N = 0;
+  m = *M;  b = *B;
+  for (j = 0; (j < Npts); j++) {
+    if (axis == 1) 
+      delta = stars[j].X - m*stars[j].Y - b;
+    else
+      delta = stars[j].Y - m*stars[j].X - b;
+    if (fabs(delta) < 2*TRAIL_WIDTH) {
+      X += stars[j].X;
+      Y += stars[j].Y;
+      X2 += stars[j].X*stars[j].X;
+      Y2 += stars[j].Y*stars[j].Y;
+      XY += stars[j].X*stars[j].Y;
+      line[N] = j;
+      path[N] = hypot (stars[j].X - stars[line[0]].X, stars[j].Y - stars[line[0]].Y);
+      N ++;
+      if (N == NLINE - 1) {
+	NLINE += 50;
+	REALLOCATE (line, int, NLINE);
+	REALLOCATE (path, double, NLINE);
+      }
+    }
+  }
+
+  if (N < NPTSINLINE) {
+    free (line);
+    return;
+  }
+
+  for (i = 0; i < N - NPTSINLINE + 1; i++) {
+    j = i + NPTSINLINE - 1;
+    dist = fabs (path[j] - path[i]);
+    if ((j - i) / dist < MIN_DENSITY) continue;
+    for (; (j < N) && (((j - i) / dist) > MIN_DENSITY); j++) {
+      dist = fabs (path[j] - path[i]);
+    }
+    if ((j == N) && (((j - i) / dist) > MIN_DENSITY)) j++;
+    j--;
+    for (; i < j; i++) {
+      mark[line[i]] = TRUE;
+      stars[line[i]].dophot = 0;
+    }
+    i--;
+  }
+  free (line);
+  free (path);
+
+  /* determine coeffs */
+  if (axis == 0) {
+    det = 1.0 / (X2*N - X*X);
+    m = det * (XY*N - X*Y);
+    b = det * (X2*Y - XY*X);
+  } else {
+    det = 1.0 / (Y2*N - Y*Y);
+    m = det * (XY*N - X*Y);
+    b = det * (Y2*X - XY*Y);
+  }
+
+  fprintf (stderr, "%f %f %d\n", m, b, N);
+
+  *M = m;
+  *B = b;
+
+}
+
+void fix_total (SMPData *stars, int Nstars, Header *header) {
+  
+  int Ngood, i;
+
+  Ngood = 0;
+  for (i = 0; i < Nstars; i++) {
+    if (stars[i].dophot != 0) Ngood ++;
+  }
+
+  gfits_modify (header, "NSTARS", "%d", 1, Ngood);
+  
+}
+  
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/find_trails.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/find_trails.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/find_trails.c	(revision 10354)
@@ -0,0 +1,23 @@
+# include "imclean.h"
+
+void find_trails (SMPData *stars, int Nstars) {
+
+  int i;
+  char *mark;
+  double Angle, axis, m, b;
+  
+  ALLOCATE (mark, char, Nstars);
+  bzero (mark, Nstars);
+  
+  for (i = 0; i < Nstars; i++) {
+    /* already marked, ignore */
+    if (mark[i]) continue;
+    if (find_group (stars, mark, Nstars, i, &Angle)) {
+      /* this point has an excess nearby concentration, find the line */
+      axis = find_line (stars, mark, Nstars, i, &m, &b, Angle);
+      find_better_line (stars, mark, Nstars, i, &m, &b, axis);
+      /* mark_trail (stars, mark, Nstars, i, m, b, axis); */
+    }
+  }
+}
+
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/imclean.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/imclean.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/imclean.c	(revision 10354)
@@ -0,0 +1,75 @@
+# include "imclean.h"
+
+int main (int argc, char **argv) {
+
+  Header header;
+  SMPData *stars;
+  int Nstars;
+
+  ConfigInit (&argc, argv);
+
+  args (argc, argv);
+
+  /* load in FITS header from image */
+  if (!gfits_read_header (argv[1], &header)) {
+    fprintf (stderr, "ERROR: can't find image file %s\n", argv[1]);
+    exit (1);
+  }
+
+  AdjustHeader (&header);
+  
+  switch (MODE) {
+  case DOPHOT:
+    stars = LoadStarsDophot (argv[2], &Nstars, &header);
+    break;
+  case CHAD:
+    stars = LoadStarsChad (argv[2], &Nstars, &header);
+    break;
+  case SEXTRACT:
+    stars = LoadStarsSex (argv[2], &Nstars, &header);
+    break;
+  default: 
+    fprintf (stderr, "unknown mode: %d\n", MODE);
+    exit (1);
+  }
+
+  sort_stars (stars, Nstars);
+  find_trails (stars, Nstars);  
+  fix_total (stars, Nstars, &header);
+
+  if (FITS_OUTPUT) {
+    wfits (argv[3], stars, Nstars, &header); 
+  } else {
+    wstars (argv[3], stars, Nstars, &header); 
+  }
+
+  fprintf (stderr, "SUCCESS\n");
+  exit (0);
+
+}
+
+/* based on fstat and markstar:
+
+   0) load config data, global parameters 
+   1) load header
+   2) load data from *.obj file
+   3) eliminate bad star types: 6, 8
+   4) identify trails
+   5) get statistics on remaining stars
+
+*/
+
+/* 
+imclean (file.fits) (file.obj) (file.cmp) 
+
+ [-p photcode]
+ [-chad]
+ [-sex]
+ [-coords RA DEC] 
+ [-astrom file]
+ [-v]
+ [-key name %f value]
+ [-key name %s value]
+ (maximum of 64 keywords can be changed)
+
+*/
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/sort_stars.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/sort_stars.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/sort_stars.c	(revision 10354)
@@ -0,0 +1,36 @@
+# include "imclean.h"
+
+void sort_stars (SMPData *X, int N) {
+
+  int l,j,ir,i;
+  SMPData tX;
+  
+  if (N < 1) return;
+  l = N >> 1;
+  ir = N - 1;
+  for (;;) {
+    if (l > 0) {
+      l--;
+      tX = X[l];
+    }
+    else {
+      tX = X[ir];
+      X[ir] = X[0];
+      if (--ir == 0) {
+	X[0] = tX;
+	return;
+      }
+    }
+    i = l;
+    j = (l << 1) + 1;
+    while (j <= ir) {
+      if (j < ir && X[j].X < X[j+1].X) j++;
+      if (tX.X < X[j].X) {
+	X[i] = X[j];
+	j += (i=j) + 1;
+      }
+      else j = ir + 1;
+    }
+    X[i] = tX;
+  }
+}
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/star_stats.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/star_stats.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/star_stats.c	(revision 10354)
@@ -0,0 +1,76 @@
+# include "imclean.h"
+
+star_stats (Header *header, SMPData *stars, int Nstars) {
+
+  int i;
+  
+  /* zero things that will sum */
+  N = A = A2 = S2 = 0;
+  gotFWHM = FALSE;
+
+  for (i = 0; i < HIST_BINS; i++) { Mhist[i] = Shist[i] = 0; }
+  for (i = 1; i <= 9; i++) { n[i] = 0; }
+
+  for (i = 0; i < Nstars; i++) {
+    n[type] ++;
+    switch (type) {
+    case 6:  /* just throw these ones out */
+    case 8:
+    case 9:
+      break;
+    case 1:
+      sscanf (&buffer[i*CHAR_LINE + AP_FIELD], "%lf", &ap);
+      if (ap < 99) {
+	apmifit = ap - m;
+	A += apmifit / SQ(df);
+	A2 += SQ(apmifit) / SQ(df);
+	S2 += 1.0 / SQ(df);
+      } 
+      if (!gotFWHM) {
+	sscanf (&buffer[i*CHAR_LINE + PSF_FIELD], "%lf %lf %lf ", &FWHMx, &FWHMy, &angle);
+	gotFWHM = TRUE;
+      }
+    case 4:
+    case 7:
+      bin = MAX (0.0, MIN (HIST_BINS, 10.0 * (m + 15.0)));  /* stick in 0.1 mag bins */
+      Mhist[bin] ++;
+    case 2:
+    case 3:
+    case 5:
+      fprintf (g, "%6.1f %6.1f %6.3f %03d %1d %3.1f\n", x, y, m+ZERO_POINT, (int)(1000*df), type, lsky);
+      N ++; 
+    }
+  }
+
+  Ap = A / S2;
+  Ap2 = sqrt(A2 / (S2) - Ap*Ap);
+  sum = 0.0;
+  for (i = 0; i < HIST_BINS; i++) {
+    sum += Mhist[i];
+    Shist[i] = sum;
+  }
+  satfound = done = FALSE;
+  for (i = 0; (i < HIST_BINS) && !done; i++) {
+    if ((!satfound) && (Mhist[i] > 0)) {
+      saturate = 0.1*(i-1) - 15.0;
+      satfound = TRUE;
+    }
+    if (Shist[i] > 0.9*Shist[HIST_BINS - 1]) {
+      complete = 0.1*i - 15.0;
+      done = TRUE;
+    }
+  }
+  
+  gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT);
+  gfits_modify (&header, "FWHM_X", "%lf", 1, FWHMx);
+  gfits_modify (&header, "FWHM_Y", "%lf", 1, FWHMy);
+  gfits_modify (&header, "ANGLE", "%lf", 1, angle);
+  gfits_modify (&header, "APMIFIT", "%lf", 1, Ap);
+  gfits_modify (&header, "dAPMIFIT", "%lf", 1, Ap2);
+  gfits_modify (&header, "FSATUR", "%lf", 1, (saturate + ZERO_POINT));
+  gfits_modify (&header, "FLIMIT", "%lf", 1, (complete + ZERO_POINT));
+  gfits_modify (&header, "NSTARS", "%d", 1, N);
+  for (i = 1; i <= 9; i++) {
+    sprintf (line, "TDOPHOT%1d\0", i);
+    gfits_modify (&header, line, "%d", 1, n[i]);
+  }
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/wfits.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/wfits.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/wfits.c	(revision 10354)
@@ -0,0 +1,27 @@
+# include "imclean.h"
+
+void wfits (char *filename, SMPData *stars, int Nstars, Header *header) {
+
+  Matrix matrix;
+  Header theader;
+  FTable table;
+
+  header[0].extend = TRUE;
+  header[0].Naxes = 0;
+  gfits_modify (header, "NAXIS",   "%d", 1, 0);
+  gfits_modify (header, "EXTEND",  "%t", 1, TRUE);
+  gfits_modify (header, "NEXTEND", "%d", 1, 1);
+
+  /* add in some keywords to specify the datatype & software version? */
+
+  /* create (empty) data matrix */
+  gfits_create_matrix (header, &matrix);
+    
+  table.header = &theader;
+  gfits_table_set_SMPData (&table, stars, Nstars);
+
+  gfits_write_header  (filename, header);
+  gfits_write_matrix  (filename, &matrix);
+  gfits_write_Theader (filename, &theader);
+  gfits_write_table   (filename, &table);
+}
Index: /tags/imclean-1-4-0/Ohana/src/imclean/src/wstars.c
===================================================================
--- /tags/imclean-1-4-0/Ohana/src/imclean/src/wstars.c	(revision 10354)
+++ /tags/imclean-1-4-0/Ohana/src/imclean/src/wstars.c	(revision 10354)
@@ -0,0 +1,43 @@
+# include "imclean.h"
+# define NCHAR 66 /* 65 char EXCLUDING return */
+
+void wstars (char *filename, SMPData *stars, int Nstars, Header *header) {
+  
+  int i, Nchar;
+  FILE *g;
+  char line[NCHAR + 3];
+
+  g = fopen (filename, "w");
+  if (g == (FILE *) NULL) {
+    fprintf (stderr, "ERROR: can't create output file %s\n", filename);
+    exit (1);
+  }
+
+  fwrite (header[0].buffer, 1, header[0].size, g);
+
+  for (i = 0; i < Nstars; i++) {
+    if (stars[i].dophot == 0) continue;
+    Nchar = snprintf (line, NCHAR, "%6.1f %6.1f %6.3f %03d %2d %3.1f %6.3f %6.3f %6.2f %6.2f %5.1f", 
+		      stars[i].X, stars[i].Y, stars[i].M, 
+		      (int)(1000*stars[i].dM), stars[i].dophot, stars[i].sky, 
+		      stars[i].Mgal, stars[i].Map, stars[i].fx, stars[i].fy, stars[i].df);
+    
+    /* this is just a little funny.  NCHAR (in) includes the trailing NULL, Nchar (out) excludes it */
+    if (Nchar != NCHAR - 1) {
+      fprintf (stderr, "funny line %d (%d)\n%s\n", i, Nchar, line);
+    } else {
+      fprintf (g, "%s\n", line);
+    }
+  }
+
+  fclose (g);
+
+}
+
+/*
+
+  63.6 2869.5 17.568 157 17.568 17.568 25.01 25.00 360.0 7 2.9
+
+  63.6 2869.5 17.568 157 7 2.9
+
+*/