Index: /trunk/Ohana/src/addstar/src/FilterStars.c
===================================================================
--- /trunk/Ohana/src/addstar/src/FilterStars.c	(revision 10897)
+++ /trunk/Ohana/src/addstar/src/FilterStars.c	(revision 10897)
@@ -0,0 +1,56 @@
+# include "addstar.h"
+
+Stars *FilterStars (Stars *instars, Image *image) {
+
+  int j, N;
+
+  /* modify resulting star list */
+  ALLOCATE (stars, Stars, image[0].nstars);
+  for (N = j = 0; j < image[0].nstars; j++) {
+    /* allow for some dynamic filtering of star list */
+    if (SNLIMIT && instars[j].dM > SNLIMIT) continue;
+    if (XMAX && (instars[j].X > XMAX)) continue;
+    if (XMIN && (instars[j].X < XMIN)) continue;
+    if (YMAX && (instars[j].Y > YMAX)) continue;
+    if (YMIN && (instars[j].Y < YMIN)) continue;
+    stars[N] = instars[j];
+
+    XY_to_RD (&stars[N].R, &stars[N].D, stars[N].X, stars[N].Y, &image[0].coords);
+    while (stars[N].R <    0.0) stars[N].R += 360.0;
+    while (stars[N].R >= 360.0) stars[N].R -= 360.0;
+    stars[N].found = -1;
+    stars[N].code = photcode;
+
+    if (SUBPIX) {
+      dMs = get_subpix (stars[N].X, stars[N].Y);
+      stars[N].M    -= dMs;
+      stars[N].Mgal -= dMs;
+      stars[N].Map  -= dMs;
+      dMs = scat_subpix (stars[N].X, stars[N].Y);
+      stars[N].dM = hypot (stars[N].dM, dMs);
+    }
+    N ++;
+  }
+  image[0].nstar = N;
+  REALLOCATE (stars, Stars, image[0].nstar);
+  free (instars);
+
+  if (VERBOSE) fprintf (stderr, "read %d stars from target file\n", image[0].nstar);
+  return (stars);
+}
+
+Stars *MergeStars (Stars *instars, int Ninstars, Stars *stars, int *Nstars) {
+
+  if (stars == NULL) {
+    ALLOCATE (stars, Stars, Ninstars);
+  } else {
+    REALLOCATE (stars, Stars, *Nstars + Ninstars);
+  }
+
+  for (j = 0, i = *Nstars; i < *Nstars + Ninstars; i++) {
+    stars[i] = instars[j];
+  }
+  
+  *Nstars += Ninstars;
+  return (stars);
+}
Index: /trunk/Ohana/src/addstar/src/LoadStars.c
===================================================================
--- /trunk/Ohana/src/addstar/src/LoadStars.c	(revision 10896)
+++ /trunk/Ohana/src/addstar/src/LoadStars.c	(revision 10897)
@@ -1,4 +1,4 @@
 
-/* read an image header corresponding to a CMF data block */
+/* read an image header corresponding to a CMF / CMP data block */
 Image *ReadImageHeader (Header *header) {
 
Index: /trunk/Ohana/src/addstar/src/ReadStarsFITS.c
===================================================================
--- /trunk/Ohana/src/addstar/src/ReadStarsFITS.c	(revision 10897)
+++ /trunk/Ohana/src/addstar/src/ReadStarsFITS.c	(revision 10897)
@@ -0,0 +1,107 @@
+# include "addstar.h"
+
+// given a file with the pointer at the start of the table block and the 
+// corresponding image header, load the stars from the table
+Stars *ReadStarsFITS (FILE *f, Header *header, Header *in_theader, unsigned int *nstars) {
+
+  int i, Nstars;
+  Header theader;
+  FTable table;
+  Stars *stars;
+  
+  if (in_theader == NULL) {
+    table.header = &theader;
+  } else {
+    table.header = in_theader;
+  }
+
+  /* load the table header (skip this?) */
+  if (in_theader == NULL) {
+    if (!gfits_fread_header (f, &theader)) {
+      fprintf (stderr, "ERROR: can't read table header\n");
+      exit (1);
+    }
+  } else {
+    Nskip = in_theader[0].size;
+    fseek (f, Nskip, SEEK_CUR); 
+  }
+
+  /* load the table data */
+  if (!gfits_fread_ftable_data (f, &table)) {
+    fprintf (stderr, "ERROR: can't read table header\n");
+    exit (1);
+  }
+
+  if (!gfits_scan (header, "EXTTYPE", "%s", 1, type)) {
+    strcpy (type, "SMPDATA");
+  }
+
+  stars = NULL;
+  if (!strcmp (type, "SMPDATA")) {
+    stars = ConvertSMPDATA (&table, &Nstars);
+  }
+  if (!strcmp (type, "PS1DATA")) {
+    stars = ConvertPS1DATA (&table, &Nstars);
+  }
+  if (stars == NULL) {
+    fprintf (stderr, "ERROR: invalid table type %s\n", type);
+    exit (1);
+  }
+  if (*nstars != Nstars) {
+    fprintf (stderr, "ERROR: inconsistent number of stars? %d vs %d\n", *nstars, Nstars);
+    exit (1);
+  }
+  return stars;
+}
+
+Stars *ConvertSMPDATA (FTable *table, int *nstars) {
+
+  int i, Nstars;
+  Stars *stars;
+  SMPData *smpdata;
+
+  smpdata = gfits_table_get_SMPData (table, &Nstars, NULL);
+
+  ALLOCATE (stars, Stars, Nstars);
+  for (i = 0; i < Nstars; i++) {
+    stars[i].X      = smpdata[i].X;
+    stars[i].Y      = smpdata[i].Y;
+    stars[i].M      = smpdata[i].M;
+    stars[i].dM     = smpdata[i].dM;
+    stars[i].dophot = smpdata[i].dophot;
+
+    stars[i].Mgal   = smpdata[i].M;
+    stars[i].Map    = smpdata[i].dM;
+    stars[i].fx     = smpdata[i].fx;
+    stars[i].fy     = smpdata[i].fy;
+    stars[i].df     = smpdata[i].df;
+  }    
+  *nstars = Nstars;
+  return (stars);
+}
+
+Stars *ConvertPS1DATA (FTable *table, int *nstars) {
+
+  int i, Nstars;
+  Stars *stars;
+  PS1Data *ps1data;
+
+  ps1data = gfits_table_get_PS1Data (table, &Nstars, NULL);
+
+  ALLOCATE (stars, Stars, Nstars);
+  for (i = 0; i < Nstars; i++) {
+    stars[i].X      = ps1data[i].X;
+    stars[i].Y      = ps1data[i].Y;
+    stars[i].M      = ps1data[i].M;
+    stars[i].dM     = ps1data[i].dM;
+    stars[i].dophot = ps1data[i].dophot;
+
+    stars[i].Mgal   = ps1data[i].M;
+    stars[i].Map    = ps1data[i].dM;
+    stars[i].fx     = ps1data[i].fx;
+    stars[i].fy     = ps1data[i].fy;
+    stars[i].df     = ps1data[i].df;
+  }    
+  *nstars = Nstars;
+  return (stars);
+}
Index: /trunk/Ohana/src/addstar/src/ReadStarsTEXT.c
===================================================================
--- /trunk/Ohana/src/addstar/src/ReadStarsTEXT.c	(revision 10897)
+++ /trunk/Ohana/src/addstar/src/ReadStarsTEXT.c	(revision 10897)
@@ -0,0 +1,85 @@
+# include "addstar.h"
+# define D_NSTARS 1000
+# define BYTES_STAR 66
+# define BLOCK 1000
+
+Stars *ReadStarsTEXT (FILE *f, unsigned int *nstars) {
+
+  int j, N, Nextra, Ninstar, Nskip, Nbytes, nbytes;
+  int done;
+  char *buffer, *c, *c2;
+  double tmp;
+  Stars *stars;
+  
+  /* load in stars by blocks of 1000 */
+  N = 0;
+  ALLOCATE (buffer, char, (BLOCK*BYTES_STAR) + 1);
+  buffer[BLOCK*BYTES_STAR] = 0;
+  Nextra = 0;
+
+  ALLOCATE (stars, Stars, *nstars);
+
+  while (N < *nstars) {
+    /* load next data block */
+    Nbytes = BYTES_STAR * BLOCK - Nextra;
+    nbytes = fread (&buffer[Nextra], 1, Nbytes, f);
+    if (nbytes == 0) {
+      *nstars = N;
+      return (stars);
+    }
+    nbytes += Nextra;
+
+    /* check line-by-line integrity */
+    c = buffer;
+    done = FALSE;
+    while ((c < buffer + nbytes) && (!done)) { 
+      for (c2 = c; *c2 == '\n'; c2++);
+      if (c2 > c) { /* extra return chars */
+	memmove (c, c2, (int)(buffer + nbytes - c2));
+	Nskip = c2 - c;
+	nbytes -= Nskip;
+	bzero (buffer + nbytes, Nskip);
+	if (VERBOSE) fprintf (stderr, "deleted %d extra return chars\n", Nskip);
+      }
+      c2 = strchr (c, '\n');
+      if (c2 == (char *) NULL) {
+	done = TRUE;	
+	continue;
+      }
+      c2++;
+      if ((c2 - c) != BYTES_STAR) { /* bad line, delete it */
+	memmove (c, c2, (int)(buffer + nbytes - c2));
+	Nskip = c2 - c;
+	nbytes -= Nskip;
+	bzero (buffer + nbytes, Nskip);
+	if (VERBOSE) fprintf (stderr, "deleted line, %d extra chars\n", Nskip);
+      } else {
+	c = c2;
+      }
+    }
+
+    /* extract data for stars */
+    Ninstar = nbytes / BYTES_STAR;
+    Nextra = nbytes % BYTES_STAR;
+    for (j = 0; (j < Ninstar) && (N < *nstars); j++, N++) {
+      dparse (&stars[N].X,  1, &buffer[j*BYTES_STAR]);
+      dparse (&stars[N].Y,  2, &buffer[j*BYTES_STAR]);
+      dparse (&stars[N].M,  3, &buffer[j*BYTES_STAR]);
+
+      /* cmp files carry dM in millimags */
+      dparse (&tmp, 4, &buffer[j*BYTES_STAR]);
+      stars[N].dM = 0.001*tmp;
+
+      dparse (&tmp,         5, &buffer[j*BYTES_STAR]);
+      stars[N].dophot = tmp;
+
+      dparse (&stars[N].Mgal, 7, &buffer[j*BYTES_STAR]);
+      dparse (&stars[N].Map,  8, &buffer[j*BYTES_STAR]);
+      dparse (&stars[N].fx,   9, &buffer[j*BYTES_STAR]);
+      dparse (&stars[N].fy,  10, &buffer[j*BYTES_STAR]);
+      dparse (&stars[N].df,  11, &buffer[j*BYTES_STAR]);
+    }
+  }
+  *nstars = N;
+  return (stars);
+}
Index: /trunk/Ohana/src/addstar/src/gstars.c
===================================================================
--- /trunk/Ohana/src/addstar/src/gstars.c	(revision 10896)
+++ /trunk/Ohana/src/addstar/src/gstars.c	(revision 10897)
@@ -1,25 +1,20 @@
 # include "addstar.h"
 
-/* break into sections, loop over images and load stars one at a time */
+enum {NONE, SIMPLE_CMP, SIMPLE_CMF, SIMPLE_MEF, MOSAIC_CMP, MOSAIC_CMF, MOSAIC_MEF, MOSAIC_PHU};
+/* note: MEF implies CMF */
 
-Stars *gstars (char *file, int *NSTARS, int photcode, Image **imageList, int Nimages) {
+Stars *gstars (char *filename, int *NSTARS, int photcode, Image **imageList, int Nimages) {
 
+  int i, Nfile, Nimage;
+  char **file;
   FILE *f;
-  int i, j, N, Nbytes, naxis;
-  int itmp, hour, min, simple;
-  char *name, *c, photname[64], line[80];
-  double tmp, sec, dMs;
-  Stars *stars, *rdstars;
-  Header header;
-  Image *images;
+  glob_t globList;
+  char **exthead, **extdata;
+  int *extnum_head, *extnum_data, *extsize;
+  Header *header, **headers;
+  Image *image;
+  Stars *inStars, *stars;
 
-  /* XXX which of these values are required, which are just interesting? */
-
-  # if (0)
-  /* file may be a glob, in which case it should represent multiple chip files
-     from a single mosaic image */
-
-  // parse the word as a glob
-  glob_t globList;
+  // parse the filename as a glob
   globList.gl_offs = 0;
   glob (file, 0, NULL, &globList);
@@ -28,123 +23,163 @@
   // otherwise save all glob matches
   if (globList.gl_pathc == 0) {
-    if (!gfits_read_header (file, &header)) {
-      fprintf (stderr, "ERROR: can't read header for %s\n", file);
-      exit (1);
+    Nfile = 1;
+    ALLOCATE (file, char *, Nfile);
+    file[0] = strcreate (filename);
+  } else {
+    Nfile = globList.gl_pathc;
+    ALLOCATE (file, char *, Nfile);
+    for (i = 0; i < Nfile; i++) {
+      file[i] = strcreate (globList.gl_pathv[i]);
+    }
+  }
+
+  // open the first file, read the PHU header
+  f = fopen (file[0], "r");
+  if (f == NULL) {
+    fprintf (stderr, "ERROR: can't read header for %s\n", file[0]);
+    exit (1);
+  }
+  ALLOCATE (header, Header, 1);
+  gfits_fread_header (f, header);
+
+  mode = GetFileMode (header);
+
+  stars = NULL;
+  Nstars = 0;
+
+  // valid for all but MEF files
+  if ((mode != SIMPLE_MEF) && (mode != MOSAIC_MEF)) {
+    Nimage = Nfile;
+    ALLOCATE (image, Image, Nimage);
+    for (i = 0; i < Nfile; i++) {
+      if (i > 0) {
+	f = fopen (file[i], "r");
+	if (f == NULL) {
+	  fprintf (stderr, "ERROR: can't read header for %s\n", file[i]);
+	  exit (1);
+	}
+	gfits_fread_header (f, header);
+      }
+
+      ReadImageHeader (header, &image[i]);
+    
+      switch (mode) {
+	case SIMPLE_CMP:
+	case MOSAIC_CMP:
+	  inStars = ReadStarsTEXT (f, header, NULL, &image[i].nstars);
+	  inStars = FilterStars (inStars, &image[i]);
+	  stars = MergeStars (inStars, NinStars, stars, &Nstars);
+	  break;
+
+	case SIMPLE_CMF:
+	case MOSAIC_CMF:
+	  inStars = ReadStarsFITS (f, header, NULL, &image[i].nstars);
+	  inStars = FilterStars (inStars, &image[i]);
+	  stars = MergeStars (inStars, NinStars, stars, &Nstars);
+	  break;
+
+	case MOSAIC_PHU:
+	  inStars = NULL;
+	  NinStars = 0;
+	  break;
+      }
+      fclose (f);
+      gfits_free_header (header);
     }
   } else {
-    for (i = 0; i < globList.gl_pathc; i++) {
-      if (!gfits_read_header (globList.gl_pathv[i], &header)) {
-	fprintf (stderr, "ERROR: can't read header for %s\n", globList.gl_pathv[i]);
-	exit (1);
+    
+    /* we need to examine the extensions to determine the headers and the data */
+    NHEADER = 10;
+    ALLOCATE (headers, Header *, NHEADER);
+
+    // the first header is already loaded
+    headers[0] = header;
+    Nskip = gfits_matrix_size (header);
+    fseek (f, Nskip, SEEK_CUR); 
+
+    // load all headers into memory
+    done = FALSE;
+    for (i = 1; !done; i++) {
+      ALLOCATE (header, Header, 1);
+      
+      status = gfits_fread_header (f, header);
+      if (!status) { 
+	done = TRUE;
+      } else {
+	headers[i] = header;
+	Nskip = gfits_matrix_size (header);
+	fseek (f, Nskip, SEEK_CUR); 
+      }
+      if (i == NHEADER - 1) {
+	NHEADER += 10;
+	REALLOCATE (headers, Header *, NHEADER);
       }
     }
-  }
-  # endif
+    Nheader = i;
+    
+    // space to store the images, indexes to the matching headers
+    Nimage = 0;
+    NIMAGE = Nheader;
+    ALLOCATE (image, Image, NIMAGE);
+    ALLOCATE (exthead, char **, NIMAGE);
+    ALLOCATE (extdata, char **, NIMAGE);
+    ALLOCATE (extnum_head, int, NIMAGE);
+    ALLOCATE (extnum_data, int, NIMAGE);
+    ALLOCATE (extsize, int, Nheader);
 
-  if (!gfits_read_header (file, &header)) {
-    fprintf (stderr, "ERROR: can't read header for %s\n", file);
-    exit (1);
-  }
+    if (mode == MOSAIC_MEF) {
+      exthead[Nimage] = strcreate ("PHU");
+      extdata[Nimage] = strcreate ("NONE");
+      extnum_data[Nimage] = -1;
+      extnum_head[Nimage] = 0;
+      Nimage ++;
+    }
 
-  /* find image rootname */
-  /* if we are reading a MEF mosaic image, attached [extname] to the name of the image */
-  name = filebasename (file);
-  strcpy (image[0].name, name);
-  free (name);
+    // now examine the headers, count the table entries, find corresponding headers
+    for (i = 0; i < Nheader; i++) {
+      extsize[i] = headers[i][0].size + gfits_matrix_size (headers[i]);
+      gfits_scan (headers[i], "EXTTYPE", "%s", 1, tmpword);
+      if (!strcmp (tmpword, "SMPDATA") ||  
+	  !strcmp (tmpword, "PS1DATA")) {
+	exttype[Nimage] = strcreate (tmpword);
+	gfits_scan (headers[i], "EXTNAME", "%s", 1, tmpword);
+	extdata[Nimage] = strcreate (tmpword);
+	gfits_scan (header[i], "EXTHEAD", "%s", 1, tmpword);
+	exthead[Nimage] = strcreate (tmpword);
+	extnum_data[Nimage] = i;
+	extnum_head[Nimage] = -1;
+	// find the matching exthead entry
+	for (j = 0; j < Nheader; j++) {
+	  gfits_scan (header[j], "EXTNAME", "%s", 1, tmpword);
+	  if (!strcmp (tmpword, exthead[Nimage])) {
+	    extnum_head[Nimage] = j;
+	  }
+	}
+	// skip or crash on table with missing matching header?
+	if (extnum_head[Nimage] == -1) {
+	  fprintf (stderr, "ERROR: can't read header for %s\n", file[0]);
+	  exit (1);
+	}
+	Nimage ++;
+      }
+    }
 
-  /* this is really lame.  get sky background from Flips data */
-  /*** the problems are: 
-       1) need to define sky entry in image table
-       2) need to have sky measurement (in reg.db) 
-       (could be derived from star data?)
-       
-    XXX EAM : make a different SKYPROBE format with sky in appropriate place? 
-  ***/
+    // now run through the images, interpret the headers and read the stars
+    for (i = 0; i < Nimage; i++) {
+      Nhead = extnum_head[i];
+      ReadImageHeader (header[Nhead], &image[i]);
 
-  /** I'm essentially using subpix to be synonymous with SKYPROBE */
-  if (SUBPIX) {
-    char *p;
-    int sky;
-    
-    sky = 0;
-    gfits_scan (&header, "HISTORY", "%S", 10, line);
-    p = strstr (line, "|I=");
-    if (p != (char *) NULL) {
-      p += 3;
-      sky = atoi (p);
+      if (!strcmp(exthead[i], "PHU")) continue;
+
+      // advance the pointer to the start of the corresponding table block
+      Ndata = extnum_data[i];
+      Nskip = 0;
+      for (j = 0; j < Ndata; j++) {
+	Nskip += extsize[j];
+      }
+      fseek (f, Nskip, SEEK_SET); 
+	 
+      inStars = ReadStarsFITS (f, header[Nhead], header[Ndata], &NinStars);
+      inStars = FilterStars (inStars, &NinStars, &image[i]);
+      stars = MergeStars (inStars, NinStars, stars, &Nstars);
     }
-    image[0].Myyyy = MIN (sky - 0x8000, 0xffff);
-  }
-
-  /* read in data section */
-  f = fopen (file, "r");
-  if (f == NULL) {
-    fprintf (stderr, "ERROR: can't read data from %s\n", file);
-    exit (1);
-  }
-  fseek (f, header.size, SEEK_SET); 
-
-  /* don't try to read stars for mosaic */
-  if (!strcmp (&image[0].coords.ctype[4], "-DIS")) {
-    rdstars = NULL;
-    image[0].nstar = 0;
-  } else {
-    /* read from FITS table or from text table */
-    naxis = 2;
-    gfits_scan (&header, "NAXIS",  "%d", 1, &naxis);
-    gfits_scan (&header, "SIMPLE", "%t", 1, &simple);
-    if ((naxis == 0) && !TEXTMODE && simple) {
-      Nbytes = gfits_matrix_size (&header);
-      fseek (f, Nbytes, SEEK_CUR); 
-      rdstars = rfits (f, &image[0].nstar);
-      if (rdstars == NULL) {
-	fprintf (stderr, "ERROR: failed to read fits table\n");
-	exit (1);
-      }
-    } else {
-      /* allocate space for stars */
-      if (!gfits_scan (&header, "NSTARS", "%d", 1, &image[0].nstar)) {
-	fprintf (stderr, "ERROR: failed to find NSTARS\n");
-	exit (1);
-      }
-      rdstars = rtext (f, &image[0].nstar);
-    }
-  }
-  fclose (f);
-
-  /* modify resulting star list */
-  ALLOCATE (stars, Stars, image[0].nstar);
-  for (N = j = 0; j < image[0].nstar; j++) {
-    /* allow for some dynamic filtering of star list */
-    if (SNLIMIT && rdstars[j].dM > SNLIMIT) continue;
-    if (XMAX && (rdstars[j].X > XMAX)) continue;
-    if (XMIN && (rdstars[j].X < XMIN)) continue;
-    if (YMAX && (rdstars[j].Y > YMAX)) continue;
-    if (YMIN && (rdstars[j].Y < YMIN)) continue;
-    stars[N] = rdstars[j];
-
-    XY_to_RD (&stars[N].R, &stars[N].D, stars[N].X, stars[N].Y, &image[0].coords);
-    while (stars[N].R <    0.0) stars[N].R += 360.0;
-    while (stars[N].R >= 360.0) stars[N].R -= 360.0;
-    stars[N].found = -1;
-    stars[N].code = photcode;
-
-    if (SUBPIX) {
-      dMs = get_subpix (stars[N].X, stars[N].Y);
-      stars[N].M    -= dMs;
-      stars[N].Mgal -= dMs;
-      stars[N].Map  -= dMs;
-      dMs = scat_subpix (stars[N].X, stars[N].Y);
-      stars[N].dM = hypot (stars[N].dM, dMs);
-    }
-    N ++;
-  }
-  image[0].nstar = N;
-  REALLOCATE (stars, Stars, image[0].nstar);
-  free (rdstars);
-
-  if (VERBOSE) fprintf (stderr, "read %d stars from target file\n", image[0].nstar);
-  *NSTARS = image[0].nstar;
-
-  return (stars);
-}