proj/mag/harp/apps/ingest_mharp.c

/*
 *  ingest_mharp.c
 *
 * Purpose:
 *  Ingest HARPs that are described in files underneath a
 *  given root directory.  The HARPs ingested are listed by 
 *  number in a set of files ("lists" argument).  Their metadata,
 *  including lat/lon bounding boxes, some space weather keywords,
 *  and bitmap outlines, are present in files and directories 
 *  underneath "root", in locations marked by the harp number 
 *  appearing in the list files.
 *  This set of files is generated by a Matlab tracker, which 
 *  is typically run by a shell script before running this
 *  ingestor.  
 *  There is also a higher-level coordination script which 
 *  stages data, runs the tracker, this ingestor, and other 
 *  programs needed for complete processing.
 *  Thus, you should not generally need to run this ingestor
 *  manually.
 * 
 * Input File Structure:
 *  A directory structure begining at the supplied "root"
 *  ("root" typically ends in Tracks/jsoc), containing:
 *    -- text file(s), named in the "lists" argument to
 *       this routine, containing HARP ID numbers.
 *  For each HARP ID HHHHHH, a directory root/track-HHHHHH, 
 *  containing:
 *    -- a coordinate-frame file with image geometry
 *    -- a statistics summary file containing region statistics
 *  Each of these two files contains one line per appearance of
 *  the HARP (each T_REC).
 *  The same directory also contains:
 *    -- a series of region-mask FITS files, with bitmap masks 
 *       for the on-region pixels, one file for each HARP 
 *       appearance.
 *
 * Output:
 *  The HARP IDs (HHHHHH) named in each of the files given
 *  in "lists" are ingested into JSOC.  For each HARP, this 
 *  includes all appearances (each T_REC), plus a given 
 *  temporal padding (tpad, in units of images)
 *  before/after the HARP.  For NRT processing, you can 
 *  use the trec argument to only ingest a given number of 
 *  the most recent records (use the default of zero to ingest
 *  all available records).
 *
 * Arguments:
 *  root       -- Input file directory, typ. ends in Track/jsoc
 *  lists [track-new.txt] -- files in "root" listing HARP numbers
 *  match [0]  -- +1/-1 to dump match-info to "root", -1 to exit after dump
 *  out        -- Output HARP data series
 *  mag [hmi.M_720s] -- Input magnetogram data series, alter for NRT
 *  mask       -- Input mask data series, e.g. hmi.Marmask_720s
 *  href       -- Reference HARP set for box geometry, default none
 *  t0cut      -- Earliest time allowed to ingest (default: none)
 *  t1cut      -- Latest time allowed to ingest (default: none)
 *  trec [0]   -- Number of T_RECs to ingest (0=all, 1=latest)
 *  tpad [120] -- Temporal padding (in images, >=0, 120=1 day)
 *  readonly[0]-- No status-update writes into root dir
 *  verb [1]   -- Verbosity: 0=errors only; 1, 2, 3 = more
 *
 * Usage:
 *  For monthly processing, only the finalized new tracks are
 *  to be ingested, because others are still subject to change.
 *  The default temporal padding of 120 (images) is standard:
 *  ingest_mharp root=/tmp22/tracking/2011/Tracks/jsoc   \
 *               out=hmi.mharp_720s                      \
 *       lists=track-new.txt                     \
 *               mask=hmi.Marmask_720s
 *
 *  For NRT processing, both new and pending tracks are ingested,
 *  but only the most recent record should be used.  Also, no 
 *  temporal padding should be used:
 *  ingest_mharp root=/tmp22/tracking/2011-nrt/Tracks/jsoc \
 *               out=hmi.mharp_720s_nrt                    \
 *       lists=track-new.txt,track-pending.txt     \
 *               mag=hmi.M_720s_nrt                        \
 *               mask=hmi.Marmask_720s_nrt                 \
 *               trec=1 tpad=0
 *
 *  Based on M. Turmon's tracking algorithm and ingest_dsds_a.c
 *
 *  Written by X. Sun, M. Turmon
 *
 *  Version history:
 *  v0.0    Feb 08 2011
 *  v0.1    Feb 26 2011
 *  v0.2    Apr 17 2011
 *  v0.5    Jun 28 2011
 *      v0.6    Aug 2011
 *      v0.7    Oct 2011
 *      v1.0    In production, March 2012
 *      v1.1    Feb-Apr 2013
 *
 *  Notes:
 *  -- In general, search for FIXME and TODO for loose ends
 *  
 *  -- v0.0
 *     patch-stats.txt not imported
 *     patch-status.txt not modified
 *  -- v0.1
 *     Changed fits reader
 *     Added file parser
 *  -- v0.2
 *     Put all data into the frame structure
 *  -- v0.5
 *     Updated SW keywords and other keywords
 *     Major restructuring:
 *     Compiling a master matrix first, with HARPNUM and T_REC as indices
 *     Ingesting according to T_REC, to minimize I/O of mask series
 *  -- v0.6
 *     fixed crash if no HARPs were present
 *     support for NRT processing: multiple input harp-list files, can specify
 *     no temporal padding and only ingest most recent T_REC
 *     patch for correct bounding-boxes for HARPS that span the equator
 *     more thorough error checking
 *     better logging and status output
 *     cleanups including format and consistent variable naming
 *     documentation and usage
 *  -- v0.7
 *     can match against NOAA ARs in a data series, output match info,
 *     read tracker parameters from root dir
 *  -- v1.0
 *     released into production
 *  -- v1.1
 *     Improved HARP/NOAA matching in case of large HARPs which match several ARs.
 *     Reorder code to load WCS just once, at start, to identify error conditions
 *     and farside passage so that metadata like T_FRST will always be correct.
 *     Add concept of ingestion filtering so that we know what portion of each 
 *     HARP will be ingested, before starting ingestion, again for metadata 
 *     correctness.
 *     Add ability to use an existing HARP series for box-size information,
 *     rather than using box-sizes in the tracker files.  Reason is to permit
 *     patching of gaps in HARP series; ingest tracker bitmaps but use existing
 *     box-sizes for compatibility with already-ingested HARP geometry.
 *  -- Added readonly mode on Monica's suggestion
 *
 */


#include "jsoc_main.h"
#include "drms_types.h"
#include <time.h>
#include <math.h>
#include <timeio.h>
#include "fitsio.h"
// lat-lon to image and back
#include "cartography.c"
// patch statistics
#include "roi_stats_mag_defs.h"
// Keyword propagation
#include "propagate_keys.c"


// allow V_printf() to use the verbosity status from anywhere
static int verbflag;
// log file pointers (initialized to NULL)
static FILE *LOGout = NULL;
static FILE *LOGerr = NULL;

// standard string length
#define STR_MAX 256
// this is one more than the max number of list files in the "lists" argument
#define LIST_MAX 10
// standard file names
// log file names (output from this code)
#define FN_INGEST_STDLOG "ingest-latest.log"
#define FN_INGEST_STDERR "ingest-latest.err"
// per-run parameter info
#define FN_TRACK_PARAM "track-param.txt"
// per-track coordinate frame info
#define FN_TRACK_FRAME "track-frame.txt"
// per-track space weather keywords
#define FN_TRACK_STATS "track-stats.txt"
// per-track status summary
#define FN_TRACK_STATUS "track-status.txt"
// per-track subdirectory name
#define FN_TRACK_DIR   "track-%06d"
// bitmap outlines, one for each appearance of each track
#define FN_TRACK_FITS  "patch-%06d.fits.gz"
// used for unspecified input arguments
#define kNOT_SPEC   "Not Specified"
// below this number of seconds, two records were taken at the same time
#define REC_EPS 10.0
// seconds per day
#define SEC_PER_DAY (24.0*60.0*60.0)

/************************************************************* 
 *
 * Error handling/logging macros
 *
 *************************************************************
 */

// V_printf: facilitate verbflag output
//   if flag is > 0, output is to stdout+LOGout, if < 0, to stderr+LOGerr
//   if flag is 0, no output is made at all
// If LOGout/LOGerr are NULL, they are skipped.
// The message is printed in the form:
//   <module_name>: <preamble><message>  (if preamble is non-NULL)
//   <message>                           (if preamble is NULL)
// Usage:
//   V_printf(verbflag > 0, "\t", "Mask(%d) = %d\n", 2048, mask[2048]);
void
V_printf(int flag, char *preamble, char *format, ...) {
  va_list args;
  extern char *module_name;
  extern FILE *LOGout;  // NULL OK to suppress
  extern FILE *LOGerr;  // NULL OK to suppress
  FILE *fp = (flag > 0) ? stdout : stderr;
  FILE *lp = (flag > 0) ? LOGout : LOGerr;

  if (flag != 0) {
    // print to stream from stdio
    va_start(args, format);
    if (preamble)
      fprintf(fp, "%s: %s", module_name, preamble);
    vfprintf(fp, format, args); // va_list version of printf
    fflush(fp);
    va_end(args);
    // optionally print to log file
    if (lp) { 
      va_start(args, format);
      if (preamble) 
    fprintf(lp, "%s: %s", module_name, preamble);
      vfprintf(lp, format, args);
      fflush(lp);
      va_end(args);
    }
  } else {
    va_start(args, format);
    // (do nothing)
    va_end(args);
  }
}
// exit with error
//  (standard trick to swallow the semicolon)
// NOTE: return() comes back from DoIt, so only call within DoIt()
#define DIE(...) do { \
    fflush(stdout); \
        V_printf(-1, "", "FATAL: %s.\n", __VA_ARGS__); \
    if (LOGout) fclose(LOGout); \
    if (LOGerr) fclose(LOGerr); \
        return 1; \
        } while (0)
// report non-fatal error
//  (standard trick to swallow the semicolon)
#define WARN(...) do { \
    fflush(stdout); \
        V_printf(-1, "", "WARNING: %s. Continuing.\n", __VA_ARGS__); \
        } while (0)


#define RADSINDEG       (M_PI/180.0)
#define RAD2ARCSEC      (648000. / M_PI)

// Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, 
// sina,cosa = sin and cos of CROTA2 resp.
// and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
// PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun 
// from disk center.
// -- These are copy-and-paste macros from JSOC code --
#define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
#define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
#define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
#define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
// instantiate local variables from a wcs_t structure
//   (note, some small changes in variable naming to suit JSOC code)
#define WCS2LOCALS(wcs) crvalx=(wcs).crval1; \
                crvaly=(wcs).crval2; \
                crpix1=(wcs).crpix1; \
                crpix2=(wcs).crpix2; \
                cdelt =(wcs).cdelt1; \
                crota2=(wcs).crota2; \
                        sina=sin(crota2 * RADSINDEG); \
            cosa=cos(crota2 * RADSINDEG)

/***********************************************************************
 *
 * Declarations
 *
 ***********************************************************************/

/*
 * WCS, in the same parameterization as in HMI keywords
 */
typedef struct {
  int ok;           // 0 if not OK
  char missing[32]; // short message indicating missing data, if !ok
  // types here are double or float according to hmi.M_720s types
  double rsun_ref; // RSUN_REF || 6.96e8
  double dsun_obs; // DSUN_OBS
  float  cdelt1;   // CDELT1
  float  crval1;   // CRVAL1 -- disc center, arcsec
  float  crval2;   // CRVAL2
  float  crpix1;   // CRPIX1 -- center, ccd, orig.
  float  crpix2;   // CRPIX2
  float  crota2;   // CROTA2
  float  crlt_obs; // CRLT_OBS (converted to radians)
  float  crln_obs; // CRLN_OBS (converted to radians)
} wcs_t;

/* 
 * Status of a patch_info_t: initialized to "invalid" upon creation
 */
typedef enum {
  Patch_Invalid = 0,  // uninitialized fields in data structure
  Patch_Normal,       // ordinary, observed image
  Patch_Padding,      // patch is temporal padding: stats not valid, lat/lon BB is
                      // padding occurs at end, and for placeholder ROIs in middle
} patch_status_t;

/*
 * Discrete tag indicating patch type that is read from a file
 *
 * Currently the frame file can have tag = 0, 1, or 2.
 * The tracker internals can have (-1) which is for a placeholder patch,
 * but these patches are not written out by the tracker.  We re-introduce
 * it as Tag_Faint to indicate the patch was too faint to see (even after 
 * trying harder).
 */
typedef enum {
  // a special value we introduce for patches not read from a file
  //   (typically padding patches)
  Patch_Tag_None   = -99,   // no tag was read for this patch
  // not read from file, but inserted by us while interpolating HARP
  Patch_Tag_Faint  =  -1,   // placeholder patch (between t0 and t1 of HARP)
  // these are as read from files and defined in Matlab
  Patch_Tag_Harder =   0,   // tried harder and matched
  Patch_Tag_Normal =   1,   // ordinary, observed image
  Patch_Tag_Merge  =   2,   // tracks were merged on this frame
} patch_tag_t;

/*
 * To-be-ingested status of each patch
 */
typedef enum {
  Patch_Ingest_Eligible = 0, // will ingest
  Patch_Ingest_Invalid,      // no patch for this harp at this t_rec
  Patch_Ingest_Gap,          // data gap in mags
  Patch_Ingest_Missing,      // data missing in masks
  Patch_Ingest_Omit,         // omit ingesting due to passed options
  Patch_Ingest_TooShort,     // harp too short, typ. cosmic ray
  Patch_Ingest_NUM,          // (number of distinct values)
} patch_ingest_t;

static const char* const 
Patch_Ingest_Names[Patch_Ingest_NUM] = {
  "eligible to ingest",
  "no patch for this harp at this t_rec [ok]",
  "data gap in mags [ok]",
  "data missing in masks [likely error]",
  "omitted ingesting due to supplied options [ok]",
  "too short, typ. cosmic ray [ok]",
};

/* 
 * Information for one patch: one HARP at one T_REC
 */
typedef struct {
  patch_status_t valid;            // see enum above
  int num;                         // HARP id number
  patch_tag_t tag;                 // tag: patch category from frame file
  patch_ingest_t ingest;           // to-be-ingested status
  int success;                     // successful ingestion?
  int x0, y0;                      // origin of HARP bitmap in FD image (0-based)
  int fits_nx, fits_ny;        // HARP bitmap size from frame file
  float lat0, lon0, lat1, lon1;    // bounding box from track-frame.txt
  float omega;                     // HARP-specific angular rotation (deg/day)
  float stats[RS_num_stats];       // statistics from track-stats.txt
  char *patchName;             // bitmap image file path -- strdup'd
  char *image;                 // bitmap HARP image itself
  int dims[2];                     // dimension of above image
  int xmin, xmax, ymin, ymax;      // calculated bounding box (0-based)
} patch_info_t;

/* 
 * Information for a HARP (spans all patches, that is, all T_REC's)
 */
// fits string = 68 chars; NOAA ID has 5, comma is 1, so 68/(5+1) = 11 NOAA IDs max
#define HARP_INFO_MATCH_MAX 11
typedef struct {
  int id;            // HARP numeric id
  int rec0, rec1;    // HARP was tracked in rec0..rec1 (index into trec array)
  int rec0p, rec1p;  // HARP ingested during rec0p..rec1p (with pad)
  int n_missing;     // number of missing-image gaps in rec0..rec1
  int n_patch;       // number of valid patches (Normal+Padding) in rec0..rec1
  int n_patchp;      // number of valid patches (Normal+Padding) in rec0p..rec1p
  int n_patch_eligible;// number of patches eligible to ingest
  int n_patch_ingested;// number of patches successfully ingested
  int n_patch_err;   // number of patches with ingestion errors
  int top_match;     // best match (index into MarpInfo array)
  int n_match;       // number of matches below (0 for none)
  int match[HARP_INFO_MATCH_MAX];
                     // list of matching MARP's (index into MarpInfo array)
} harp_info_t;

/* 
 * Information for a T_REC (spans all HARPs)
 */
typedef enum {
  Trec_Data_Unknown = 0, // have not looked at data presence
  Trec_Data_Gap,         // no mag present (data gap -- assume no mask either)
  Trec_Data_Missing,     // mag present, but no mask (should be very rare)
  Trec_Data_OK,          // mag and mask both present
  Trec_Data_NUM,         // (number of distinct values)
} trec_data_t;

static const char* const 
Trec_Data_Names[Trec_Data_NUM] = {
  "un-inspected [ok: data record not examined]",
  "no mag present [ok: data gap]",
  "mag present, but no mask [likely error: missing data]",
  "mag and mask both present [ok]",
};

typedef struct {
  TIME t;            // T_REC as a number
  char str[24];      // T_REC as a string
  trec_data_t data;  // is (mag,mask) data present and readable?
  wcs_t wcs;         // wcs info, directly from DRMS keywords
  int nharp;         // number of HARP patches to be ingested at this time
} trec_info_t;

/*
 * singletons for ingestion results
 */
typedef struct {
  int try_ok;
  int try_err;
  int hist[Patch_Ingest_NUM];
} patch_ingest_tab_t;

typedef struct {
  int trec_patch_err;
  int trec_patch_valid;
  int hist[Trec_Data_NUM];
} trec_ingest_tab_t;

/*
 * Whole-run info
 * (run start/end, tracker parameters)
 */
typedef struct {
  int n_par;              // number of parameters
  char **par_names;       // parameter names
  char **par_vals;        // parameter values
} run_info_t;

/* 
 * Information for a match AR (spans all patches, that is, all T_REC's)
 *
 * (rec0, rec1 not set now)
 */
// no need for innovation here
typedef struct {
  int id;            // NOAA AR numeric id
  int rec0, rec1;    // AR tracked in rec0..rec1 (index into trec array) -- unset!
  int top_match;     // best match (index into HarpInfo array)
  int n_match;       // number of times MARP was matched to a HARP
} marp_info_t;

/*
 * Used in a decorate-sort-undecorate idiom for qsort to rank matches
 */
typedef struct  {
  int id;         // NOAA AR id
  int marp_index; // index into array of marp_info's
  double count;
  double size;
} marp_sort_t;

/* 
 * Information for one match-patch: one (NOAA) AR at one T_REC
 */
typedef struct {
  patch_status_t valid;  // enum above
  int id;                // not needed, but makes life easier
  TIME t;                // time
  float lat, lon;        // center lat/lon (heliographic)
  float lon_carr;        // carrington longitude
  float lon_wid;         // width in degrees longitude
  float area;            // AR area microhemispheres
  int spot;              // spot count
  char class_mag[32];    // magnetic classification (Alpha, etc. -- truncated if need be)
  char class_zur[4];     // zurich classification (3 chars -- truncated otherwise)
} match_info_t;

static const char *tracker_run_info_keys[] = {
  // the most important per-run params
  "TKP_KWID", "kwid",
  "TKP_KLAT", "klat",
  "TKP_TAU",  "tau",
  "TKP_TAU2", "tau2",
  "TKP_ACTV", "active",
  "TKP_FNUM", "final_num",
  "TKP_FTIM", "final_time",
  "TKP_MAPR", "maprate",
  // other per-run parameters
  "TKP_RUNN", "run_name",
  // this is present in params, but not currently in JSD, and is unimportant
  // "TKP_RUNC", "run_number", 
  "TKP_RUNT", "run_time",
  // end sentinel -- required
  NULL, NULL,
};

/***********************************************************************
 *
 * Function signatures
 *
 * load_x:  load x from a file
 * set_x:   set up x from already-loaded data
 *
 ***********************************************************************/

// load tracker params from text file
int load_tracker_params(const char *rootDir, run_info_t *ri);
// load master index array info from file
int load_run_info(const char *rootDir, char **, double, int *nHarp, int *nRec_all, 
          TIME *, TIME *);
// get HARP ID portion of harp info array
int load_harp_ids(const char *rootDir, char **, harp_info_t *harpInfo);
// main info parser/loader
int load_all_patch_info(const char *rootDir, const char *refHarp,
            patch_info_t *, harp_info_t *, trec_info_t *, 
            int nHarp, int, int, double, int);
// Frame file parser, one line per time, return 0 if successful
int load_frame(FILE *fp, patch_info_t *pInfo, TIME *t);
// Stats file parser, one line per time, return 0 if successful
int load_stats(FILE *fp, float *stats);
// mark harp as ingested within disk file
int mark_ingested_harp(const char *, run_info_t *, harp_info_t *);
// print one patch to file
void patch_print(FILE *fp, char *s, trec_info_t *tI, patch_info_t *pI);
// set up fields in trec info array
void set_trec_times(trec_info_t *tRec, TIME t0, int nRec_all, int nRec_pad, double cadence);
// set up has-data field in trec info array
void set_trec_hasdata(trec_info_t *tRec, char *magQuery, char *maskQuery, int *npix, int nRec_all);
// mark patches as ingestable
void mark_eligible_patches(patch_info_t *, harp_info_t *, trec_info_t  *tRec, 
               int nHarp, int nRec_all,  int nRec_eat, 
               int rec0cut, int rec1cut, int nRec_min);
// ingest one entire HARP (all records)
int ingest_harp(patch_info_t *, harp_info_t *, trec_info_t *trec1, trec_info_t *trec,
        run_info_t *,
        marp_info_t *,
        DRMS_Record_t *magRec, DRMS_Record_t *maskRec, 
        char *maskImg, int mask_stride, char *outQuery);
// ingest one record (T_REC) of one HARP
int ingest_record(patch_info_t *, harp_info_t *, trec_info_t *trec1, trec_info_t *trec, 
          run_info_t *,
          marp_info_t *,
          DRMS_Record_t *magRec, DRMS_Record_t *maskRec, char *outQuery);
// tabulate results of ingestion
int tabulate_patches(patch_ingest_tab_t *, trec_ingest_tab_t *, 
             patch_info_t *, harp_info_t *, trec_info_t *, int, int);
// Update bitmap image
int update_bitmap(patch_info_t *pInfo, char *maskImg, int maskStride);
// Keywords
char *get_instrument_mode(DRMS_Record_t *);
int set_keys_code_info(char *, DRMS_Record_t *);
int set_keys_runinfo(DRMS_Record_t *rec, run_info_t *runInfo);
int set_keys_stats(char *inst_mode, DRMS_Record_t *rec, DRMS_Record_t *magRec, DRMS_Record_t *maskRec, 
           patch_info_t *, harp_info_t *, trec_info_t *);
int set_keys_match(DRMS_Record_t *rec, harp_info_t  *hInfo, marp_info_t  *mInfo);
// convert DRMS "TIME" to integer t_rec offset
int time_to_recnum(TIME t, TIME t0, int rec0, int nRec, double);
// Find pixel bounding box from lat-lon box
int compute_boundary(patch_info_t *, int, wcs_t);
// Distance of one point from disk center
void distance2center(wcs_t wcs, double x, double y, double *xDist, double *yDist);
// load drms record, segment, wcs
char *load_drms_rec_seg(DRMS_RecordSet_t **, DRMS_Record_t **, DRMS_Array_t **, 
            wcs_t *, char *query, char *trec, char *segname, DRMS_Type_t);
// Ad-hoc FITS reader
int load_fits(char *filename, char **image, int *dims);


/************************************************************* 
 *
 * Timing
 * from T. Larson
 *
 *************************************************************
 */
static double getwalltime(void)
{
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return tv.tv_sec * 1000.0 + tv.tv_usec/1000.0;
}

static double getcputime(double *utime, double *stime)
{
  struct rusage ru;
  getrusage(RUSAGE_SELF, &ru);
  *utime = ru.ru_utime.tv_sec * 1000.0 + ru.ru_utime.tv_usec / 1000.0;
  *stime = ru.ru_stime.tv_sec * 1000.0 + ru.ru_stime.tv_usec / 1000.0;
  return *utime + *stime;
}

/***********************************************************************
 *
 * Match ARs: NOAA AR Loader
 *
 ***********************************************************************/

// this is currently just loading info, not matching
#include "mharp_match.c"

/*
 * export harp/marp match distance matrix to files for inspection
 *
 * if rootdir is NULL, skip the exports
 */
static
int
match_harp_export(const char *rootDir,
          double *dist,
          match_info_t *matchInfo,
          marp_info_t  *marpInfo,
          int nMarp, 
          patch_info_t *patchInfo,
          harp_info_t *harpInfo, 
          trec_info_t *tRec,
          int nRec_all,
          int nRec_min,
          int nHarp)
{
  const int print_if_invalid = 1;  // print a line even if patch is invalid
  FILE *fp, *fp2;
  const char *base = "Param";      // base for filenames of parameters
  const char *extn = "dat";        // extension for above
  char fn[STR_MAX];
  char root[STR_MAX];
  int rv;                          // return value
  int h, m;
  int i, inx, inx1;

  if (!rootDir)
    return 0; // skip it
  snprintf(root, sizeof(root), "%s/%s", rootDir, "Match");
  if (mkdir(root, 0775) != 0) {
    if (errno == EEXIST) {
      V_printf(-1, "", "Match dump: `%s' already exists, skipping dump.\n", root);
    } else {
      V_printf(-1, "", "Match dump: could not create `%s'.\n", root);
      perror(root);
      return 1;
    }    
  } else {
    V_printf(verbflag > 0, "", "Match dump: dumping to `%s'.\n", root);
  }
  /*
   *  HARP/MARP summary info
   */ 
  // distance
  snprintf(fn, sizeof(fn), "%s/%s-dist.%s", root, base, extn);
  if ((fp = fopen(fn, "w")) == NULL) return 1;
  for (h = 0; h < nHarp; h++) {
    for (m = 0; m < nMarp; m++)
      fprintf(fp, "%g ", dist[h*nMarp+m]);
    fprintf(fp, "\n");
  }
  fclose(fp);
  // harp info
  snprintf(fn, sizeof(fn), "%s/%s-harp.%s", root, base, extn);
  if ((fp = fopen(fn, "w")) == NULL) return 1;
  for (h = 0; h < nHarp; h++) {
    fprintf(fp, "%d\n", harpInfo[h].id);
  }
  fclose(fp);
  // marp info
  snprintf(fn, sizeof(fn), "%s/%s-marp.%s", root, base, extn);
  if ((fp = fopen(fn, "w")) == NULL) return 1;
  for (m = 0; m < nMarp; m++) {
    fprintf(fp, "%d\n", marpInfo[m].id);
  }
  fclose(fp);
  // match-patch info
  snprintf(fn, sizeof(fn), "%s/%s-match.%s", root, base, extn);
  if ((fp = fopen(fn, "w")) == NULL) return 1;
  for (m = 0; m < nMarp; m++) {
    for (i = 0; i < nRec_all; i++) {
      inx = nRec_all * m + i; // offset into match-patch array (marp #m, time #i)
      if (print_if_invalid || matchInfo[inx].valid == Patch_Normal)
    match_ar_print(fp, "", matchInfo+inx);
    }
  }
  fclose(fp);
  // patch info
  snprintf(fn, sizeof(fn), "%s/%s-patch.%s", root, base, extn);
  if ((fp = fopen(fn, "w")) == NULL) return 1;
  for (h = 0; h < nHarp; h++) {
    for (i = 0; i < nRec_all; i++) {
      inx = nRec_all * h + i; // offset into patchInfo array (harp #h, time #i)
      if (print_if_invalid || patchInfo[inx].valid == Patch_Normal)
    patch_print(fp, "", tRec+i, patchInfo+inx);
    }
  }
  /*
   *  Synoptic map info
   */ 
  // ar rotation rate in deg/sec -- crude
  double degsec = 360.0/(27.2753*24*3600); // 360 degrees in one Carrington period
  double t0 = tRec[0].t; // time of origin
  double d_lon;
  const int Harp_smpl = 1;  // keep at 1 to always get singletons
  const int Marp_smpl = 10; // once per 2 hours

  // overall list, for harps and marps
  snprintf(fn, sizeof(fn), "%s/Patch-list.csv", root);
  if ((fp2 = fopen(fn, "w")) == NULL) return 1;
  // HARPs
  for (h = 0; h < nHarp; h++) {
    // don't export to synoptic map if too short
    // check: length=1  =>  rec1-rec0+1 = 1  => if() succeeds for nRec_min = 2
    if ((harpInfo[h].rec1 - harpInfo[h].rec0 + 1) < nRec_min)
      continue;
    snprintf(fn, sizeof(fn), "%s/harp-%06d.csv", root, harpInfo[h].id);
    if ((fp = fopen(fn, "w")) == NULL) return 1;
    for (inx1 = -1, i = 0; i < nRec_all; i += Harp_smpl) {
      inx = nRec_all * h + i; // offset into patchInfo array (harp #h, time #i)
      if (patchInfo[inx].valid == Patch_Normal) {
    inx1 = inx; // >= 0
    d_lon = (tRec[i].t - t0) * degsec; // longitude offset at current tRec
    fprintf(fp, "%.2f,%.2f,%.16s\n",  // YYYY.MM.DD_HH:MM, strip :00_TAI
        //  Mlat = matchInfo[Minx].lat;
        //  Mlon = matchInfo[Minx].lon;
        patchInfo[inx].stats[RS_ar_area_lat], 
        patchInfo[inx].stats[RS_ar_area_lon] - d_lon, 
        tRec[i].str);
      }
    }
    fclose(fp);
    // record HARP only if present
    if (inx1 >= 0) {
      double lat_wid = patchInfo[inx1].lat1 - patchInfo[inx1].lat0;
      double lon_wid = patchInfo[inx1].lon1 - patchInfo[inx1].lon0;
      fprintf(fp2, "%06d,%s,%.2f,%.2f\n", harpInfo[h].id, "harp", lat_wid, lon_wid);
    }
  }
  // MARPs
  for (h = 0; h < nMarp; h++) {
    snprintf(fn, sizeof(fn), "%s/marp-%06d.csv", root, marpInfo[h].id);
    if ((fp = fopen(fn, "w")) == NULL) return 1;
    for (inx1 = -1, i = 0; i < nRec_all; i += Marp_smpl) {
      inx = nRec_all * h + i; // offset into patchInfo array (marp #h, time #i)
      if (matchInfo[inx].valid == Patch_Normal) {
    inx1 = inx; // >= 0
    d_lon = (tRec[i].t - t0) * degsec; // longitude offset at current tRec
    fprintf(fp, "%.2f,%.2f,%.16s\n", // YYYY.MM.DD_HH:MM, strip :00_TAI
        matchInfo[inx].lat, 
        matchInfo[inx].lon - d_lon, 
        tRec[i].str);
      }
    }
    fclose(fp);
    // record MARP only if present
    if (inx1 >= 0) {
      double lat_wid = 0.0;
      double lon_wid = matchInfo[inx1].lon_wid;
      fprintf(fp2, "%06d,%s,%.2f,%.2f\n", marpInfo[h].id, "marp", lat_wid, lon_wid);
    }
  }
  fclose(fp2);

  // DONE
  V_printf(verbflag > 0, "", "Match dump: dump to `%s' complete.\n", root);
  return 0;
}

/*
 * match harps up with marps by computing distance matrix based on overlap
 *
 * this is obsolete.
 */
static
int
match_harp_marp_old(const char *rootDir,         // place for distance output, or NULL
            match_info_t *matchInfo, 
            marp_info_t  *marpInfo,
            int nMarp, 
            patch_info_t *patchInfo, 
            harp_info_t *harpInfo, 
            trec_info_t *tRec,
            int nRec_all, 
            int nRec_min, 
            int nHarp)
{
  int rec;
  int h, m;
  int Hinx, Minx;
  double d1, wt1, ct1;
  double diff, wt;
  double Hlat, Hlon, Mlat, Mlon;
  double *dist;

  // find the 2d overlap matrix
  if ((dist = calloc(nMarp*nHarp, sizeof(*dist))) == NULL) {
    V_printf(-1, "", "Could not allocate %d x %d match-dist matrix\n",
             nMarp, nHarp);
    return 1;
  }
  for (h = 0; h < nHarp; h++)
    for (m = 0; m < nMarp; m++) {
      d1 = 0;
      wt1 = 0;
      ct1 = 0;
      for (rec = 0; rec < nRec_all; rec++) {
    Hinx = h*nRec_all+rec;
    Minx = m*nRec_all+rec;
    if (patchInfo[Hinx].valid != Patch_Normal)
      continue;
    ct1++;
    if (matchInfo[Minx].valid != Patch_Normal)
      continue;
    // measure distance between the two patches
    Hlat = patchInfo[Hinx].stats[RS_ar_area_lat];
    Hlon = patchInfo[Hinx].stats[RS_ar_area_lon];
    Mlat = matchInfo[Minx].lat;
    Mlon = matchInfo[Minx].lon;
    diff = sqrt(1.00 * (Hlat - Mlat)*(Hlat - Mlat) + 
                0.25 * (Hlon - Mlon)*(Hlon - Mlon)); // weight lon less
    wt = patchInfo[Hinx].stats[RS_rgn_size]; // projected area
    d1  += diff * wt;
    wt1 += wt;
      } // end for rec
      // NB: no overlap => wt1 == 0 => dist = NaN
      dist[h*nMarp+m] = d1/wt1; // area-weighted distance
    } // end for h,m
  // export distance metric, etc. (rootDir == NULL means don't)
  if (match_harp_export(rootDir, dist, 
            matchInfo, marpInfo, nMarp, 
            patchInfo, harpInfo, tRec, 
            nRec_all, nRec_min, nHarp) != 0) {
    V_printf(-1, "", "Failed to export (fopen problem?), continuing.\n");
  }
  // assign ARs to HARPs to maximize this matrix
  double d_best;
  const double d_best_enough = 5.0;
  int h_best = 0;
  for (m = 0; m < nMarp; m++) {
    // match marp #m -- it is OK if there is no match
    for (d_best = d_best_enough*1e6, h = 0; h < nHarp; h++) {
      // note: NaNs do not satisfy this condition, which is correct
      if (d_best > dist[h*nMarp+m]) {
    d_best = dist[h*nMarp+m];
    h_best = h;
      }
    }
    if (d_best < d_best_enough) {
      // match harp h to marp m
      V_printf(verbflag > 1, "\t", "Matching HARP %d to MARP %d (%.2f)\n",
           harpInfo[h_best].id, marpInfo[m].id, d_best);
      // don't worry if two harps match the same marp
      harpInfo[h_best].match[harpInfo[h_best].n_match++] = m;
    } else {
      V_printf(verbflag > 1, "\t", "No match for MARP %d (%g)\n",
           marpInfo[m].id, d_best);
    }
  }
  free(dist); // forgot this in earlier versions
  return 0;
}

/*
 * match harps up with marps the new way: at each T_REC, find the 
 * harp-patch that the marp lands in. 
 *
 * The algorithm is roughly: 
 * 
 */
/* helper function: at one t_rec, find two harp-marp distances:
 *   dist1 = L2 distance of marp to patch box (0 if marp is within patch)
 *   dist2 = L2 distance of marp to patch centroid
 * If dist1 == 0, can use dist2 to further localize MARP w.r.t. HARP
 * Can have dist2 as a NaN if the HARP is a placeholder
 */
#define max(a,b) (((a)>(b))?(a):(b))

static
void
match_harp_marp_distance(match_info_t *matchI, 
             patch_info_t *patchI, 
             double *dist1, 
             double *dist2)
{
  double Mlat  = matchI->lat;
  double Mlon  = matchI->lon;
  double Hlat  = patchI->stats[RS_ar_area_lat]; // can be NaN
  double Hlon  = patchI->stats[RS_ar_area_lon];
  // *not* the patchI->stats[] box; it is tighter than the lat/lon box below
  double Hlat0 = patchI->lat0;
  double Hlon0 = patchI->lon0;
  double Hlat1 = patchI->lat1;
  double Hlon1 = patchI->lon1;

  // D{lat,lon}{0,1} is distance of MARP from the corresponding lat/lon line
  // e.g., Dlat0 is the distance of MARP from bottom of HARP box, or 0 if MARP 
  // is above bottom edge of HARP box.  If Dlat0 == Dlat1 == 0, MARP lands in 
  // the lat range of the HARP box.  At most one of Dlat0/Dlat1 can be nonzero, 
  // and similarly for Dlon0/Dlon1.
  double Dlat0 = max(0, Hlat0 - Mlat); // 0 if M above lat0
  double Dlon0 = max(0, Hlon0 - Mlon); // 0 if M right of lon0
  double Dlat1 = max(0, Mlat - Hlat1); // 0 if M below lat1
  double Dlon1 = max(0, Mlon - Hlon1); // 0 if M left of lon1
  double Dlat  = Dlat0 + Dlat1; // distance in latitude from MARP to HARP box
  double Dlon  = Dlon0 + Dlon1; // distance in lngitude from MARP to HARP box
  // weighted L2 distance of MARP from HARP bbox
  *dist1 = sqrt(1.00 * Dlat*Dlat + 0.25 * Dlon*Dlon);
  // weighted L2 distance of MARP from HARP centroid
  *dist2 = sqrt(1.00 * (Hlat - Mlat)*(Hlat - Mlat) + 
        0.25 * (Hlon - Mlon)*(Hlon - Mlon)); // weight lon less
}

// tidy namespace
#undef max

/*
 * comparison utilities for qsort() calls in match_harp_marp
 */
// decorate-sort-undecorate idiom: sort by size, descending
static int 
match_sort_cmp_size(const void *av, const void *bv) {
  const marp_sort_t *a = (const marp_sort_t *) av;
  const marp_sort_t *b = (const marp_sort_t *) bv;
  return b->size - a->size; // flipped sense for descending
}
// decorate-sort-undecorate idiom: sort by MARP id, ascending
static int 
match_sort_cmp_id(const void *av, const void *bv) {
  const marp_sort_t *a = (const marp_sort_t *) av;
  const marp_sort_t *b = (const marp_sort_t *) bv;
  return a->id - b->id; // regular sense
}

static
int
match_harp_marp(const char *rootDir,     // place for distance output, or NULL
        match_info_t *matchInfo, 
        marp_info_t  *marpInfo,
        int nMarp, 
        patch_info_t *patchInfo, 
        harp_info_t *harpInfo, 
        trec_info_t *tRec,
        int nRec_all, 
        int nRec_min,           // unused
        int nHarp)
{
  int rec;
  int h, m, h_best, n_match;
  int Hinx, Minx;
  double d_box, d_box_best;
  double d_ctr, d_ctr_best;
  double ct_best;
  // OLD:
  // const double d_best_enough = 5.0; // degrees
  // more strict:
  const double d_best_enough = 0.5; // degrees
  double *HMct, *Msize;
  marp_sort_t match1, *match_list; // hold MARPs matching a given HARP

  if ((HMct = calloc(nMarp*nHarp, sizeof(*HMct))) == NULL) {
    V_printf(-1, "", "Could not allocate %d x %d match-count matrix\n",
             nMarp, nHarp);
    return 1;
  }
  if ((Msize = calloc(nMarp, sizeof(*Msize))) == NULL) {
    V_printf(-1, "", "Could not allocate %d-length match-size vector\n",
             nMarp);
    return 1;
  }
  if ((match_list = calloc(nMarp, sizeof(*match_list))) == NULL) {
    V_printf(-1, "", "Could not allocate %d-length match-sort vector\n",
             nMarp);
    return 1;
  }
  // fill in HMct and Msize
  for (rec = 0; rec < nRec_all; rec++) {
    for (m = 0; m < nMarp; m++) {
      // does marp m exist at this time?
      Minx = m*nRec_all+rec;
      if (matchInfo[Minx].valid != Patch_Normal) continue;
      Msize[m] += (matchInfo[Minx].area > 0) ? matchInfo[Minx].area : 1; // can be <0
      // find the best HARP match to the MARP
      //   e.g.: MARP lies in HARP bbox, and bonus: MARP/HARP centers close
      h_best = -1; // flag for "no match" (e.g., no HARPs present)
      for (h = 0; h < nHarp; h++) {
    // does harp h exist at this time?
    Hinx = h*nRec_all+rec;
    if (patchInfo[Hinx].valid != Patch_Normal) continue;
    // d_box: is MARP within HARP?  d_ctr: how close are centers?
    match_harp_marp_distance(matchInfo+Minx, patchInfo+Hinx, &d_box, &d_ctr);
    // find best match
    if ((h_best < 0) || (d_box < d_box_best)) {
      d_box_best = d_box;
      d_ctr_best = d_ctr;
      h_best     = h;
    } else if (( d_box == d_box_best) && 
           ((d_ctr  < d_ctr_best) || isnan(d_ctr_best))) {
      d_ctr_best = d_ctr;
      h_best     = h;
    }
      } // end for h
      // was there a match, and was it good enough
      if ((h_best >= 0) &&
      (d_box_best < d_best_enough)) {
    HMct[h_best*nMarp + m]++; // one more matching T_REC
      }
    }
  }

  // clean up poor matches, which are probably erroneous
  // (note, we do this for a fixed MARP, across all HARPs, because
  // we expect a MARP to match one HARP.  but, we don't do it for a 
  // fixed HARP, across MARPs, because we expect HARPs to sometimes
  // match many MARPs)
  for (m = 0; m < nMarp; m++) {
    // find the best HARP match to the MARP
    // (note: ending with ct_best == 0 is OK)
    h_best = -1;
    for (ct_best = 0, h = 0; h < nHarp; h++) {
      if (HMct[h*nMarp+m] > ct_best) {
    ct_best = HMct[h*nMarp+m]; // find max ct[]
    h_best = h;
      }
    }
    marpInfo[m].top_match = h_best; // just keep the best one, not the whole list
    // zero-out matches worse than ct_best * 0.25 
    //  (e.g., matched on 80 T_RECs vs. 20 T_RECs)
    for (h = 0; h < nHarp; h++)
      if (HMct[h*nMarp+m] < ct_best*0.25)
    HMct[h*nMarp+m] = 0;
  }
  
  // (optionally write out HMct)
  // export distance metric, etc. (rootDir == NULL means don't)
  if (match_harp_export(rootDir, HMct, 
            matchInfo, marpInfo, nMarp, 
            patchInfo, harpInfo, tRec, 
            nRec_all, nRec_min, nHarp) != 0) {
    V_printf(-1, "", "Failed to export (fopen problem?), continuing.\n");
  }

  // for each HARP, find all matching MARPs
  //  (this loop order reverses the above MARP-to-HARP association)
  for (h = 0; h < nHarp; h++) {
    for (n_match = m = 0; m < nMarp; m++) {
      if (HMct[h*nMarp+m] > 0) {
    // summary of the match
    match1.marp_index = m; // 0..nMarp-1
    match1.id         = marpInfo[m].id; // AR number
    match1.count      = HMct[h*nMarp+m];
    match1.size       = Msize[m];
    // insert MARP m in match_list
    match_list[n_match++] = match1;
    // add this match to marp summary info
    marpInfo[m].n_match++;
      }
    }
    // sort matches by size, to find the largest ones
    qsort(match_list, n_match, sizeof(marp_sort_t), match_sort_cmp_size);
    if (n_match > 0)
      harpInfo[h].top_match = match_list[0].marp_index; // 0..nMarp-1
    else
      harpInfo[h].top_match = -1; // illegal value -- no match
    // cap on the number of good matches, according to size
    if (n_match > HARP_INFO_MATCH_MAX)
      n_match = HARP_INFO_MATCH_MAX;
    // re-sort top matches according to MARP id
    qsort(match_list, n_match, sizeof(marp_sort_t), match_sort_cmp_id);
    // insert these matches into HARP match-list
    harpInfo[h].n_match = n_match;
    for (m = 0; m < n_match; m++) {
      // match harp h to marp m
      V_printf(verbflag > 1, "\t", "Matching HARP %d to MARP %d (ct=%.0f,siz=%f)\n",
           harpInfo[h].id, 
           match_list[m].id,
           match_list[m].count, 
           match_list[m].size);
      // don't worry if two harps match the same marp
      harpInfo[h].match[m] = match_list[m].marp_index; // marp_index in 0..nMarp-1
    }
  }
  // summary info for each MARP
  char matchStr[20];
  for (m = 0; m < nMarp; m++) {
    // the best match to the Marp, as a string
    snprintf(matchStr, sizeof(matchStr), "[best: %d]", 
         (marpInfo[m].top_match >= 0) ? harpInfo[marpInfo[m].top_match].id : -1);
    // usually 0 or 1 matches, so just print best if it exists
    V_printf(verbflag > 0, "\t", "Matches found for MARP %d: %d %s\n", 
         marpInfo[m].id, marpInfo[m].n_match,
         (marpInfo[m].n_match > 0) ? matchStr : "");
  }
  free(HMct);
  free(Msize);
  free(match_list);
  return 0;
}

/***********************************************************************
 *
 * Module input arguments
 *
 ***********************************************************************/

char *module_name = "ingest_mharp";
char *hmi_mharp_version = "1.1";  // used for keywords

ModuleArgs_t module_args[] = {
  {ARG_STRING, "root",  kNOT_SPEC,   "Input file directory"},
  {ARG_STRING, "lists", "track-new.txt",   "File(s) listing HARP numbers (comma-separated)"},
  {ARG_INT,    "match", "0",         "Dump match-AR summary info? (<0 for post-dump exit)"},
  {ARG_STRING, "out",   kNOT_SPEC,   "Output HARP data series"},
  {ARG_STRING, "mag",   "hmi.M_720s","Input magnetogram data series"},
  {ARG_STRING, "mask",  kNOT_SPEC,   "Input mask data series"},
  {ARG_STRING, "href",  kNOT_SPEC,   "Reference HARP bounding box source (default: none)"},
  {ARG_STRING, "t0cut", kNOT_SPEC,   "Earliest time allowed to ingest (default: none)"},
  {ARG_STRING, "t1cut", kNOT_SPEC,   "Latest time allowed to ingest (default: none)"},
  {ARG_INT,    "trec",  "0",         "Number of T_RECs to ingest (0=all, 1=latest)"},
  {ARG_INT,    "tpad",  "120",       "Temporal padding (in images, >=0, 120=1 day)"},
  {ARG_INT,    "tmin",  "1",         "Minimum T_REC's needed to ingest a patch (1=all)"},
  {ARG_DOUBLE, "cadence","720.0",    "Image cadence in seconds"},
  {ARG_INT,    "readonly", "0",      "Read-only with respect to rootdir (if ro=1, match=0)"},
  {ARG_INT,    "verb",  "1",         "Verbosity: 0=errs only; 1, 2 = more"},
  {ARG_END}
};

/*
 *    
 * The nested loops in the ingestion step below follow this form:
 *
 *   for t in T_RECs:
 *     if magnetogram(t) is unavailable:
 *       continue // all HARPs at t fail -- not an error
 *     if mask(t) is unavailable:
 *       continue // all HARPs at t fail -- is an error
 *     for h in harps(t):
 *       if ingest(patches(h, t)):
 *         patch.success = 1;
 */

int DoIt(void)
{
  // input parameters
  const char *listNamesIn;
  const char *rootDir;
  int matchStatus;    // 0 -> no match output; +1/-1 -> output; -1 for exit after
  int readOnly;       // 0 -> regular, 1 -> no update of files in "root"
  char *outQuery, *maskQuery, *magQuery;
  const char *refHarp;    // reference HARP bboxes
  char *t0cutIn, *t1cutIn;
  char fn_log[STR_MAX], fn_err[STR_MAX];
  // records, segments, arrays
  DRMS_RecordSet_t *inRecSet, *magRecSet, *maskRecSet;
  DRMS_Record_t *inRec, *maskRec, *magRec;
  DRMS_Array_t *maskData;
  char *maskImg;
  int maskStride;
  // harp-list files
  int listNum;
  char *listNames[LIST_MAX]; // heads of all list file names, NULL terminated
  char *listBuf;             // buffer pointed in to by above
  // HARP bookkeeping arrays
  patch_info_t *patchInfo = NULL;  // nHarp X nRec_all array of tracks
  harp_info_t *harpInfo = NULL;    // length = nHarp
  trec_info_t *tRec = NULL;        // length = nRec_all
  run_info_t runInfo;              // overall run info (tracker params)
  // matching-region bookkeeping
  int nMarp;                       // number of match-regions (NOAA ARs)
  match_info_t *matchInfo = NULL;  // patch info of above, like patchInfo
  marp_info_t *marpInfo = NULL;    // per-match-region info (like harpInfo)
  // HARP times, counts, etc.
  TIME t0, t1;               // times bracketing HARPs listed to ingest
  char t0_S[32], t1_S[32];   // time strings for above
  TIME t0cut, t1cut;         // limit times for ingestion (eg, subset of [t0,t1])
  char t0cut_S[32], t1cut_S[32]; // time strings for above
  int rec0cut, rec1cut;      // T_REC indexes for above
  int nHarp;                 // number of HARPs being ingested
  int nRec;                  // number of T_REC slots in union of *seen* HARPs
  int nRec_pad;              // padding (context) on either *unseen* end of HARP
  int nRec_all;              // total number of T_REC slots = nRec + 2*nRec_pad
  int nRec_eat;              // number of T_REC slots to ingest (eat)
  int nRec_min;              // number of T_REC slots required before ingestion
  double cadence;            // image cadence in seconds
  int npix_fulldisk;         // (square) fulldisk image size
  int Harp_try;              // Valid HARPs at current T_REC
  int Harp_ok, Harp_err;     // ingested OK vs. error
  int rec;                   // record counter
  int i;                     // harp counter
  int rv, inx;               // misc ints
  char *err;
  // Time measuring -- ugly, but functional
  double wt0, wt;
  double ut0, ut;
  double st0, st;
  double ct0, ct;
    
  // timings
  wt0 = getwalltime();
  ct0 = getcputime(&ut0, &st0);
  // Key command line parameter -- root directory
  rootDir = cmdparams_get_str(&cmdparams, "root", NULL);
  if (strcmp(rootDir, kNOT_SPEC) == 0) DIE("Root directory must be specified");
  // know root dir, so can open log files
  snprintf(fn_log, sizeof(fn_log), "%s/%s", rootDir, FN_INGEST_STDLOG);
  snprintf(fn_err, sizeof(fn_err), "%s/%s", rootDir, FN_INGEST_STDERR);
  if ((LOGout = fopen(fn_log, "w")) == NULL)
    DIE("Could not open regular log file `%s' for writing", fn_log);
  if ((LOGerr = fopen(fn_err, "w")) == NULL)
    DIE("Could not open error log file `%s' for writing", fn_err);
  // Remaining command line parameters
  outQuery = (char *) cmdparams_get_str(&cmdparams, "out", NULL);
  if (strcmp(outQuery, kNOT_SPEC) == 0) DIE("Output HARP series must be specified");
  listNamesIn = cmdparams_get_str(&cmdparams, "lists", NULL);
  if (strcmp(listNamesIn, kNOT_SPEC) == 0) DIE("HARP list file(s) must be specified");
  matchStatus  = cmdparams_get_int(&cmdparams, "match", NULL);
  readOnly     = cmdparams_get_int(&cmdparams, "readonly", NULL);
  if (readOnly) matchStatus = 0; // do not write match info
  magQuery  = (char *) cmdparams_get_str(&cmdparams, "mag",  NULL);
  maskQuery = (char *) cmdparams_get_str(&cmdparams, "mask", NULL);
  if (strcmp(maskQuery, kNOT_SPEC) == 0) DIE("mask series must be specified");
  refHarp  = cmdparams_get_str(&cmdparams, "href", NULL);
  if (strcmp(refHarp, kNOT_SPEC) == 0) refHarp = NULL; // no reference HARPs: OK
  nRec_eat = cmdparams_get_int(&cmdparams, "trec", NULL);
  nRec_pad = cmdparams_get_int(&cmdparams, "tpad", NULL);
  nRec_min = cmdparams_get_int(&cmdparams, "tmin", NULL);
  cadence  = cmdparams_get_double(&cmdparams, "cadence", NULL);
  // cut times: strings -> numeric times -> normalized strings
  t0cutIn = (char *) cmdparams_get_str(&cmdparams, "t0cut", NULL);
  t1cutIn = (char *) cmdparams_get_str(&cmdparams, "t1cut", NULL);
  verbflag = cmdparams_get_int(&cmdparams, "verb", NULL); // an int, not a flag
  // check pad/ingest counts for sanity
  if (nRec_min <= 0)
    DIE("tmin must be a positive integer");
  if (nRec_eat < 0)
    DIE("trec must be a nonnegative integer");
  if (nRec_pad < 0)
    DIE("tpad must be a nonnegative integer");
  if ((nRec_eat > 0) && (nRec_eat <= nRec_pad))
    DIE("number of trec to ingest (trec) must be larger than padding (tpad)");
  // check and apply logic to cut times
  if (strcmp(t0cutIn, kNOT_SPEC) == 0) {
    t0cut = NAN;
    strcpy(t0cut_S, "no cut-time before");
  } else {
    t0cut = sscan_time(t0cutIn);
    if (time_is_invalid(t0cut)) DIE("Could not parse given t0cut");
    sprint_time(t0cut_S, t0cut, "TAI", 0);
  }
  if (strcmp(t1cutIn, kNOT_SPEC) == 0) {
    t1cut = NAN;
    strcpy(t1cut_S, "no cut-time after");
  } else {
    t1cut = sscan_time(t1cutIn);
    if (time_is_invalid(t1cut)) DIE("Could not parse given t1cut");
    sprint_time(t1cut_S, t1cut, "TAI", 0);
  }
  if (t0cut > t1cut) DIE("Need t0cut <= t1cut");
  if (!isnan(t0cut) || !isnan(t1cut))
    V_printf(1, "", "Chosen T_REC cut interval: %s - %s\n", t0cut_S, t1cut_S);
  // separate list files into pieces
  listNum = 0;
  listBuf = strdup(listNamesIn); // below, we write \0's into list names, so must dup
  bzero(listNames, sizeof(listNames)); // ensure it's all zero
  for (char **listp = listNames; (*listp = strsep(&listBuf, ",")) != NULL; listNum++)
    if (**listp != '\0')
      if (++listp >= &listNames[LIST_MAX])
    break; // no more room
  if (listNum >= LIST_MAX || listNames[listNum] != NULL)
    DIE("Got too many list files, reduce or recompile");
  // summary info for input args
  V_printf(verbflag > 0, "", "Ingesting from `%s'\n", rootDir);
  V_printf(verbflag > 0, "", "  Got %d list file%s from `%s'\n", 
       listNum, listNum == 1 ? "" : "s", listNamesIn);
  V_printf(verbflag > 0, "", "Ingesting into `%s'\n", outQuery);
  V_printf(verbflag > 0, "", "Mask input series `%s'\n", maskQuery);
  if (refHarp)
    V_printf(verbflag > 0, "", "Reference harp series is `%s'\n", refHarp);

  /*
   * First pass: scan all tracks to get dimensions of patch info
   */
  if (load_run_info(rootDir, listNames, cadence, &nHarp, &nRec, &t0, &t1))
    DIE("Error finding basic run info (bad rootDir, listfile, or format?)");
  if (load_tracker_params(rootDir, &runInfo))
    WARN("Failed to load tracker parameter file.  Some keys will not be set");
  if (nRec == 0) {
    V_printf(1, "", "Got nHarp = 0.  Using nominal time range.  Setting pad = 0.\n");
    nRec_pad = 0;
  }
  nRec_all = nRec + 2*nRec_pad; // number of T_REC slots needed to cover the HARPs
  // progress update
  sprint_time(t0_S, t0, "TAI", 0);
  sprint_time(t1_S, t1, "TAI", 0);
  V_printf(verbflag > 0, "", "nHarp = %d; nRec = %d; nRec_all = %d; pad = %d.\n", 
       nHarp, nRec, nRec_all, nRec_pad);
  V_printf(verbflag > 0, "", "t0 = %s\n", t0_S); 
  V_printf(verbflag > 0, "", "t1 = %s\n", t1_S); 
  // set up "ingest-most-recent" record numbers (NRT mode)
  if (nRec_eat == 0) 
    nRec_eat = nRec_all; // fill in for the shorthand 0
  if (nRec_eat > nRec_all) {
    nRec_eat = nRec_all;  // cap the number to be ingested
    // this is nonfatal, but warn anyway
    WARN("Specified trec to ingest was larger than total trecs + padding; capping it");
  }
  // set up "cut" record numbers (gap-filling mode)
  if (t0cut < t0 || t0cut > t1) WARN("Specified t0cut outside [t0,t1].");
  if (t1cut < t0 || t1cut > t1) WARN("Specified t1cut outside [t0,t1].");
  rec0cut = time_to_recnum(t0cut, t0, nRec_pad, -1, cadence); // t0cut is NaN -> rec0cut = -1
  rec1cut = time_to_recnum(t1cut, t0, nRec_pad, -1, cadence); // t1cut is NaN -> rec1cut = -1
  // cap rec0cut; the complementary case (rec1 < 0) will look like NaN, not worth fixing.
  if (rec0cut > nRec_all) rec0cut = nRec_all; 
  // summary -- might use help?
  if (rec0cut >= 0 || rec1cut >= 0) 
    // gap-fill
    V_printf(verbflag > 0, "", "Cut times: rec0cut = %d, rec1cut = %d\n", rec0cut, rec1cut);
  else 
    // usual case
    V_printf(verbflag > 0, "", 
         "Compiling %d HARPs over %d distinct T_REC's among %d total T_REC's.\n", 
         nHarp, nRec_eat, nRec_all); 
  
  /*
   * Second pass: Insert all HARP information into the master arrays
   */
  // Allocate the informational arrays
  patchInfo = calloc(nHarp * nRec_all, sizeof(*patchInfo));
  harpInfo  = calloc(nHarp,            sizeof(*harpInfo));
  tRec      = calloc(        nRec_all, sizeof(*tRec));
  if (!tRec || !harpInfo || !patchInfo)
    DIE("Memory allocation failure");
  // initialize T_REC array -- basic time info
  V_printf(verbflag > 0, "", "Setting T_REC times\n");
  set_trec_times(tRec, t0, nRec_all, nRec_pad, cadence);
  // insert image availability + get WCS and image dimension
  V_printf(verbflag > 0, "", "Checking mag/mask availability and WCS across T_REC\n");
  set_trec_hasdata(tRec, magQuery, maskQuery, &npix_fulldisk, nRec_all);
  // load HARP ID array from file
  V_printf(verbflag > 0, "", "Getting HARP id numbers\n");
  if (load_harp_ids(rootDir, listNames, harpInfo) != 0)
    DIE("Error creating HARP id list (bad listfile or format?)");
  // fill in the HARP x T_REC master index
  V_printf(verbflag > 0, "", "Getting detailed HARP info\n");
  if (load_all_patch_info(rootDir, refHarp,
              patchInfo, harpInfo, 
              tRec, nHarp, nRec_all, nRec_pad, cadence, npix_fulldisk))
    DIE("Error fetching HARP info (bad listfiles/format)");
  // mark patches eligible to ingest; count eligible patches at each T_REC,HARP
  mark_eligible_patches(patchInfo, harpInfo, tRec, nHarp, 
            nRec_all, nRec_eat, rec0cut, rec1cut, nRec_min);
  // elapsed time so far
  wt = getwalltime();
  ct = getcputime(&ut, &st);
  V_printf(1, "", "elapsed time after indexing: %.3f s wall, %.3f s cpu\n",
       (wt - wt0)*1e-3, (ct - ct0)*1e-3);

  /*
   * Third step: Compute matching ARs
   */
  V_printf(verbflag > 0, "", "Loading match regions\n");
  // load the matching ARs
  if (match_load_resample(nRec_all, tRec, &nMarp, &marpInfo, &matchInfo) != 0) {
    V_printf(-1, "", "Could not load and resample Match-AR metadata.\n");
    DIE("Unable to read MATCH AR information");
  }
  if (nMarp == 0)
    V_printf(1, "", "No match ARs supplied/found, continuing.\n");
  V_printf(verbflag > 0, "", "Matching %d HARPs to %d Match-regions\n", nHarp, nMarp);
  // feb 2013: this is the new matcher
  match_harp_marp((matchStatus == 0) ? NULL : rootDir, // NULL -> no output
          matchInfo, marpInfo, nMarp, 
          patchInfo, harpInfo, tRec, 
          nRec_all, nRec_min, nHarp);
  free(matchInfo);  // no longer needed, but marpInfo is
  matchInfo = NULL;
  // elapsed time so far
  wt = getwalltime();
  ct = getcputime(&ut, &st);
  V_printf(1, "", "elapsed time after matching: %.3f s wall, %.3f s cpu\n",
       (wt - wt0)*1e-3, (ct - ct0)*1e-3);
  // allow exit without ingestion to facilitate match testing
  if (matchStatus < 0) {
    V_printf(-1, "", "Early exit due to given match option (match = %d)\n", matchStatus);
    // close log files
    if (LOGout) fclose(LOGout);
    if (LOGerr) fclose(LOGerr);
    return DRMS_SUCCESS;
  }

  /* 
   * Fourth step: Ingest into JSOC, with outer loop done across T_REC
   *
   *   (this should be factored out as a separate routine)
   */
  // T_REC loop covers all times
    // This double loop (over records, then HARPs) allocates memory for images (patchInfo[i].image).
    // These images are not used outside this loop.
  for (rec = 0; rec < nRec_all; rec++) {
    if (tRec[rec].nharp == 0) {
      V_printf(verbflag > 0, "\t", "[%3d/%3d] %s: Nothing to ingest.\n", 
           rec+1, nRec_all, tRec[rec].str);
      continue; 
    }
    V_printf(verbflag > 1, "\t", "Ingesting at %s [%d/%d]\n", tRec[rec].str, rec+1, nRec_all);
    // Open mag record -- this should succeed
    V_printf(verbflag > 2, "\t\t", "Opening mag: %s[%s]\n", magQuery, tRec[rec].str);
    err = load_drms_rec_seg(&magRecSet, &magRec, NULL, NULL,
                magQuery, tRec[rec].str, NULL, DRMS_TYPE_RAW);
    if (err) {
      V_printf(-1, "", "%s (%s[%s]) -- skipping this T_REC.\n", err, magQuery, tRec[rec].str);
      continue;
    }
    // Open mask record + data segment
    V_printf(verbflag > 2, "\t\t", "Opening mask: %s[%s]\n", maskQuery, tRec[rec].str);
    err = load_drms_rec_seg(&maskRecSet, &maskRec, &maskData, NULL, 
                maskQuery, tRec[rec].str, "mask", DRMS_TYPE_CHAR);
    if (err) {
      V_printf(-1, "", "%s (%s[%s]) -- skipping this T_REC.\n", err, maskQuery, tRec[rec].str);
      drms_close_records(magRecSet, DRMS_FREE_RECORD); magRecSet  = NULL;
      continue;
    }
    maskImg    = (char *) maskData->data;
    maskStride = maskData->axis[0];

    Harp_try = Harp_ok = Harp_err = 0;
    for (i = 0; i < nHarp; i++) {
      inx = i*nRec_all + rec;
      if (patchInfo[inx].valid == Patch_Invalid || patchInfo[inx].ingest != Patch_Ingest_Eligible)
    continue;
      V_printf(verbflag > 2, "\t\t", "  HARP[%d][%s] begin\n", harpInfo[i].id, tRec[rec].str);
      // ingest the HARP at the current time
      Harp_try++;
      rv = ingest_harp(patchInfo+inx, harpInfo+i, tRec+rec, 
               tRec,             // all time info (index->time map for keys)
               &runInfo,         // run info, also for keys
               marpInfo,         // matching ARs (index->id map for keys)
               magRec, maskRec,  // for keys and data links
               maskImg, maskStride, // mask and stride (NAXIS1)
               outQuery);
      if (rv != 0) {
    V_printf(-1, "", "Error ingesting HARP[%d][%s], skipped\n", harpInfo[i].id, tRec[rec].str);
    Harp_err++;
    continue;
      }
      Harp_ok++;
      patchInfo[inx].success = 1; // it is left at zero in all other cases
      V_printf(verbflag > 1, "\t", "  HARP[%d][%s] done\n", harpInfo[i].id, tRec[rec].str);
    }
    V_printf(verbflag > 0, "\t", "[%3d/%3d] %s:\t%2d attempted = %2d OK + %2d errors%s\n", 
         rec+1, nRec_all, tRec[rec].str, Harp_try, Harp_ok, Harp_err,
         Harp_err > 0 ? " ***" : "");
    drms_free_array(maskData);
    drms_close_records(magRecSet,  DRMS_FREE_RECORD); magRecSet  = NULL;
    drms_close_records(maskRecSet, DRMS_FREE_RECORD); maskRecSet = NULL;
  } // end double loop over records and HARPs.
  /* 
   * Fifth step: Record ingested HARPs
   */
  // loop over HARPs, marking them as ingested (suppress marking if readOnly)
  for (i = 0; i < nHarp; i++) {
    if (mark_ingested_harp(readOnly ? NULL : rootDir, &runInfo, harpInfo+i) != 0)
      V_printf(-1, "", "HARP %d was ingested into JSOC but not noted in status file.\n", 
           harpInfo[i].id);
  }

  /* 
   * Find and print summary info
   */
  trec_ingest_tab_t tRec_tab;
  patch_ingest_tab_t patch_tab;
  int clean_exit = 1;

  tabulate_patches(&patch_tab, &tRec_tab, patchInfo, harpInfo, tRec, nHarp, nRec_all);
  
  // summary across HARPs
  rv = 0;
  V_printf(1, "", "Ingestion summary -- across HARPs:\n");
  for (i = 0; i < nHarp; i++) {
    // this happens for CRs, and if a time cutoff is given
    if (harpInfo[i].n_patch_eligible == 0) continue;
    rv = 1;
    // per-HARP summary
    V_printf(1, "", "  HARP[%d]: %4d T_RECs eligible = %4d ok + %d errs.\n", 
         harpInfo[i].id, harpInfo[i].n_patch_eligible, 
         harpInfo[i].n_patch_ingested, harpInfo[i].n_patch_err);
    // rare case where an ingest was attempted, and failed
    if (harpInfo[i].n_patch_err > 0)
      V_printf(-1, "", "WARNING: HARP[%d] had %d ingestion errors in %d eligible T_RECs\n",
           harpInfo[i].id, harpInfo[i].n_patch_err, harpInfo[i].n_patch_eligible);
  }
  if (rv == 0)
    V_printf(1, "", "  (No HARPs ingested)\n");

  // summary across T_RECs
  V_printf(1, "", "Ingestion summary -- across T_RECs:\n");
  V_printf(1, "", "  Considered NT = %d T_RECs.\n", nRec_all);
  V_printf(1, "", "  Ingested across %d T_RECs (< NT is OK).\n", tRec_tab.trec_patch_valid);
  V_printf(1, "", "  Ingest error in %d T_RECs.\n", tRec_tab.trec_patch_err);
  for (i = 0; i < Trec_Data_NUM; i++) {
    V_printf(1, "", "    Of NT, %d T_RECs were: %s\n", tRec_tab.hist[i], Trec_Data_Names[i]);
  }
  if (tRec_tab.hist[Trec_Data_Missing] > 0) {
    clean_exit = 0;
    for (i = -1; i <= 1; i += 2)
      V_printf(i, "", "WARNING: %d mask records were missing (not due to simple mag gap).\n",
           tRec_tab.hist[Trec_Data_Missing]);
  }

  // summary across patches
  V_printf(1, "", "Ingestion summary -- across all patches (HARP appearances):\n");
  V_printf(1, "", "  Processed NP = %d patches over all HARPs and all T_RECs.\n",
       nRec_all*nHarp - patch_tab.hist[Patch_Ingest_Invalid]);
  for (i = 0; i < Patch_Ingest_NUM; i++) {
    if (i == (int) Patch_Ingest_Invalid) continue; // skip these
    V_printf(1, "", "    Of NP, %d patches were: %s\n", patch_tab.hist[i], Patch_Ingest_Names[i]);
  }
  V_printf(1, "", "  Attempted to ingest NP1 = %d patches over all T_RECs.\n",
       patch_tab.hist[Patch_Ingest_Eligible]);
  V_printf(1, "", "    Of NP1, %d patches were ingested successfully.\n", patch_tab.try_ok);
  V_printf(1, "", "    Of NP1, %d patches had ingestion errors.\n", patch_tab.try_err);
  if (patch_tab.try_err > 0) {
    clean_exit = 0;
    for (i = -1; i <= 1; i += 2)
      V_printf(i, "", "WARNING: %d patch ingestion errors were NOT due to missing mags.\n",
           patch_tab.try_err);
  }

  // NULLs OK here
  free(marpInfo);
  free(listBuf);
    
  // Free alloc'd strings.
  for (i = 0; i < nHarp * nRec_all; i++)
  {
      if (patchInfo[i].patchName)
      {
          free(patchInfo[i].patchName);
          patchInfo[i].patchName = NULL;
      }
  }
    
  free(patchInfo);
  free(tRec);
  free(harpInfo);

  // elapsed time
  wt = getwalltime();
  ct = getcputime(&ut, &st);
  V_printf(1, "", "total time used: %.3f s wall, %.3f s cpu\n",
       (wt - wt0)*1e-3, (ct - ct0)*1e-3);
  if (clean_exit) {
    V_printf(1, "", "Exiting OK.\n");
  } else {
    for (i = -1; i <= 1; i += 2)
      V_printf(i, "", "Exiting with some warnings.\n");
  }
  // close log files
  if (LOGout) fclose(LOGout);
  if (LOGerr) fclose(LOGerr);
  return DRMS_SUCCESS;
}

// end of main routine

/*
 * tabulate statistics of ingestion success for each patch and t-rec
 *   returns nonzero if bad code (internal error)
 */
static
int
tabulate_patches(patch_ingest_tab_t *pTab,
         trec_ingest_tab_t  *tTab,
         patch_info_t *patchInfo, 
         harp_info_t  *harpInfo, 
         trec_info_t  *tRec, 
         int nHarp,
         int nRec_all)
         
{
  int rec, i, inx;           // counters
  int something_failed;
  patch_ingest_t p_code;
  trec_data_t d_code;

  // zero out accumulators -- cannot assume they were calloc'd
  pTab->try_ok  = 0;
  pTab->try_err = 0;
  bzero(pTab->hist, sizeof(pTab->hist));
  tTab->trec_patch_err   = 0;
  tTab->trec_patch_valid = 0;
  bzero(tTab->hist, sizeof(tTab->hist));
  for (i = 0; i < nHarp; i++)
    harpInfo[i].n_patch_eligible   = 
      harpInfo[i].n_patch_ingested =
      harpInfo[i].n_patch_err      = 0;
  // accumulate
  for (rec = 0; rec < nRec_all; rec++) {
    /*
     * histogram of trec[rec].data
     */
    d_code = tRec[rec].data;
    if (d_code < 0 || d_code >= Trec_Data_NUM) {
      V_printf(-1, "", 
           "Internal error: Bad trec code %d at %s in tabulate_patches.\n", 
           d_code, tRec[rec].str);
      return 1;
    }
    tTab->hist[(int) d_code]++;
    /*
     * inspect all patches at this T_REC
     */
    something_failed = 0;
    for (i = 0; i < nHarp; i++) {
      inx = i*nRec_all + rec;
      p_code = patchInfo[inx].ingest;
      if (p_code < 0 || p_code >= Patch_Ingest_NUM) {
    V_printf(-1, "", 
         "Internal error: Bad patch code %d at %s, harp %d, in tabulate_patches.\n", 
         p_code, tRec[rec].str, patchInfo[inx].num);
    return 1;
      }
      // make a histogram
      pTab->hist[(int) p_code]++;
      // for the ones we tried to ingest -- tabulate success and failure
      if (p_code == Patch_Ingest_Eligible) {
    harpInfo[i].n_patch_eligible++;
    if (patchInfo[inx].success) {
      harpInfo[i].n_patch_ingested++;
      pTab->try_ok++;
    } else {
      harpInfo[i].n_patch_err++;
      pTab->try_err++;
      something_failed = 1;
    }
      }
    } // end for i
    /*
     * per-T_REC counts
     */
    if (something_failed)
      tTab->trec_patch_err++;  // some patch failed at T_REC
    if (tRec[rec].nharp > 0)
      tTab->trec_patch_valid++; // there was something to ingest at T_REC
  } // end for rec
  return 0;
}



/***********************************************************************
 *
 * Utility routines
 *
 ***********************************************************************/

/*
 * Convert time (drms "TIME") to record number (integer)
 *   - rounds the record number off
 *   - optionally clips it to 0..nRec-1 if nRec >= 0
 *   - the key is that at t0, the record number is rec0
 *   - NaN time returns (-1) record number
 */
static
int
time_to_recnum(TIME t, TIME t0, int rec0, int nRec, double cadence)
{
  if (isnan(t)) return -1;
  double rec_diff = (t - t0) / cadence;
  int rec = rec0 + ((int) rint(rec_diff));
  // optionally force in-range rec
  if (nRec >= 0) {
    if (rec < 0) rec = 0;
    if (rec >= nRec) rec = nRec-1;
  }
  return rec;
}

/*
 * get ephemeris into *wcs from a DRMS record
 *   returns 0 for success, nonzero for failed getkey or NaN in result
 *   note that all fields in *wcs are named exactly as the corresponding
 *   keyword, except in lower-case
 */

static
int 
wcs_load(wcs_t *wcs, DRMS_Record_t *rec)
{
  int status, retval;
  const size_t Nmiss = sizeof(wcs->missing);

  // RSUN_REF and DSUN_OBS are doubles in the JSD, rest are floats
  status = 0;
  wcs->rsun_ref = drms_getkey_double(rec, "RSUN_REF", &status);
  if (status) wcs->rsun_ref = 6.96e8;
  status = 0;
  wcs->dsun_obs = drms_getkey_double(rec, "DSUN_OBS", &status);
  wcs->cdelt1   = drms_getkey_float (rec, "CDELT1",   &status);  // arcsec, if dx=dy
  wcs->crval1   = drms_getkey_float (rec, "CRVAL1",   &status);  // disc center, arcsec
  wcs->crval2   = drms_getkey_float (rec, "CRVAL2",   &status);
  wcs->crpix1   = drms_getkey_float (rec, "CRPIX1",   &status);  // center, ccd, orig.
  wcs->crpix2   = drms_getkey_float (rec, "CRPIX2",   &status);
  wcs->crota2   = drms_getkey_float (rec, "CROTA2",   &status);
  wcs->crlt_obs = drms_getkey_float (rec, "CRLT_OBS", &status);
  wcs->crln_obs = drms_getkey_float (rec, "CRLN_OBS", &status);
  // assume set valid status...
  strcpy(wcs->missing, "");
  retval = 0;
  // ...and alter as necessary
  if (status != 0) {
    snprintf(wcs->missing, Nmiss, "Failed getkey()");
    retval = 1;
  } else {
    if (isnan(wcs->dsun_obs)) snprintf(wcs->missing, Nmiss, "NaN %s", "DSUN_OBS");
    if (isnan(wcs->cdelt1  )) snprintf(wcs->missing, Nmiss, "NaN %s", "CDELT1");
    if (isnan(wcs->crval1  )) snprintf(wcs->missing, Nmiss, "NaN %s", "CRVAL1");
    if (isnan(wcs->crval2  )) snprintf(wcs->missing, Nmiss, "NaN %s", "CRVAL2");
    if (isnan(wcs->crpix1  )) snprintf(wcs->missing, Nmiss, "NaN %s", "CRPIX1");
    if (isnan(wcs->crpix2  )) snprintf(wcs->missing, Nmiss, "NaN %s", "CRPIX2");
    if (isnan(wcs->crota2  )) snprintf(wcs->missing, Nmiss, "NaN %s", "CROTA2");
    if (isnan(wcs->crlt_obs)) snprintf(wcs->missing, Nmiss, "NaN %s", "CRLT_OBS");
    if (isnan(wcs->crln_obs)) snprintf(wcs->missing, Nmiss, "NaN %s", "CRLN_OBS");
    // if any of the above conditions fire, wcs->missing will be non-empty
    if (strlen(wcs->missing) > 0) 
      retval = 2;
  }
  if (retval) V_printf(1, "wcs_load: ", "WCS fail: %s\n", wcs->missing);
  return retval;
}

/*
 * open drms record, optionally load WCS, look-up segment, load segment
 *   i.e., general-purpose drms data-getter/checker
 * Usage: Call, for example, with query = 'hmi.M_720s' and trec = '2012.05.01_TAI'.
 *
 * Always: opens the corresponding record set, makes sure it has exactly one 
 *   matching record, and checks its quality for missing-ness (0x80000000).
 * Typically returns the RecordSet in *recSetP (it is returned because it must 
 *   be freed) and the single record in *recP (which becomes invalid once *recSetP 
 *   is freed).
 * Option 0a: If recP is NULL, finds but does not return recP.  
 * Option 0b: If, further to (0a) above, recSetP is NULL, finds but does not 
 *   return recSetP.
 * In either option 0 case, the error status is returned, allowing data-present 
 *   checks.  Note, you don't use 0b without 0a.
 * Option 1: If wcs is non-null, looks up and fills in wcs_t structure.  If this
 *   fails or WCS are partly missing, it's an error.
 * Option 2a: If segname is non-null, looks up the named segment 
 *   within the record, to verify it is there.  Nothing extra is returned.
 * Option 2b: If, further to (2a) above, dataP is non-null, loads 
 *   the segment according to d_typ, storing it in *dataP.  Again, *dataP 
 *   becomes invalid once *recSetP is freed.  
 *   The data in the *dataP is in (*dataP)->data.  
 *   If unneeded, you can let d_typ = DRMS_TYPE_RAW, which is noncommittal.
 *
 * Note:
 *   If you want the data (dataP) or the record (recP), you *have* to get the 
 *   recSet, because that's what you free.  If you want the WCS, you don't
 *   need to get the recSet.
 *   
 * If success, returns NULL.  If failure, returns a descriptive string, and in this
 *   case, nothing needs to be freed.  If success, be sure to:
 *     drms_close_records(*recSetP, DRMS_FREE_RECORD);
 *   when you are done, if you requested the recSet.
 */

static
char *
load_drms_rec_seg(DRMS_RecordSet_t **recSetP, 
          DRMS_Record_t       **recP, 
          DRMS_Array_t       **dataP, 
          wcs_t *wcs,
          // remainder are inputs
          char *query, 
          char *trec,
          char *segname,
          DRMS_Type_t d_typ)
{
  char rec_query[STR_MAX];
  char tag[STR_MAX];
  int  status, return_code;
  DRMS_RecordSet_t *the_recSet;
  DRMS_Record_t    *the_rec;
  DRMS_Array_t     *the_data;
  DRMS_Segment_t   *the_seg;
  
  // exclude mistakes
  if (!query || !trec)
    return "Illegal argument to load_drms_rec_seg: need query and trec";
  if (!recSetP && (recP || dataP)) 
    return "Illegal argument to load_drms_rec_seg: recSet required here";
  if (!segname && dataP)
    return "Illegal argument to load_drms_rec_seg: segname required for data";
  // null out output pointers as permitted, in case of early exit
  if (recSetP) *recSetP = NULL;
  if (recP)    *recP    = NULL;
  if (dataP)   *dataP   = NULL;
  if (wcs)     {wcs->ok = 0; strcpy(wcs->missing, ""); } // no message for now
  // descriptive tag
  snprintf(tag, sizeof(tag), "%s %s", 
       segname ? "and segment" : "(no segment)",
       dataP ? "and array" : "");
  // query for data
  snprintf(rec_query, sizeof(rec_query), "%s[%s]", query, trec);
  V_printf(verbflag > 2, "\t\t", "Opening record %s %s\n", rec_query, tag);
  status = 0;
  the_recSet = drms_open_records(drms_env, rec_query, &status);
  if (status != 0) {
    // (nothing to close)
    return "Could not open record at this T_REC";
  }
  // (below here, we close recSet before error out)
  if (the_recSet->n != 1) {
    drms_close_records(the_recSet, DRMS_FREE_RECORD);
    return "Could not open single record at this T_REC";
  }
  the_rec = the_recSet->records[0];
  // quality check to identify missing data 
  int quality = drms_getkey_int(the_rec, "QUALITY", &status);
  if (status != 0 || (quality & 0x80000000)) {
    drms_close_records(the_recSet, DRMS_FREE_RECORD);
    return "Missing data at this T_REC (quality)";
  }
  // option 1: load WCS
  if (wcs) {
    return_code = wcs_load(wcs, the_rec);
    if (return_code != 0) {
      drms_close_records(the_recSet, DRMS_FREE_RECORD);
      return "Missing WCS at this T_REC";
    }
  }
  // option 2a: look up segname
  if (segname) {
    the_seg = drms_segment_lookup(the_rec, segname);
    if (the_seg == NULL) {
      drms_close_records(the_recSet, DRMS_FREE_RECORD);
      return "Failed data segment lookup at this T_REC";
    }
  }
  // option 2b: data
  if (segname && dataP) {
    // (know the_seg is ok from above)
    the_data = drms_segment_read(the_seg, d_typ, &status);
    if (the_data == NULL || status != 0) {
      drms_close_records(the_recSet, DRMS_FREE_RECORD);
      return "Failed data segment read at this T_REC";
    }
  }
  // all OK: set up outputs as allowed
  if (recSetP) *recSetP = the_recSet;
  if (recP)    *recP    = the_rec;
  if (dataP)   *dataP   = the_data;
   
  // OTOH: close the record set if it was not returned
  if (!recSetP)
    drms_close_records(the_recSet, DRMS_FREE_RECORD);
  return NULL;
}

/* 
 * Fill in times in T_REC array
 *   key is: t0 is the time of tRec[nRec_pad]
 */
static
void
set_trec_times(trec_info_t *tRec, 
           TIME t0, 
           int nRec_all, int nRec_pad, 
           double cadence)
{
  TIME t;
  int rec;

  for (rec = 0; rec < nRec_all; rec++) {
    t = t0 + (rec - nRec_pad) * cadence;
    // fill in time fields
    tRec[rec].t = t;
    sprint_time(tRec[rec].str, t, "TAI", 0);
    // zero out some other fields
    tRec[rec].data = Trec_Data_Unknown;
    tRec[rec].nharp = 0;
  }
}


/*
 * mark each patch according to whether it is eligible to be ingested; 
 * then do the same for the tRec's, by looking over all patches at each time.
 */
static
void
mark_eligible_patches(patch_info_t *patchInfo, 
              harp_info_t  *harpInfo, 
              trec_info_t  *tRec, 
              int nHarp,
              int nRec_all, 
              int nRec_eat, 
              int rec0cut, 
              int rec1cut,
              int nRec_min)

{
  int rec, the_rec;
  int i, hits;

  V_printf(verbflag > 0, "\t", "Patch scan: finding patches eligible for ingest.\n");
  // mark everything eligible
  for (i = 0; i < nHarp; i++)
    for (rec = 0; rec < nRec_all; rec++)
      patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_Eligible;
  /* The order of consideration below determines which of several causes
   * is (later on) listed as "the" cause something was not ingested -- 
   * in each loop we check if the patch is eligible, and set it to 
   * ineligible (for whatever reason) if it fails the test.  Schematically:
   *   if (eligible && fails_foo_test) patch.ingest = ineligible_foo;
   * Thus, the first reason encountered becomes the winner.  So, the order
   * of the lines below is carefully chosen.
   */
  // "omit" case 1: nRec_eat limitation
  //    in NRT mode, everything but the last T_REC will be omitted
  if (nRec_eat < nRec_all)
    V_printf(verbflag > 0, "\t", "  Implementing trec-to-ingest limitation (e.g., NRT mode)\n");
  for (i = 0; i < nHarp; i++)
    for (rec = 0; rec < nRec_all-nRec_eat; rec++)
      if (patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible)
    patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_Omit;
  // "omit" case 2: t0cut limitation (if t0cut was nan, rec0cut = -1)
  //    now used in gap-filling
  if (rec0cut >= 0) {
    V_printf(verbflag > 0, "\t", "  Implementing t0cut limitation\n");
    // basically: mark from 0...rec0cut-1 as Omitted
    for (i = 0; i < nHarp; i++) {
      // if the padding straddles rec0cut, ingest the pad too
      // otherwise, omit ingesting at rec0cut as usual
      if (harpInfo[i].rec0p <= rec0cut && harpInfo[i].rec0 > rec0cut)
    the_rec = harpInfo[i].rec0p; // omit from 0..(here-1)
      else {
    the_rec = rec0cut; // omit from 0..(here-1)
    if (the_rec > nRec_all) the_rec = nRec_all; // happens if rec0cut is large
      }
      for (rec = 0; rec < the_rec; rec++)
    if (patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible)
      patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_Omit;
    }
  }
  // "omit" case 2: t1cut limitation (if t1cut was nan, rec1cut = -1)
  if (rec1cut >= 0) {
    V_printf(verbflag > 0, "\t", "  Implementing t1cut limitation\n");
    // basically: mark from rec1cut+1..nRec_all-1 as Omitted
    for (i = 0; i < nHarp; i++) {
      // if the padding straddles rec1cut, ingest the pad too
      // otherwise, omit ingesting at rec1cut as usual
      if (harpInfo[i].rec1p > rec1cut && harpInfo[i].rec1 <= rec1cut)
    the_rec = harpInfo[i].rec1p+1; // omit from here..top
      else
    the_rec = rec1cut+1; // omit from here..top
      // loop is OK even if the_rec > nRec_all
      for (rec = the_rec; rec < nRec_all; rec++)
    if (patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible)
      patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_Omit;
    }
  }
  // mark no-data patches (i.e., no mag or no mask)
  //  (patch is marked invalid already, so this must appear before 
  //   invalid section below)
  //  there are two kinds of data gap/missing conditions
  //     (1) no-data patch has rec0p <= rec <= rec1p --> list as gap/missing
  //     (2) no-data, and outside this interval --> list as Patch_Ingest_Invalid, below
  //  Note that valid == Patch_Invalid signals gap/missing AND outside-interval, so
  //  looking at Patch_Invalid is not enough
  for (rec = 0; rec < nRec_all; rec++) {
    if (tRec[rec].data == Trec_Data_Gap) { 
      // "gap" => mag missing
      for (hits = i = 0; i < nHarp; i++) {
    if ((patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible) &&
        ((rec >= harpInfo[i].rec0p) && (rec <= harpInfo[i].rec1p))) {
      patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_Gap;
      hits++;
    }
      }
      if (hits > 0)
    V_printf(verbflag > 0, "\t", "  Skipping %s (with %d HARP%s), no mag.\n", 
         tRec[rec].str, hits, hits == 1 ? "" : "s");
    } else if (tRec[rec].data == Trec_Data_Missing) {
      // "missing" => mag present, mask missing
      for (hits = i = 0; i < nHarp; i++) {
    if ((patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible) &&
        ((rec >= harpInfo[i].rec0p) && (rec <= harpInfo[i].rec1p))) {
      patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_Missing;
      hits++;
    }
      }
      if (hits > 0)
    V_printf(verbflag > 0, "\t", "  Skipping %s (with %d HARP%s), no mask.\n", 
         tRec[rec].str, hits, hits == 1 ? "" : "s");
    }
  }
  // mark in-valid patches (i.e., HARP did not extend that far)
  //  - note, padding/placeholder patches are Patch_Padding, not Patch_Invalid
  //  - needs to be done before CR exclusion, because a CR results in all HARP
  //    patches being marked ineligible
  for (rec = 0; rec < nRec_all; rec++)
    for (i = 0; i < nHarp; i++)
      if ((patchInfo[i*nRec_all + rec].valid  == Patch_Invalid) &&
      (patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible))
    patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_Invalid;
  // mark too-short patches (cosmic rays, mostly)
  for (i = 0; i < nHarp; i++) {
    // check: length=1  =>  rec1-rec0+1 = 1  => if() succeeds for nRec_min = 2
    if ((harpInfo[i].rec1 - harpInfo[i].rec0 + 1) < nRec_min) {
      V_printf(verbflag > 0, "\t", "  Skipping HARP[%d], too short (CR)\n", harpInfo[i].id);
      for (rec = 0; rec < nRec_all; rec++)
    if (patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible)
      patchInfo[i*nRec_all + rec].ingest = Patch_Ingest_TooShort;
    }
  }
  // count Eligible patches in each T_REC, given the above restrictions
  //   (wait until after ingestion to count per-harp patches)
  for (rec = 0; rec < nRec_all; rec++) {
    tRec[rec].nharp = 0;
    for (i = 0; i < nHarp; i++)
      if (patchInfo[i*nRec_all + rec].ingest == Patch_Ingest_Eligible) {
    tRec[rec].nharp++;
      }
  }
  V_printf(verbflag > 0, "\t", "Patch scan: done.\n");
}


/*
 *  for each T_REC, determine if there is data present; fill in wcs
 *  also finds and sets the full-disk image dimensions (i.e., NAXIS1), in *npix
 *  (since we examine all records, it this will be found if any record exists)
 */
static
void
set_trec_hasdata(trec_info_t *tRec, char *magQuery, char *maskQuery, int *npix, int nRec_all)
{
  int rec;
  wcs_t the_wcs;
  int Rec_valid, Rec_miss, Rec_err, Rec_ok;
  char *err;
  int npix_fulldisk = 0;

  V_printf(verbflag > 0, "\t", "T_REC scan: scanning image series for missing data\n");
  Rec_valid = Rec_miss = Rec_err = Rec_ok = 0;
  // loop over all records: we need status on everything to set HARP extent (e.g., T_FRST)
  // for each T_REC (of interest), determine if data is present
  for (rec = 0; rec < nRec_all; rec++) {
    Rec_valid++;
    // check for mag, load WCS
    err = load_drms_rec_seg(NULL, NULL, NULL, // NULL: do not want recSet or data
                &the_wcs,
                magQuery, tRec[rec].str, 
                NULL, DRMS_TYPE_RAW);
    if (err) {
      // mag data missing -- ordinary data gap, not an error, so no stderr output
      V_printf(1, "\t  ", "%s%s (%s[%s]) -- skipping this T_REC.\n", err, 
           (strlen(the_wcs.missing) > 0) ? the_wcs.missing : "",
           magQuery, tRec[rec].str);
      tRec[rec].data = Trec_Data_Gap; // no mag or mag wcs
      Rec_miss++;
      continue; 
    }
    // check for mask
    err = load_drms_rec_seg(NULL, NULL, NULL, // NULL: do not want recSet or data
                NULL,             // NULL: do not want WCS
                maskQuery, tRec[rec].str, 
                "mask", DRMS_TYPE_RAW);
    if (err) {
      // log to stdout and stderr
      V_printf( 1, "\t  ", "%s (%s[%s]) -- skipping this T_REC.\n", err, maskQuery, tRec[rec].str);
      V_printf(-1, "",   "%s (%s[%s]) -- skipping this T_REC.\n", err, maskQuery, tRec[rec].str);
      tRec[rec].data = Trec_Data_Missing; // mag, but no mask
      Rec_err++;
      continue; 
    } 
    if (npix_fulldisk == 0) {
      // find full-disk image dims -- we know the mask and mag are OK
      DRMS_RecordSet_t *mRecSet;
      DRMS_Array_t     *mArray = NULL;
      err = load_drms_rec_seg(&mRecSet, NULL,     // need recSet
                  &mArray,  NULL,     // do not want WCS
                  maskQuery, tRec[rec].str, 
                  "mask", DRMS_TYPE_RAW);
      // (do not expect errors, since the query just succeeded above)
        if (!err) {
            npix_fulldisk = mArray->axis[0];
            // Free the mArray.
            if (mArray)
            {
                drms_free_array(mArray);
                mArray = NULL;
            }
            
            drms_close_records(mRecSet, DRMS_FREE_RECORD);
        }
    }
    // data present
    tRec[rec].data = Trec_Data_OK;
    tRec[rec].wcs = the_wcs;
    Rec_ok++;
  }
  // set up the axis dimensions
  *npix = npix_fulldisk;
  V_printf(verbflag > 0, "\t", "T_REC scan: %d valid = %d OK + %d missing + %d errors\n", 
       Rec_valid, Rec_ok, Rec_miss, Rec_err);
  V_printf(verbflag > 1, "\t\t", "T_REC scan found NAXIS1 = %d\n", npix_fulldisk);
  // missing mask at a T_REC, when mag is present, is an error
  if (Rec_err > 0)
    V_printf(-1, "", "Found missing masks in %s at %d T_REC(s).  Continuing.\n", maskQuery, Rec_err);
}

/*
 * Summarize a patch by printing it -- diagnostic aid
 *
 * tI is a single tRec, pI is a single patch
 * s is a prefix string, "" for none
 */
static
void
patch_print(FILE *fp, char *s, trec_info_t *tI, patch_info_t *pI)
{
  const int human_friendly = 0;
  char trec_str[24];

  if (human_friendly) {
    if (pI->valid != Patch_Normal) {
      fprintf(fp, "%sPatch not valid\n", s);
    } else {
      sprint_time(trec_str, tI->t, "TAI", 0);
      fprintf(fp, 
          "%sID = %d (%s) @ (%.2f,%.2f) A = %.1f\n", 
          s, pI->num, trec_str, 
          pI->stats[RS_ar_area_lat],
          pI->stats[RS_ar_area_lon],
          pI->stats[RS_rgn_size]);
    }
  } else {
    // machine-readable output format
    // (time is #seconds since epoch, 1977/01/01)
    fprintf(fp, 
        "%d\t%d\t%.1f\t"
        "%.2f\t%.2f\t"
        "%.2f\t%.2f\t"
        "%.2f\t%.2f\t"
        "%.1f\n", 
        pI->num, pI->valid, (double) tI->t, 
        pI->stats[RS_ar_area_lat],
        pI->stats[RS_ar_area_lon],
        pI->lat0, pI->lon0,
        pI->lat1, pI->lon1,
        pI->stats[RS_rgn_size]);
  }
}

/*
 * load tracker parameters
 */
int
load_tracker_params(const char *rootDir, run_info_t *ri)
{
  FILE *fp;
  int p, np, llen;
  char fn[STR_MAX];
  char line[STR_MAX];
  char *val, *name_end;

  snprintf(fn, sizeof(fn), "%s/%s", rootDir, FN_TRACK_PARAM);
  if ((fp = fopen(fn, "r")) == NULL) {
    V_printf(-1, "", "Failed to open `%s'.\n", fn);
    return 1;
  }
  // find the number of parameters to be stored
  np = 0;
  while (fgets(line, sizeof(line), fp) != NULL) {
    llen = strlen(line);
    // we must have read all the way to the newline
    if (line[llen-1] != '\n') {
      V_printf(-1, "", "Line: `%.20s ...' in param file `%s' is too long.\n", line, fn);
      return 1;
    }
    if (line[0] == '#' || line[0] == '\n') continue;
    np++;
  }
  ri->n_par = np;
  if (np == 0) {
    fclose(fp);
    return 0;
  }
  // allocate space
  ri->par_names = calloc(np, sizeof(*(ri->par_names)));
  ri->par_vals  = calloc(np, sizeof(*(ri->par_vals )));
  if (ri->par_names == NULL || ri->par_vals == NULL) {
    V_printf(-1, "", "Calloc fail in load_tracker_params for np = %d\n", np);
    fclose(fp);
    return 1;
  }
  // read params
  rewind(fp);
  p = 0;
  while (fgets(line, sizeof(line), fp) != NULL) {
    llen = strlen(line);
    if (line[0] == '#' || line[0] == '\n') continue;
    // we know it's newline terminated, so drop the newline
    line[--llen] = '\0';
    // drop any other trailing whitespace
    while (line[llen-1] == ' ' || line[llen-1] == '\t')
      line[--llen] = '\0';
    // find the :
    val = strchr(line, (int) ':');
    if (val == NULL) {
      V_printf(-1, "", "load_tracker_params: no colon in line `%s' in file `%s'\n", 
           line, fn);
      fclose(fp);
      return 1;
    }
    // trim trailing whitespace before : (if any)
    for (name_end = val-1; *name_end == ' ' || *name_end == '\t'; name_end--)
      *name_end = '\0';
    *val++ = '\0'; // null the colon, thus terminating *name
    // null any following whitespace
    while (*val == ' ' || *val == '\t')
      *val++ == '\0';
    // copy params
    ri->par_names[p] = strdup(line);
    ri->par_vals [p] = strdup( val);
    p += 1; // next param
  }
  fclose(fp);
  V_printf(verbflag > 0, "", "Got %d tracker params from `%s'\n", np, fn);
  return 0; // OK
}

/* 
 * Fetch master index array summary info, return 0 if successful 
 */
static
int 
load_run_info(const char *rootDir, char **lists, 
          double cadence,
          int *nHarp, 
          int *nRec, 
          TIME *p_t0,  // pointer to t0 
          TIME *p_t1)  // pointer to t1
{
  char *list1;
  char full_list[STR_MAX], currDir[STR_MAX], frameFile[STR_MAX];
  char line[STR_MAX];
  char *endptr;
  FILE *fp, *frameFp;
  int harp_id;       // external HARP ID number
  int harp_count;    // number of HARPs read (good and bad)
  int line_num;      // friendly error messages
  int rv, status;    // error codes
  patch_info_t pInfo;
  TIME t, t0, t1;    // time, and end times
  int *harp_ids;
  int nHarp_top = 1024; // size of above array
  int dupe;
  int j;

  harp_ids = calloc(nHarp_top, sizeof(*harp_ids));
  if (harp_ids == NULL) {
    V_printf(-1, "", "Error: Unable to allocate space for harp-id list\n");
    return 1;
  }
  harp_count = 0;
  t0 = sscan_time("2999.12.31_23:59:28Z");  // earliest time seen
  t1 = sscan_time("-4712.01.01_11:59:28Z"); // latest time seen
  while ((list1 = *lists++) != NULL) {
    // read from list1
    if (*list1 == '/')
      snprintf(full_list, sizeof(full_list), "%s", list1); // rooted filename
    else
      snprintf(full_list, sizeof(full_list), "%s/%s", rootDir, list1);
    fp = fopen(full_list, "r");
    if (!fp) {
      V_printf(-1, "", "Error: Could not open HARP list file %s\n", full_list);
      return 1;
    }
    status = 0; // OK
    while (fgets(line, sizeof(line), fp)) {
      if (*line == '#') continue; // comment line
      harp_id = strtol(line, &endptr, 10);
      // check for errors here
      if (harp_id <= 0 || harp_id > INT_MAX) {
    // 0 is the error code for strtol
    V_printf(-1, "", "Bad HARP number format in %s: %s\n", full_list, line);
    status = 1;
    break;
      } 
      // skip if duplicated
      // makes it O(n^2), but n is small
      for (dupe = j = 0; !dupe && j < harp_count; j++)
    if (harp_ids[j] == (int) harp_id)
      dupe = 1;
      // updates count (which is returned) as well as id list (which is not)
      if (!dupe)
    harp_ids[harp_count++] = (int) harp_id;
      else
    V_printf(verbflag > 1, "  ", "ID count: Skipped dupe HARP #%d (this is ok).\n", 
         (int) harp_id);
      // resize if needed
      if (harp_count == nHarp_top) {
    nHarp_top *= 4;
    if ((harp_ids = realloc(harp_ids, nHarp_top*sizeof(*harp_ids))) == NULL) {
      V_printf(-1, "", "Failed realloc of harp ID buffer to size %d\n", nHarp_top);
      status = 1;
      break;
    }
      }
      // HARP was OK, now scan for its time range
      // using implicit string concatenation below to construct format
      snprintf(frameFile, sizeof(frameFile), 
           "%s/" FN_TRACK_DIR "/%s", rootDir, harp_id, FN_TRACK_FRAME);
      // open the frame info for this HARP
      frameFp = fopen(frameFile, "r");
      if (!frameFp) {
    V_printf(-1, "", "Error: frame list file %s not found for HARP %06d (prescan).\n", 
         frameFile, harp_id);
    status = 1;
    break;
      }
      // read the frame file, just to find the time range
      for (line_num = 1; 
       (rv = load_frame(frameFp, &pInfo, &t)) == 0; 
       line_num++) {
    if (t < t0) t0 = t;
    if (t > t1) t1 = t;
      }
      fclose(frameFp);
      // this value signals format error
      if (rv == 1) {
    V_printf(-1, "", "Error: Bad frame line %d in %s.\n", line_num, frameFile);
    status = 1; // push error upward
    break;
      }
    }
    fclose(fp);
    if (status)
      return 1; // Error: fail out after fclose()
  } // end loop over lists
  free(harp_ids);

  // No errors: set up outputs
  // put correct nRec, and a sensible time in, even if there were no HARPs
  if (harp_count > 0) {
    float nRec_check = ((t1 - t0) / cadence) + 1; // number of records
    *p_t0 = t0;
    *p_t1 = t1;
    *nHarp = harp_count;
    *nRec = (int) nRec_check;
    if (abs(nRec_check - *nRec) > REC_EPS/cadence) {
      V_printf(-1, "", "Error: cadence (%f) does not divide HARP TREC range (%f).\n", 
           cadence, t1-t0);
      return 1;
    }
  } else {
    *nHarp = 0;
    *nRec = 0;
    *p_t0 = *p_t1 = sscan_time("2011.02.14_00:00:00_TAI");
  }
  // printf("nHarp = %d, nRec = %d, t0 = %f, t1 = %f\n", *nHarp, *nRec, t0, t1);
  return 0;
}


/* 
 * Get HARP ID field in HarpInfo structure, return 0 if successful 
 */
static 
int 
load_harp_ids(const char *rootDir, char **lists, harp_info_t *harpInfo)
{
  char *list1;
  char full_list[STR_MAX];
  char line[STR_MAX];
  char *endptr;
  FILE *fp;
  long harp_id; // external HARP ID number
  int i, j;     // internal HARP index (starts at 0)
  int dupe;     // duplicate harp ID found?
    
  i = 0; // retain across list files
  while ((list1 = *lists++) != NULL) {
    // read list1
    if (*list1 == '/')
      snprintf(full_list, sizeof(full_list), "%s", list1); // rooted filename
    else
      snprintf(full_list, sizeof(full_list), "%s/%s", rootDir, list1);
    fp = fopen(full_list, "r");
    if (!fp) {
      V_printf(-1, "", "Could not open HARP summary file %s\n", full_list);
      return 1;
    }
    while (fgets(line, sizeof(line), fp)) {
      if (*line == '#') continue; // comment line
      // this error check was already made by getRunInfo
      harp_id = strtol(line, &endptr, 10);
      if (harp_id <= 0 || harp_id > INT_MAX) {
    // 0 is the error code
    V_printf(-1, "", "Bad HARP ID format in %s: %s\n", full_list, line);
    fclose(fp);
    return 1;
      } 
      // skip if duplicated
      // makes it O(n^2), but n is small
      for (dupe = j = 0; !dupe && j < i; j++)
    if (harpInfo[j].id == (int) harp_id)
      dupe = 1;
      if (!dupe)
    harpInfo[i++].id = (int) harp_id;
      else
    V_printf(verbflag > 2, "  ", "ID load: Skipped dupe HARP #%d (this is ok)\n", 
         (int) harp_id);
    }
    fclose(fp);
  }
  return 0;
}

/*
 * load reference bounding box from HARP
 * tri-valued return code:
 *    (1) NULL if box was found and is OK,
 *    (2) "" if no reference box found, or 
 *    (3) error message if reference box was missing/erroneous
 */
static
char *
load_harp_reference(const char *refHarp, int harp_id, TIME *t0, 
            double *lat0, double *lon0, double *lat1, double *lon1, double *omega)
{
  char harpRecQuery[STR_MAX];
  DRMS_RecordSet_t *harpRecSet;
  DRMS_Record_t *harpRec;
  int status = 0;
  
  // e.g., hmi.Mharp_720s[1234][^]
  snprintf(harpRecQuery, sizeof(harpRecQuery), "%s[%d][%s]", refHarp, harp_id, "^");
  V_printf(verbflag > 2, "\t\t", "Opening HARP %s (box match)\n", harpRecQuery);
  harpRecSet = drms_open_records(drms_env, harpRecQuery, &status);
  if (status != 0)
    return ""; // reference not found: not an error
  if (harpRecSet->n != 1) 
    return "Found more than one HARP record";
  harpRec = harpRecSet->records[0];
  *lat0  = drms_getkey_double(harpRec, "LATDTMIN", &status);
  *lon0  = drms_getkey_double(harpRec, "LONDTMIN", &status);
  *lat1  = drms_getkey_double(harpRec, "LATDTMAX", &status);
  *lon1  = drms_getkey_double(harpRec, "LONDTMAX", &status);
  *omega = drms_getkey_double(harpRec, "OMEGA_DT", &status);
  *t0    = drms_getkey_time  (harpRec, "T_REC",    &status);
  if (status != 0)
    return "Could not get HARP keys";
  if (isnan(*lat0) || isnan(*lon0) || isnan(*lat1) || isnan(*lon1) || isnan(*omega))
    return "HARP bbox was invalid/missing";
  if (time_is_invalid(*t0))
    return "HARP time was invalid/missing";
  // everything OK
  return NULL;
}


/*
 * patch_set_to_harp_reference: set patch bounding box information, given a reference
 * box and rotational velocity
 */
static
void
patch_set_to_harp_reference(patch_info_t *pI, TIME t, TIME t0, 
                double lat0, double lon0, double lat1, double lon1, double omega)
{
  // units: a "TIME" is a double, in seconds
  TIME dt = t - t0; // always >= 0
  // units: basic omega is deg/day -- omegaP is deg/sec
  double omegaP = omega / SEC_PER_DAY;

  pI->lat0 = lat0;
  pI->lon0 = lon0 + omegaP * dt; 
  pI->lat1 = lat1;
  pI->lon1 = lon1 + omegaP * dt; 
  pI->omega = omega;
}


/*
 * Loads the patches from a single HARP
 *   Returns nonzero for failure
 */
static
int 
load_harp_patches(int harp_id,
          const char *refHarp,      // if non-null: source for HARP boxes
          char *currDir,
          char *frameFile,
          char *statFile,
          patch_info_t *patchInfo,  // head of patches for this HARP
          harp_info_t  *harpInfo,   // info for just this HARP
          trec_info_t  *tRec,       // head of full tRec array
          int nRec_all,
          double cadence)
{
  // reference bounding box info
  double lat0, lon0, lat1, lon1, omega; // not float, for more digits
  TIME t0;
  char t0_S[32];  // time string
  int ref_ready;
  // iteration and place-keeping
  FILE *frameFp, *statFp;
  int rec, rec0, rec1;    // first and last rec#, plus counter
  TIME t, t_rec0, t_rec1; // time for rec0, rec1
  patch_info_t pInfo;     // place for one frame's info
  int rv1, rv2, status;   // return values, error status
  int line_num, nPatch;   // counters
  char the_patchName[STR_MAX];

  // open files
  frameFp = fopen(frameFile, "r");
  if (!frameFp) {
    V_printf(-1, "", "Frame list file `%s' not found for HARP %06d (full scan).\n", 
         frameFile, harp_id);
    return 1;
  }
  statFp = fopen(statFile, "r");
  if (!statFp) {
    fclose(frameFp);
    V_printf(-1, "", "Stat list file `%s' not found for HARP %06d\n", 
         statFile, harp_id);
    return 1;
  }

  // attempt to get existing HARP box info, if it was asked for
  ref_ready = 0;  // assume it was not there
  if (refHarp) {
    char *msg;
    msg = load_harp_reference(refHarp, harp_id, &t0, &lat0, &lon0, &lat1, &lon1, &omega);
    if (msg == NULL) {
      ref_ready = 1;  // box found and was OK
      sprint_time(t0_S, t0, "TAI", 0); // time string
      V_printf(verbflag > 2, "\t\t", 
           "Reference box: %s[%d][%s]: (%.4f,%.4f) to (%.4f,%.4f)\n", 
           refHarp, harp_id, t0_S, lat0, lon0, lat1, lon1);
    } else {
      strcpy(t0_S, ""); // clear time string
      if (*msg == '\0')
    // simple failure to open, not alarming
    V_printf(verbflag > 2, "\t\t", 
         "Could not open reference harp %s[%d] -- no such harp.\n", 
         refHarp, harp_id);
      else
    // an erroneous box, a little more surprising
    V_printf(verbflag > 1, "\t\t", 
         "Error opening reference harp %s[%d] -- bad box (%s).\n", 
         refHarp, harp_id, msg);
    }
  }

  // read files line-by-line in lockstep
  for (status = 0, nPatch = 1, line_num = 1; 
       status == 0 &&
     (rv1 = load_frame(frameFp, &pInfo, &t)) == 0 &&
     (rv2 = load_stats(statFp, pInfo.stats)) == 0; 
       nPatch++, line_num++) {
    // generate filename, start from 000001
    // format is embedded using literal string concatenation
    snprintf(the_patchName, sizeof(the_patchName), 
         "%s/" FN_TRACK_FITS, currDir, nPatch);
    if ((pInfo.patchName = strdup(the_patchName)) == NULL) {
      V_printf(-1, "", "Error: failed strdup (%s), harp %d\n", the_patchName, harp_id);
      status = 1;
      continue; // will error out
    }
    double rec_check = (t - tRec[0].t) / cadence;  // should be integral
    rec = (int) rec_check;
    // check for in-range rec
    if (rec < 0 || rec >= nRec_all) {
      V_printf(-1, "", "Error: time at line %d out-of-range\n", line_num);
      status = 1;
      continue; // will error out
    }
    // rec must be an even multiple of cadence
    if (abs(rec - rec_check) > REC_EPS/cadence) {
      V_printf(-1, "", "Error: Mismatch between time at line %d, "
           "initial time, and cadence %f\n", line_num, cadence);
      status = 1;
      continue; // will error out
    }
    // set bounding box from reference BB, if desired and possible
    if (ref_ready)
      patch_set_to_harp_reference(&pInfo, t, t0, lat0, lon0, lat1, lon1, omega);

    // insert the patch into its place in the master index
    //  (gaps are allowed and will be left as patch_invalid)
    //  (this "=" duplicates pInfo into patchInfo)
    patchInfo[rec] = pInfo;
  }
  fclose(frameFp);
  fclose(statFp);
  // error out: file trouble
  if (rv1 == 1) {
    V_printf(-1, "", "Error: Bad frame line %d in `%s'.\n", line_num, frameFile);
    status = 1;
  } else if (rv2 == 1) {
    V_printf(-1, "", "Error: Bad stats line %d in `%s'.\n", line_num, statFile);
    status = 1;
  } else if (status) {
    V_printf(-1, "", "Error: Bad line %d in `%s'.\n", line_num, frameFile);
  }
  return status;
}

/* 
 * Fill in blanks in patches of one HARP -- padding and placeholders
 *   (1) interpolates placeholder ROIs into HARP proper
 *   (2) extrapolates padding on either end of HARP
 *   Only does (1) and (2) where data is present
 *   returns 0 if successful, nonzero if not
 *   (currently, always successful)
 */ 
static
int
set_patch_info_blanks(patch_info_t *patchInfo,    // head of array this HARP's patches
              harp_info_t  *harpInfo,     // this HARP's info
              trec_info_t  *tRec, 
              int nRec_all,
              double cadence)
{
  int rec, sn;                         // counters
  const int rec0p = harpInfo->rec0p;   // first rec within pad
  const int rec0  = harpInfo->rec0;    // first rec with a real patch
  const int rec1  = harpInfo->rec1;    // last such rec
  const int rec1p = harpInfo->rec1p;   // last rec within pad
  const float omega =                  // rotation (deg/frame) -- same for all rec
    patchInfo[rec0].omega / SEC_PER_DAY * cadence;

  // printf("\tHarp %d: recs=[%d, %d]\n", harpInfo->id, rec0, rec1);
  // put in lat/lon for placeholder ROIs, which always occur between rec0 and rec1
  for (rec = rec0; rec <= rec1; rec++) {
    // install patchInfo[rec] if we have data, but saved no patch
    if (tRec[rec].data == Trec_Data_OK && patchInfo[rec].valid == Patch_Invalid) {
      // it will be invalid if it was never initialized
      patchInfo[rec].valid = Patch_Padding;   // placeholder ROI
      patchInfo[rec].tag   = Patch_Tag_Faint; // indicate HARP not observed
      patchInfo[rec].num   = patchInfo[rec0].num;
      patchInfo[rec].lat0  = patchInfo[rec0].lat0;
      patchInfo[rec].lat1  = patchInfo[rec0].lat1;
      patchInfo[rec].lon0  = patchInfo[rec0].lon0 + (rec - rec0)*omega;
      patchInfo[rec].lon1  = patchInfo[rec0].lon1 + (rec - rec0)*omega;
      patchInfo[rec].omega = patchInfo[rec0].omega;
      // Don't make a pointer to static memory - previously this was malloc'd.
      // patchInfo[rec].patchName = "";   // not a real patch
      if (patchInfo[rec].patchName)
      {
          free(patchInfo[rec].patchName);
          patchInfo[rec].patchName = calloc(1, sizeof(char));   // not a real patch
      }

      for (sn = 0; sn < RS_num_stats; sn++)
    patchInfo[rec].stats[sn] = DRMS_MISSING_FLOAT;
    }
  }
  // extend backward -- base lon offset on rec0
  for (rec = rec0p; rec < rec0; rec++) {
    // set up patchInfo[rec] if we have mask/mag data
    if (tRec[rec].data == Trec_Data_OK) {
      patchInfo[rec].valid = Patch_Padding;  // extrapolated
      patchInfo[rec].tag   = Patch_Tag_None; // no tag
      patchInfo[rec].num   = patchInfo[rec0].num;
      patchInfo[rec].lat0  = patchInfo[rec0].lat0;
      patchInfo[rec].lat1  = patchInfo[rec0].lat1;
      patchInfo[rec].lon0  = patchInfo[rec0].lon0 + (rec - rec0)*omega;
      patchInfo[rec].lon1  = patchInfo[rec0].lon1 + (rec - rec0)*omega;
      patchInfo[rec].omega = patchInfo[rec0].omega;
      // Don't make a pointer to static memory - previously this was malloc'd.
      // patchInfo[rec].patchName = "";   // not a real patch
      if (patchInfo[rec].patchName)
      {
          free(patchInfo[rec].patchName);
          patchInfo[rec].patchName = calloc(1, sizeof(char));   // not a real patch
      }
        
      for (sn = 0; sn < RS_num_stats; sn++)
    patchInfo[rec].stats[sn] = DRMS_MISSING_FLOAT;
    }
  }
  // extend forward (count down) -- base lon offset on rec1
  for (rec = rec1p; rec > rec1; rec--) {
    // set up patchInfo[rec] if we have mask/mag data
    if (tRec[rec].data == Trec_Data_OK) {
      patchInfo[rec].valid = Patch_Padding;  // extrapolated
      patchInfo[rec].tag   = Patch_Tag_None; // no tag
      patchInfo[rec].num   = patchInfo[rec1].num;
      patchInfo[rec].lat0  = patchInfo[rec1].lat0;
      patchInfo[rec].lat1  = patchInfo[rec1].lat1;
      patchInfo[rec].lon0  = patchInfo[rec1].lon0 + (rec - rec1)*omega;
      patchInfo[rec].lon1  = patchInfo[rec1].lon1 + (rec - rec1)*omega;
      patchInfo[rec].omega = patchInfo[rec1].omega;
        
      // Don't make a pointer to static memory - previously this was malloc'd.
      // patchInfo[rec].patchName = "";   // not a real patch
      if (patchInfo[rec].patchName)
      {
          free(patchInfo[rec].patchName);
          patchInfo[rec].patchName = calloc(1, sizeof(char));   // not a real patch
      }
      for (sn = 0; sn < RS_num_stats; sn++)
    patchInfo[rec].stats[sn] = DRMS_MISSING_FLOAT;
    }
  }
  return 0;  // always OK
}


/* 
 * set up one HARP's temporal extent info: rec0,rec1 and rec0p,rec1p
 *   the rec0p/rec1p settings are idealized; we refine them later after we account
 *     for farside geometry
 *   does not set up the number of valid patches in-between; we do that after
 *     the farside check
 *   note, HARP id field is already set
 *   returns 0 if successful, nonzero if not
 */ 
static
int
set_harp_info_extent(patch_info_t *patchInfo, // offset to this HARP's patches
             harp_info_t *hInfo,      // offset to this HARP
             trec_info_t *tRec,       // all t_rec's
             int nRec_all, 
             int nRec_pad)
{
  int rec;                // T_REC counter
  int rec0,  rec1;        // first and last rec# (normal patches)
  int rec0p, rec1p;       // first and last rec# (padding patches)

  // find range of rec numbers, {rec0...rec1} and {rec0p...rec1p}, for this harp
  rec0 = nRec_all - 1;
  rec1 = 0; 
  for (rec = 0; rec < nRec_all; rec++) {
    if (patchInfo[rec].valid == Patch_Normal) {
      if (rec < rec0) rec0 = rec;
      if (rec > rec1) rec1 = rec;
    }
  }
  // the "empty HARP" case will cause trouble later
  if (rec0 > rec1) {
    V_printf(-1, "", "Error: No patches present in this HARP.\n");
    return 1;
  }
  // account for padding at front
  rec0p = rec0 - nRec_pad;
  if (rec0p < 0) rec0p = 0;
  // advance rec0p until it points to a T_REC having data
  while (tRec[rec0p].data != Trec_Data_OK) 
    rec0p++;
  // account for padding at end
  rec1p = rec1 + nRec_pad;
  if (rec1p >= nRec_all) rec1p = nRec_all-1;
  // drop rec1p until it points to a T_REC having data
  while (tRec[rec1p].data != Trec_Data_OK) 
    rec1p--;
  // Cannot count missing patches because extrapolation has not been done yet
  V_printf(verbflag > 2, "\t\t", 
       "HARP %d: recs = {%d..%d} recspad = {%d..%d}\n", 
       hInfo->id, rec0, rec1, rec0p, rec1p);
  hInfo->rec0  = rec0;
  hInfo->rec1  = rec1;
  hInfo->rec0p = rec0p;
  hInfo->rec1p = rec1p;
  return 0;
}

/* 
 * set up one HARP's pixel-box info: {x,y}{min,max} from existing lat-lon info
 *
 *   the pixel-box is used later to make a bitmap to ingest
 *   additionally, at this time we determine if the HARP is on the far side;
 *     if so we mark it as invalid, and adjust rec0p/rec1p if the invalid
 *     times cause a shrinkage of the interval
 *   rec0/rec1 will never need to be adjusted because they are not on far side
 *   pix_max is the full-disk image size
 */ 
static
void
set_patch_info_pixbox(patch_info_t *patchInfo, // offset to this HARP's patches
              harp_info_t *hInfo,      // offset to this HARP
              trec_info_t *tRec,       // all t_rec's
              int pix_max, 
              int nRec_all)
{
  int rec;
  const int rec0p = hInfo->rec0p;   // first rec within pad
  const int rec0  = hInfo->rec0;    // first rec with a real patch
  const int rec1  = hInfo->rec1;    // last such rec
  const int rec1p = hInfo->rec1p;   // last rec within pad
  int rv, farside, pad_lo, pad_hi;

  // convert all lat/lon boxes to x/y boxes, collect farside info
  farside = 0;
  for (rec = rec0p; rec <= rec1p; rec++) {
    // Invalid's in this range can occur for missing mag or mask
    if (patchInfo[rec].valid == Patch_Invalid)
      continue;
    // Find x/y boundary of HARP at T_REC using lat0/lon0, lat1/lon1 info
    //    rv = 0 for OK, 1 for all pix on farside
    V_printf(verbflag > 2, "\t\t", "HARP[%d][%s]: bbox update\n",
         hInfo->id, tRec[rec].str);
    rv = compute_boundary(&(patchInfo[rec]), pix_max, tRec[rec].wcs);
    // check farside case
    if (rv) {
      // entire HARP was on farside (should only happen for extrapolated HARPs)
      farside++;
      // set to in-valid, to ignore the patch
      //   note, once this happens, rec0p/rec1p may not point 
      //   to the true end of the padding -- we correct them below
      patchInfo[rec].valid = Patch_Invalid;
    } 
  }

  // account for (potential) vanishing padding at front
  //   advance rec until it points to a valid patch 
  //   (second && should not be needed because patch at rec0 is always valid)
  pad_lo = 0;
  for (rec = rec0p; patchInfo[rec].valid == Patch_Invalid && rec < rec0; rec++)
    pad_lo++;
  hInfo->rec0p = rec;

  // account for vanishing padding at end
  //   drop rec until it points to a valid patch
  pad_hi = 0;
  for (rec = rec1p; patchInfo[rec].valid == Patch_Invalid && rec > rec1; rec--)
    pad_hi++;
  hInfo->rec1p = rec;

  if (pad_lo > 0 || pad_hi > 0)
    V_printf(verbflag > 2, "\t\t", 
         "HARP %d: recspad shrinks to {%d..%d} after finding %d farside patches\n", 
         hInfo->id, hInfo->rec0p, hInfo->rec1p, farside);
}


/* 
 * set up one HARP's missing patch info
 *   note, HARP boundary info is already set
 *   takes data availability in tRec into account
 *   returns 0 if successful, nonzero if not
 */ 
static
int
set_harp_info_miss(patch_info_t *patchInfo, // offset to this HARP's patches
           harp_info_t  *hInfo,     // offset to this HARP
           trec_info_t  *tRec,      // not offset, of course
           int nRec_all)
{
  int rec;                          // T_REC counter
  int n_missing;                    // count missing patches
  int n_patch, n_patchp;            // count present patches
  const int rec0p = hInfo->rec0p;   // first rec within pad
  const int rec0  = hInfo->rec0;    // first rec with a real patch
  const int rec1  = hInfo->rec1;    // last such rec
  const int rec1p = hInfo->rec1p;   // last rec within pad

  // count present and missing T_REC's within the range {rec0..rec1}
  for (n_patch = n_missing = 0, rec = rec0; rec <= rec1; rec++) {
    if (tRec[rec].data == Trec_Data_OK && patchInfo[rec].valid != Patch_Invalid)
      n_patch++;
    if (tRec[rec].data != Trec_Data_OK || patchInfo[rec].valid == Patch_Invalid)
      n_missing++;
  }
  // count all present T_REC's within the range {rec0p..rec1p}
  for (n_patchp = 0, rec = rec0p; rec <= rec1p; rec++)
    if (tRec[rec].data == Trec_Data_OK && patchInfo[rec].valid != Patch_Invalid)
      n_patchp++;
  V_printf(verbflag > 2, "\t\t", 
       "HARP %d: #patches-seen = %d, #missing = %d, #patches-with-padding = %d\n", 
       hInfo->id, n_patch, n_missing, n_patchp);
  hInfo->n_missing = n_missing;
  hInfo->n_patch   = n_patch;
  hInfo->n_patchp  = n_patchp;
  return 0;
}

/* 
 * Mark a HARP as ingested: returns 0 if OK
 */
static
int 
mark_ingested_harp(const char *rootDir, 
           run_info_t *runInfo,
           harp_info_t *hInfo)
{
  char statusFile[STR_MAX];
  FILE *statusFp;
  int harp_id;
  time_t tloc;
  char timestr[32]; // enough for an ISO date
  const char *run_key = "run_name"; // in runInfo, from tracker_params
  char *run_name; // "" if not found in runInfo
  int p;

  if (rootDir == NULL) return 1; // "readOnly" option -- will cause log message
  harp_id = hInfo->id;
  // ISO date, like: 20120217T102626
  time(&tloc);
  strftime(timestr, sizeof(timestr), "%Y%m%dT%H%M%S", localtime(&tloc));
  // run name: find a parameter p matching the name in run_key
  run_name = "unnamed_run"; // in case it's not found
  for (p = 0; p < runInfo->n_par; p++) {
    if (strcasecmp(run_key, runInfo->par_names[p]) == 0) {
      run_name = runInfo->par_vals[p];
      break;
    }
  }
  // Name of HARP status file
  // using implicit string concatenation below to construct format
  snprintf(statusFile, sizeof(statusFile), 
       "%s/" FN_TRACK_DIR "/%s", rootDir, harp_id, FN_TRACK_STATUS);
  // append to status info for this HARP
  statusFp = fopen(statusFile, "a");
  if (!statusFp) {
    V_printf(-1, "", 
         "Warning: unable to append to status file %s for HARP %06d (post-ingest).\n", 
         statusFile, harp_id);
    return 1; // probably non-fatal
  }
  fprintf(statusFp, "%s\t%s\t%s\n", "ingested", run_name, timestr);
  fclose(statusFp);
  return 0; // OK
}


/* 
 * Most complete routine for loading patch information
 * Runs after we know how many T_REC's and HARPS we have
 *   returns 0 if successful, nonzero if not
 */ 
static
int 
load_all_patch_info(const char *rootDir, 
            const char *refHarp,
            patch_info_t *patchInfo, 
            harp_info_t  *harpInfo, 
            trec_info_t  *tRec, 
            int nHarp, 
            int nRec_all, int nRec_pad, 
            double cadence,
            int pix_max)
{
  // file info
  char currDir[STR_MAX];
  char frameFile[STR_MAX], statFile[STR_MAX];
  // harp and record info
  int i;                  // harp counter
  int harp_id;            // external harp id number
  int rec;                // rec# counter
  patch_info_t pInfo;     // place for one frame's info
  int status;             // return values, error status
  
  for (i = 0; i < nHarp; i++) {
    //
    // Read the per-HARP info for the i'th track given to us
    //
    // assemble file names
    harp_id = harpInfo[i].id; // its external ID
    // using string concatenation to insert dirname format
    snprintf(currDir, sizeof(currDir), "%s/" FN_TRACK_DIR, rootDir, harp_id);
    snprintf(frameFile, sizeof(frameFile), "%s/%s", currDir, FN_TRACK_FRAME);
    snprintf(statFile,  sizeof(statFile),  "%s/%s", currDir, FN_TRACK_STATS);
    // Load all patches for this HARP
    status = 
      load_harp_patches(harp_id,
            refHarp,
            currDir, frameFile, statFile,
            patchInfo + i*nRec_all, // offset to HARP #i
            harpInfo + i,           // offset to HARP #i
            tRec,                   // not offset
            nRec_all, cadence);
    if (status != 0) {
      V_printf(-1, "", "Error: Could not load HARP %d (id = %06d).\n", i, harp_id);
      return status;
    }
    // Put per-harp temporal extent information into harpInfo[i] -- rec0/rec1, etc.
    //   rec0p/rec1p are refined later, after we have farside info
    status = set_harp_info_extent(patchInfo+nRec_all*i, harpInfo+i, tRec, nRec_all, nRec_pad);
    if (status != 0) {
      V_printf(-1, "", "Error: Could not find boundaries: HARP %d (id = %06d).\n", 
           i, harp_id);
      return status;
    }
    // Fill in blank patches according to harp boundary above -- padding and placeholders
    set_patch_info_blanks(patchInfo+nRec_all*i, harpInfo+i, tRec, nRec_all, cadence);

    // Set up pixel bounding box of each patch in the harp
    //   - uses per-tRec geometry and patch lat/lon bounding box determined above
    //   - as a side-effect, determines which patches were on far side, sets
    //     them back to invalid, and shrinks the rec0p/rec1p interval as needed
    set_patch_info_pixbox(patchInfo+nRec_all*i, harpInfo+i, tRec, pix_max, nRec_all);

    // Put per-harp missing-patch information into harpInfo[i] -- n_patch, etc.
    set_harp_info_miss(patchInfo+nRec_all*i, harpInfo+i, tRec, nRec_all);
  }
  // testing
  if (0)
    for (i = 0; i < 1 && i < nHarp; i++)
      for (rec = 0; rec < nRec_all; rec++) {
    if (patchInfo[i * nRec_all + rec].valid) {
      patch_info_t tmp = patchInfo[i * nRec_all + rec];
      if (tmp.valid != Patch_Normal) continue;
      printf("\trec# = %d, x0 = %d, y0 = %d\n", rec, tmp.x0, tmp.y0);
      printf("\ttrack stats: #pix=%.2f, area1=%.2f, area2=%.2f\n", 
         tmp.stats[0], tmp.stats[1], tmp.stats[2]);
      fflush(stdout);
      break;
    }
      }
  return 0;
}



/* 
 * ingest one harp at a particular T_REC
 * Returns 0 if successful, 1 if not
 * Some aspects of ingestion may fail in minor ways.  This still results
 * in a "success" exit code.  Two examples: failure to insert some keywords,
 * failure to insert link.
 * Note, HARP rotating off the disk used to be a minor failure, but we catch
 * this condition earlier now.
 */
static
int 
ingest_harp(patch_info_t *pInfo,     // the patch to ingest, function of (HARP,T_REC)
        harp_info_t  *hInfo,     // info for the HARP
        trec_info_t  *tRec1,     // info for the T_REC
        trec_info_t  *tRec,      // info for *all* T_REC's
        run_info_t   *runInfo,   // info for whole run
        marp_info_t  *mInfo,     // matching ARs
        DRMS_Record_t *magRec,   // mag[T_REC]
        DRMS_Record_t *maskRec,  // mask[T_REC]
        char *maskImg,           // fulldisk mask at T_REC
        int maskStride,          // stride length in maskImg (NAXIS1)
        char *outQuery)          // e.g., hmi.Mharp[H][T_REC]
{
  int rv;

  if (pInfo->valid == Patch_Invalid) return 0; // not an error
  // 1: Read the bitmap mask of the HARP
  if (pInfo->valid == Patch_Normal) {
    // read fits bitmap
    if (!pInfo->patchName) {
      V_printf(-1, "", "Internal error: Bad patchName during ingest %s.\n", outQuery);
      return 1; // error
    }
    if (load_fits(pInfo->patchName, &(pInfo->image), pInfo->dims) != 0) {
      V_printf(-1, "", "Patch read failed, could not ingest %s.\n", outQuery);
      return 1; // error
    }
    // detect dimension mismatch
    if (pInfo->dims[0] != pInfo->fits_nx || 
    pInfo->dims[1] != pInfo->fits_ny) {
      free(pInfo->image); 
      pInfo->image = NULL;
      V_printf(-1, "", "Patch dims mismatch: fits = (%d,%d) frame = (%d,%d).\n", 
           pInfo->dims[0], pInfo->dims[1],
           pInfo->fits_nx, pInfo->fits_ny);
      V_printf(-1, "", "Patch read failed, could not ingest %s.\n", outQuery);
      return 1; // error
    }
  } else {
    // "padding" images will not have a bitmap
    pInfo->image = NULL;
  }
  // 2: Update HARP outline using the full-disk mask
  // (this implies a size change to pInfo->image, which is reallocated)
  if (update_bitmap(pInfo, maskImg, maskStride) != 0) {
    V_printf(-1, "", "Could not compute full HARP mask at %s\n", outQuery);
    return 1; // error
  }

  // 3: Write the HARP at T_REC
  // Also frees pInfo->image
  //    rv = 0 for success, 1 if some errors, 2 if record not inserted
  rv = ingest_record(pInfo, hInfo, tRec1, tRec, runInfo, 
             mInfo, magRec, maskRec, outQuery);
  if (rv == 2) {
    V_printf(-1, "", "Failed to insert %s\n", outQuery);
    return 1;  // failure
  } else if (rv == 1) {
    // keyword/link errors, not noteworthy
    V_printf(verbflag > 1, "", "Some nonfatal errors with %s\n", outQuery);
  }

  return 0; // OK
}

/* 
 * Update bitmap image so it is correct size, and contains mask overlay
 * Input: 
 *    inner HARP bitmap image, as read from file, in pInfo->image
 *    full-disk mask image
 *    desired bitmap image bounding box in pInfo->{xmin, xmax, ymin, ymax}
 *    the number of pixels in a row of the *full-disk* mask (pix_max, like NAXIS1)
 * Output: 
 *    updates pInfo->image with enlarged version according to bbox
 * 
 * Returns 0 on success, 1 on failure
 */

static 
int
update_bitmap(patch_info_t *pInfo, char *maskImg, int maskStride)
{
  char *patch_map = pInfo->image;  // HARP mask image
  char  patch_val;                 // one value from above
  int nx = pInfo->xmax - pInfo->xmin + 1; // desired HARP dimensions
  int ny = pInfo->ymax - pInfo->ymin + 1;
  int i, j, i0, j0; // counters

  if (pInfo->valid == Patch_Invalid)
    return 0; // not an error
  // compute the new HARP cutout by merging in the mask to the HARP outline
  // note: merging can cause some off-HARP pixels in pInfo->image to be NaN, 
  // which we correct here by using the fulldisk mask image
  // note dims[] is still the old fits file dimensions
  // Old condition:
  // if (pInfo->xmin != pInfo->x0 || 
  //     pInfo->ymin != pInfo->y0 ||
  //     nx != pInfo->dims[0] || 
  //     ny != pInfo->dims[1]) {

  // space for the updated HARP mask image
  pInfo->image = calloc(nx*ny, sizeof(*(pInfo->image)));
  if (pInfo->image == NULL) {
    pInfo->image = patch_map; // may as well keep the old one
    V_printf(-1, "", "Unable to allocate space for new patch bitmap\n");
    return 1;
  }
  // copy maskImg into new cutout first (it is the "background")
  // note: xmin, ymin are 0-based offsets: (0,0) is the top corner
  for (j = 0; j < ny; j++)
    for (i = 0; i < nx; i++)
      pInfo->image[j*nx + i] = maskImg[(j + pInfo->ymin) * maskStride + 
                       (i + pInfo->xmin)];
  V_printf(verbflag > 2, "\t\t    ", "shifting mask (nx,ny) = (%d, %d) -> (%d, %d)\n", 
       pInfo->dims[0], pInfo->dims[1], nx, ny);
  // overlay the HARP bitmap, if it was found, onto the mask
  //   (valid == Patch_Padding -> no HARP bitmap)
  if (patch_map) {
    for (j = 0; j < ny; j++) {
      for (i = 0; i < nx; i++) {
    // compute old bitmap index (i0, j0) from expanded bitmap index (i,j)
    // (keep these lines together for clarity, ignore minor speed penalty)
    // note: xmin,ymin and x0,y0 are 0-based offsets
    i0 = pInfo->xmin + i - pInfo->x0;
    j0 = pInfo->ymin + j - pInfo->y0;
    // ensure (i0, j0) is in the old bitmap
    if (i0 < 0 || i0 >= pInfo->dims[0]) continue;
    if (j0 < 0 || j0 >= pInfo->dims[1]) continue;
    patch_val = patch_map[j0 * (pInfo->dims[0]) + i0];
    // plug in new value only if it is nonzero
    // (the unknown portion caused due to merges is zero)
    if (patch_val != 0)
      pInfo->image[j*nx+i] = patch_val;
      }
    }
    // done with initial HARP overlay
    free(patch_map);
  }
  // change the patch metadata to reflect the new size
  pInfo->x0 = pInfo->xmin; 
  pInfo->y0 = pInfo->ymin;
  pInfo->dims[0] = nx; 
  pInfo->dims[1] = ny;
  return 0;
}

/* 
 * Create and write out record: one HARP at one T_REC
 *
 * Return codes:
 *   0 = all OK
 *   1 = some trouble (keyword, link), but record was inserted
 *   2 = fatal, could not insert this record
 *
 * As a side effect, also frees pInfo->image (if it is non-NULL)
 */

static
int 
ingest_record(patch_info_t *pInfo, 
          harp_info_t  *hInfo,
          trec_info_t  *tRec1,     // this particular T_REC
          trec_info_t  *tRec,      // all T_REC's, for keys
          run_info_t   *runInfo,   // run metadata, for keys
          marp_info_t  *mInfo,     // matching ARs, for keys
          DRMS_Record_t *magRec, 
          DRMS_Record_t *maskRec, 
          char *outQuery)
{
  DRMS_Record_t *outRec = NULL;
  DRMS_Segment_t *outSeg = NULL;
  DRMS_Array_t *outData = NULL;
  DRMS_Link_t *link;
  int datavals, totvals;
  double xDist, yDist;
  char *inst_mode; // hmi or mdi
  int status = 0;  // this routine's error status: 0/1/2
  int rv;          // error code for subroutine we call ("return value")
  int ok, not_ok;  // keyword error count
  int i;
    
  rv = 0; // the generic error code return, not the same as "status"
  outRec = drms_create_record(drms_env, outQuery, DRMS_PERMANENT, &rv);
  if (rv) {
    V_printf(-1, "", "Could not create record for patch bitmap.\n");
    // 2 is not fatal. Free image.
    if (pInfo->image)
    {
        free(pInfo->image);
        pInfo->image = NULL;
    }
    return 2;
  }
  outSeg = drms_segment_lookup(outRec, "bitmap");
  if (!outSeg) {
    V_printf(-1, "", "Could not find bitmap segment in patch.\n");
    drms_close_record(outRec, DRMS_FREE_RECORD);
    // 2 is not fatal. Free image.
    if (pInfo->image)
    {
        free(pInfo->image);
        pInfo->image = NULL;
    }
    return 2;
  }
  // pInfo->image will be freed when the DRMS array is freed, below
  outData = drms_array_create(DRMS_TYPE_CHAR, 2, pInfo->dims, pInfo->image, &rv);
  if (rv) {
    V_printf(-1, "", "Could not create array for patch bitmap segment.\n");
    if (outData)
    {
        drms_free_array(outData); // Frees pInfo->image.
        pInfo->image = NULL; // drms_free_array() doesn't set the pointer to NULL.
    }
    return 2; // fatal
  }
  // link outData into outSeg
  outSeg->axis[0] = outData->axis[0];
  outSeg->axis[1] = outData->axis[1];
  outData->parent_segment = outSeg;
    
  /*
   * Keys
   */
  inst_mode = get_instrument_mode(magRec); // use HMI vs. MDI keywords
  // prime keys
  not_ok = 0;  // count failed keys
  ok = drms_setkey_int   (outRec, "HARPNUM", pInfo->num);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_time  (outRec, "T_REC", tRec1->t);
  if (ok != DRMS_SUCCESS) not_ok++;
  // JSOC information
  ok = drms_setkey_string(outRec, "BLD_VERS", jsoc_version);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_time  (outRec, "DATE", CURRENT_SYSTEM_TIME);
  if (ok != DRMS_SUCCESS) not_ok++;
  // code version
  not_ok += set_keys_code_info(inst_mode, outRec);

  // per tracker run keys
  not_ok += set_keys_runinfo(outRec, runInfo);
  // patch summary statistics
  not_ok += set_keys_stats(inst_mode, outRec, magRec, maskRec, pInfo, hInfo, tRec);
  // match summary
  not_ok += set_keys_match(outRec, hInfo, mInfo);
  // other keys
  totvals = (pInfo->dims[0])*(pInfo->dims[1]);
  for (datavals = i = 0; i < totvals; i++) 
    if (pInfo->image[i] != 0) datavals++;  // TODO: not enough info to check if on-disk
  ok = drms_setkey_int(outRec, "TOTVALS",  totvals);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_int(outRec, "DATAVALS", datavals);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_int(outRec, "MISSVALS", totvals-datavals);
  if (ok != DRMS_SUCCESS) not_ok++;
  // use WCS in tRec1 to get new wcs
  distance2center(tRec1->wcs, pInfo->x0, pInfo->y0, &xDist, &yDist);
  ok = drms_setkey_float(outRec, "CRVAL1", xDist);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_float(outRec, "CRVAL2", yDist);
  if (ok != DRMS_SUCCESS) not_ok++;
  // summary
  if (not_ok > 0) {
    V_printf(-1, "\t\t", "Failed to set up %d patch keys: %s\n", not_ok, tRec1->str);
    // commented out to make a simple key insertion failure not notable
    // status = 1; 
  }

  // Links
  rv = drms_setlink_static(outRec, "MDATA", magRec->recnum);
  if (rv) {
    status = 1; // is a nonfatal error
    V_printf(-1, "\t", "MDATA link failed\n");
  }
  rv = drms_setlink_static(outRec, "MASK", maskRec->recnum);
  if (rv) {
    status = 1; // is a nonfatal error
    V_printf(-1, "\t", "MASK link failed\n");
  }
  // Write the segment
  rv = drms_segment_write(outSeg, outData, 0);
  if (rv) {
    V_printf(-1, "", "Failed to write segment for patch.\n");
    drms_free_array(outData);
    return 2;  // fatal
  }
  drms_free_array(outData);  // frees pInfo->image
  pInfo->image = NULL; // drms_free_array() doesn't set the pointer to NULL.
  drms_close_record(outRec, DRMS_INSERT_RECORD);
  return status;  // OK
}


/*
 * utility: return 1 if tag is invalid, 0 otherwise
 * (allow Faint = -1 here, even though don't expect to read it)
 */
static
int 
tag_is_invalid(int tag)
{
  if ((tag == (int) Patch_Tag_Faint ) ||
      (tag == (int) Patch_Tag_Harder) ||
      (tag == (int) Patch_Tag_Normal) ||
      (tag == (int) Patch_Tag_Merge ))
    return 0;  // valid
  else
    return 1;  // invalid
}

/* 
 * Coordinate-Frame file parser, one line per call
 * return 0 if successful, -1 for EOF, +1 for format error
 */
static
int 
load_frame(FILE *fp, patch_info_t *pInfo, TIME *t)
{
  const int rv_expect = 12; // tied to format string below
  char t_rec[24]; // time
  int tag;        // tag
  int rv;

  pInfo->valid = Patch_Invalid; // set to invalid
  // reading rv_expect items
  rv = fscanf(fp, "%d %23s %d %d %d %d %d %f %f %f %f %f", 
          &(pInfo->num), 
          t_rec, // string
          &tag,
          &(pInfo->fits_nx), 
          &(pInfo->fits_ny),
          &(pInfo->x0), 
          &(pInfo->y0),
          &(pInfo->lat0), 
          &(pInfo->lon0), 
          &(pInfo->lat1), 
          &(pInfo->lon1),
          &(pInfo->omega));
  if (rv == EOF) return -1; // regular EOF
  if (rv != rv_expect) {
    V_printf(-1, "", "Warning: could not read frame line, bad file?\n");
    return 1;
  }
  if (tag_is_invalid(tag)) {
    V_printf(-1, "", "Warning: could not read frame tag (%d), bad value?\n", tag);
    return 1;
  } else {
    pInfo->tag = (patch_tag_t) tag;
  }
  // convert inputs, which were from Matlab and have origin=1 (x0,y0), 
  // to origin=0 in keeping with typical C usage
  pInfo->x0 -= 1;
  pInfo->y0 -= 1;
  // convert the time to a rec number
  *t = sscan_time(t_rec);
  if (time_is_invalid(*t)) {
    V_printf(-1, "", "Warning: could not read frame time (%s), bad file?\n", t_rec);
    return 1;
  }
  if (isnan(pInfo->lat0) || isnan(pInfo->lat1) ||
      isnan(pInfo->lon0) || isnan(pInfo->lon1) || isnan(pInfo->omega)) {
    V_printf(-1, "", "Warning: found a NaN in frame line, bad file?\n");
    return 1;
  }

  // now it's OK
  pInfo->valid = Patch_Normal; // regular patch
  return 0;
}

/* 
 * Stats file parser
 * return 0 if successful, -1 for EOF, +1 for format error
 */
static
int 
load_stats(FILE *fp, float *stats)
{
  int i;
  int rv;
  float x;
    
  for (i = 0; i < RS_num_stats; i++) {
    rv = fscanf(fp, "%f", &x); // EOF, 1, or 0
    if (rv == EOF) return -1; // regular end-of-file
    if (rv != 1) {
      V_printf(-1, "", "Warning: load_stats could not read stat #%d, bad file?\n", i+1);
      return 1; // error
    }
    stats[i] = x;
    // printf("%g ", stats[i]);
  }
  return 0; // OK
}

/*
 * Fits reader
 * returns 0 for success, 1 for error
 */

static
int 
load_fits(char *filename, char **image, int *dims)
{
  char *data = NULL;     // must init to null
  fitsfile *fptr = NULL; // fitsio input file pointer
  int fstatus = 0;       // fitsio status variable: all ok now
  int mem_err = 0;       // no allocation error now
  int bitpix;            // Data type of fits file
  int naxis;             // number of axes from fits file
  long datasize;         // number of pixels
  long dims_long[2];     // dimension of image
  long fpixel[2] = {1, 1}; // start reading at image pixel (1,1)
  
  fits_open_file(&fptr, filename, READONLY, &fstatus);
  if (fstatus != 0) goto done;
  fits_get_img_param(fptr, 2, &bitpix, &naxis, dims_long, &fstatus);
  if (fstatus != 0) goto done;
  // sizing
  datasize = dims_long[0] * dims_long[1];
  if ((data = calloc(datasize, sizeof(*data))) == NULL) { 
    mem_err = 1; // allocation error, not a fits error
    goto done;
  }
  // load image array
  fits_read_pix(fptr, TSBYTE, fpixel, datasize, NULL, data, NULL, &fstatus);
  if (fstatus != 0) goto done;
  // note, fstatus = 0
 done:
  if (fptr) 
    fits_close_file(fptr, &fstatus);
  if (fstatus || mem_err) {
    // report the error
    char err_msg[FLEN_STATUS];
    if (mem_err)
      strncpy(err_msg, "failed calloc for image data", sizeof(err_msg));
    else
      fits_get_errstatus(fstatus, err_msg);
    V_printf(-1, "", "Error opening fits file `%s': %s\n", filename, err_msg);
    // free the image if it was allocated (data == NULL is OK)
    free(data);
    // clear outputs
    *image = NULL;
    dims[0] = dims[1] = 0;
  } else {
    // set up outputs
    *image = data;
    dims[0] = dims_long[0]; 
    dims[1] = dims_long[1];
  }
  return (fstatus || mem_err) ? 1 : 0;
}

/*
 * what flavor of keywords do we use 
 *    determined using `TELESCOP' KW
 *    currently just HMI vs. MDI 
 *    (note, INSTRUME is not present in synoptic MDI mags)
 */
static
char *
get_instrument_mode(DRMS_Record_t *xmRec)
{
  char * const default_mode = "HMI";
  int status = 0;
  char *kw;

  kw = drms_getkey_string(xmRec, "TELESCOP", &status);
  if (status == 0) {
    // get the instrument, and normalize it to a standard set
    if (strcasestr(kw, "SDO") || strcasestr(kw, "HMI"))
      return "HMI";
    else if (strcasestr(kw, "SOHO") || strcasestr(kw, "MDI"))
      return "MDI";
    else
      return default_mode;
  } 
  return default_mode;
}

/*
 * set up code information keywords 
 */
static int
set_keys_code_info(char *inst_mode, DRMS_Record_t *outRec)
{
  extern char *hmi_mharp_version; // defined globally
  char *doc_url;
  char keystr[80]; // plausible FITS keyword length
  int ok;
  int not_ok = 0;

  // mode switch for documentation
  //   url points to descriptive text in wiki
  if (strcmp(inst_mode, "HMI") == 0)
    doc_url = "http://jsoc.stanford.edu/jsocwiki/HARPDataSeries";
  else if (strcmp(inst_mode, "MDI") == 0)
    doc_url = "http://jsoc.stanford.edu/jsocwiki/MTARPDataSeries";
  else
    doc_url = "No known documentation";
  ok = drms_setkey_string(outRec, "HRPDOCU", doc_url);
  if (ok != DRMS_SUCCESS) not_ok++;
  // version number (of the present code)
  snprintf(keystr, sizeof(keystr), "%s ver: %s", module_name, hmi_mharp_version);
  ok = drms_setkey_string(outRec, "HRPCODEV", keystr);
  if (ok != DRMS_SUCCESS) not_ok++;
  return not_ok;
}

/*
 * set_keys_runinfo: set up drms keys using per-tracker-run keys
 * which are stored in runInfo
 */
static
int
set_keys_runinfo(DRMS_Record_t *rec, 
         run_info_t *runInfo)
{
  const char **keys = tracker_run_info_keys;
  const char *drms_key;
  const char *tracker_key;
  int not_ok, iKey;    // outer loop
  int p, found, ok;    // inner loop

  for (not_ok = iKey = 0; keys[iKey] != NULL; iKey += 2) {
    // map drms_key -> tracker_key
    drms_key    = keys[iKey];
    tracker_key = keys[iKey+1];
    // find a parameter p matching tracker_key
    for (found = p = 0; p < runInfo->n_par; p++) {
      if (strcasecmp(tracker_key, runInfo->par_names[p]) == 0) {
    found = 1;
    ok = drms_setkey_string(rec, drms_key, runInfo->par_vals[p]);
    if (ok != DRMS_SUCCESS) not_ok++;
      }
    }
    if (!found) not_ok++; // could not find a tracker param for the key
  }
  return not_ok;
}

/* 
 * Set HARP keywords for region-matches (NOAA_AR, etc.)
 * returns: number of failed insertions (0 if all OK)
 */
static
int
set_keys_match(DRMS_Record_t *rec, 
           harp_info_t  *hInfo, 
           marp_info_t  *mInfo)   // NB: head of *all* marpInfo array
{
  int ok, status;
  int not_ok = 0; // count failed insertions
  const int nMatch = hInfo->n_match;
  int best_match, m;
  char str[20];  // enough for an int in decimal
  char all_match[STR_MAX];

  // number of matching ARs (0 OK)
  ok = drms_setkey_int(rec, "NOAA_NUM", nMatch);
  if (ok != DRMS_SUCCESS) not_ok++;
  // NOAA ID of best matching AR (0 OK)
  if (nMatch == 0)
    best_match = 0; // no match
  else
    best_match = mInfo[hInfo->top_match].id;
  ok = drms_setkey_int(rec, "NOAA_AR", best_match);
  if (ok != DRMS_SUCCESS) not_ok++;
  // string of all matches (more than one OK)
  *all_match = '\0';
  for (m = 0; m < nMatch; m++) {
    snprintf(str, sizeof(str), "%s%d", (m == 0) ? "" : ",", mInfo[hInfo->match[m]].id);
    strncat(all_match, str, sizeof(all_match));
  }
  ok = drms_setkey_string(rec, "NOAA_ARS", all_match);
  if (ok != DRMS_SUCCESS) not_ok++;
  return not_ok;
}

/* 
 * Set HARP keywords
 *
 * returns: number of failed insertions (0 if all OK)
 */
static
int
set_keys_stats(char *inst_mode,
           DRMS_Record_t *rec, 
           DRMS_Record_t *magRec, 
           DRMS_Record_t *maskRec, 
           patch_info_t *pInfo, 
           harp_info_t  *hInfo, 
           trec_info_t  *tRec)   // NB: head of *all* trec array
{
  float *stats = pInfo->stats;
  int ok, status;
  int not_ok = 0; // count failed insertions

  status = propagate_keys_harp(inst_mode, magRec, maskRec, rec);
  if (status == -1) {
    // this should never happen
    V_printf(-1, "", "set_keys: propagate_keys failed (null record)\n");
  } else {
    not_ok += status;
  }
  // remainder of WCS
  // x0, y0 are 0-based offsets
  ok = drms_setkey_float(rec, "CRPIX1", pInfo->x0);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_float(rec, "CRPIX2", pInfo->y0);
  if (ok != DRMS_SUCCESS) not_ok++;
  // sizes (not part of WCS, but related to CRPIX* above)
  ok = drms_setkey_float(rec, "CRSIZE1", pInfo->dims[0]);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_float(rec, "CRSIZE2", pInfo->dims[1]);
  if (ok != DRMS_SUCCESS) not_ok++;
  // harp-extent metadata
  ok = drms_setkey_time(rec, "T_FRST",  tRec[hInfo->rec0p].t);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_time(rec, "T_FRST1", tRec[hInfo->rec0 ].t);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_time(rec, "T_LAST1", tRec[hInfo->rec1 ].t);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_time(rec, "T_LAST",  tRec[hInfo->rec1p].t);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_int(rec, "N_PATCH",  hInfo->n_patchp);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_int(rec, "N_PATCH1", hInfo->n_patch);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_int(rec, "N_PATCHM", hInfo->n_missing);
  if (ok != DRMS_SUCCESS) not_ok++;
  // other metadata, from frame file
  // tag, from ROI_s(patch).tag
  ok = drms_setkey_int(rec, "H_MERGE", (pInfo->tag == Patch_Tag_Merge ) ? 1 : 0);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_int(rec, "H_FAINT", 
               ((pInfo->tag == Patch_Tag_Harder) || 
            (pInfo->tag == Patch_Tag_Faint )) ? 1 : 0);
  if (ok != DRMS_SUCCESS) not_ok++;
  // Note, these are for the whole disk transit ("DT")
  ok = drms_setkey_float(rec, "LONDTMIN", pInfo->lon0);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_float(rec, "LATDTMIN", pInfo->lat0);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_float(rec, "LONDTMAX", pInfo->lon1);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_float(rec, "LATDTMAX", pInfo->lat1);
  if (ok != DRMS_SUCCESS) not_ok++;
  ok = drms_setkey_float(rec, "OMEGA_DT", pInfo->omega);
  if (ok != DRMS_SUCCESS) not_ok++;
  // insert "stats" keywords
  int sn;
  for (sn = 0; sn < RS_num_stats; sn++) {
    char *name1 = (char *) RS_index2keyname[sn];
    //      printf("%s ", name1+1);
    // insert if HMI keyword name is non-null
    if (name1) {
      switch (*name1) {
      case 'i':
    ok = drms_setkey_int(rec, name1+1, (int) stats[sn]);
    break;
      case 'f':
    ok = drms_setkey_float(rec, name1+1, (float) stats[sn]);
    break;
      default:
    // should never happen
    V_printf(-1, "", "set_keys: Unknown key type (first char of keyname %s)\n", name1);
    ok = DRMS_SUCCESS + 1; // i.e., != DRMS_SUCCESS
    break;
      }
      if (ok != DRMS_SUCCESS) not_ok++;
      V_printf(ok != DRMS_SUCCESS, "", "set_keys: failed to set %s -> %f\n", name1, stats[sn]);
    }
  }
  return not_ok;
}


/* 
 * Find bounding box from given lat/lon range 
 *
 * Computes: pInfo->{xmin, xmax, ymin, ymax}
 *
 * return value:
 *   1 if all pixels were off-disk (not an error)
 *   0 otherwise
 *
 * Algorithm: Scan along the edge, sampling in the lat/lon space
 *   with given resolution, nominally ~ 0.03 degree
 * Not sure if the caller cares about +1 value, but someone should 
 * handle this case.
 * We handle the case where the hi/lo pixel site is not on the 
 * lat/lon bounding box (e.g., lat1 = -lat0 = 20, lon0 = 45,
 * lon1 = 45+90).
 */
// macro to update local variables
#define UPDATE_RANGE  do { \
                       ondisk = 1; \
                       if (x < xmin) xmin = x; \
                       if (x > xmax) xmax = x; \
                       if (y < ymin) ymin = y; \
                       if (y > ymax) ymax = y; \
                      } while (0)

static
int 
compute_boundary(patch_info_t *pInfo, int pix_max, wcs_t wcs)
{
  int farside, ondisk;
  double xmin, xmax, ymin, ymax;
  double x, y, lat, lon, lat0, lon0, lat1, lon1, dl;        // loop control
  /*
   *  ephemeris
   */
  // standard local variables for WCS
  double crvalx, crvaly, crpix1, crpix2, cdelt, crota2, sina, cosa;
  WCS2LOCALS(wcs); // initialize the above vars
  // center of disk, starting at 0 -- these macros use the above variables
  double xcen   = PIX_X(0.0,0.0) - 1.0;
  double ycen   = PIX_Y(0.0,0.0) - 1.0;
  // extras for img2sphere()/sphere2img()
  double latc   = wcs.crlt_obs * RADSINDEG;
  double lonc   = wcs.crln_obs * RADSINDEG;
  double peff   = - wcs.crota2 * RADSINDEG;
  double ang_r  = asin(wcs.rsun_ref / wcs.dsun_obs); // rSun in radians -- about .26 deg
  double rsun   = asin(wcs.rsun_ref / wcs.dsun_obs) * RAD2ARCSEC / wcs.cdelt1;
    
  V_printf(verbflag > 2, "\t\t  ", "cen = (%.2f, %.2f) R=%f latc=%.2f lonc=%.2f\n", 
       pInfo->num, xcen, ycen, rsun, latc/RADSINDEG, lonc/RADSINDEG);

  // Step 0: The initial lat/lon bounding box in degrees
  lat0 = pInfo->lat0 * RADSINDEG;
  lat1 = pInfo->lat1 * RADSINDEG;
  lon0 = pInfo->lon0 * RADSINDEG + lonc;
  lon1 = pInfo->lon1 * RADSINDEG + lonc;
  dl = 0.03 * RADSINDEG; // lat and lon step size -- a design parameter
  
  // Step 1: Traverse lat/lon BB edges to find pixel BB
  /*  sphere2img synopsis:
   *  Perform a mapping from heliographic coordinates latitude and longitude
   *    (in radians) to plate location on an image of the sun.  The plate
   *    location is in units of the image radius and is given relative to
   *    the image center.  The function returns 1 if the requested point is
   *    on the far side (>90 deg from disc center), 0 otherwise.
   */
  farside = 0; // do any pixels in the lat/lon BB map to the far side?
  // the five variables below are maintained by the UPDATE_RANGE macro
  ondisk = 0;  // do any pixels in the lat/lon BB map to the visible disk?
  xmin = ymin =  HUGE_VAL; // will be reset in the loop
  xmax = ymax = -HUGE_VAL;
  // Scan latitude
  for (lat = lat0; lat <= lat1; lat += dl) {
    if (!sphere2img(lat, lon0, latc, lonc, &x, &y, xcen, ycen,
            rsun, peff, 0., 0., 0., 0.))
      UPDATE_RANGE; 
    else 
      farside = 1;
    if (!sphere2img(lat, lon1, latc, lonc, &x, &y, xcen, ycen,
            rsun, peff, 0., 0., 0., 0.))
      UPDATE_RANGE; 
    else 
      farside = 1;
  }
  // Scan longitude
  for (lon = lon0; lon <= lon1; lon += dl) {
    if (!sphere2img(lat0, lon, latc, lonc, &x, &y, xcen, ycen,
            rsun, peff, 0., 0., 0., 0.))
      UPDATE_RANGE;
    else 
      farside = 1;
    if (!sphere2img(lat1, lon, latc, lonc, &x, &y, xcen, ycen,
            rsun, peff, 0., 0., 0., 0.))
      UPDATE_RANGE;
    else 
      farside = 1;
  }
  V_printf(verbflag > 2, "\t\t  ", 
       "BB     : (%.3f, %.3f) X (%.3f, %.3f), %s, %s\n", xmin, xmax, ymin, ymax, 
       ondisk  ? "some on-disk pixels" : "no visible pixels", 
       farside ? "some farside pixels" : "no farside pixels");
  // off-disk, skip it
  if (!ondisk)
    return 1; // whole HARP, not just one pixel, was on farside
  // below this if (!ondisk), the BB has been set up OK

  // Step 2: If the box crossed partly to the farside, expand all sides.
  // We may have missed the precise boundary as we traversed the edge.
  if (farside) {
    xmin -= 1.0;
    xmax += 1.0;
    ymin -= 1.0;
    ymax += 1.0;
  }

  // Step 3: Check if the equator/poles were in the lat/lon BB
  // If they were, the pixel BB must be adjusted because the extremum might
  // not be on any edge of the lat/lon BB.
  // int img2sphere (double x, double y, 
  //                 double ang_r, double latc, double lonc, double pa, 
  //                 double *rho, double *lat, double *lon, 
  //                 double *sinlat, double *coslat, 
  //                 double *sig, double *mu, double *chi)
  double rho, sinlat, coslat, sig, mu, chi;  // dummy output arguments
  int rv;
  for (double xpos = -1; xpos <= 1; xpos += 1)
    for (double ypos = -1; ypos <= 1; ypos += 1) {
      if (!(xpos == 0 || ypos == 0)) continue;
      if (xpos == 0 && ypos == 0) continue;
      rv = img2sphere(xpos, ypos, 
              ang_r, latc, lonc, peff,
              &rho, &lat, &lon, &sinlat, &coslat, &sig, &mu, &chi);
      if (rv)
    V_printf(verbflag > 3, "\t\t\t", "(%f,%f) not on sphere\n", xpos, ypos);
      else
    V_printf(verbflag > 3, "\t\t\t", 
         "(x,y) = (% .1f,% .1f) -> (lat,lon) = (% 6.1f,% 6.1f)\n", 
         xpos, ypos, lat/RADSINDEG, lon/RADSINDEG);
      // Determine whether (lat,lon) is in the lat-lon BB
      while (lon < lon0) lon += 2*M_PI;  /* wrap lon to best fit in [lon0,lon1] */
      while (lon > lon1) lon -= 2*M_PI;
      if (lat0 <= lat && lat <= lat1 && lon0 <= lon && lon <= lon1) {
    // this is rare, so announce it
    V_printf(verbflag > 2, "\t\t  ", "Note: Extremum found off the box boundary.\n");
    // patch up the pixel BB to include the pixel:
    //  (xpos, ypos) * rsun + (xc,yc)
    x = xcen + rsun*xpos;
    y = ycen + rsun*ypos;
    UPDATE_RANGE;
      }
    }

  // Step 4: Make pixel BB into an integer box and clip to image size
  // clip to legal range {0..Npix - 1}
  V_printf(verbflag > 2, "\t\t  ", "BB pre : (%.3f, %.3f) X (%.3f, %.3f)\n", 
       xmin, xmax, ymin, ymax);
  if (xmin <  0   )    xmin = 0;
  if (xmax >= pix_max) xmax = pix_max - 1;
  if (ymin <  0   )    ymin = 0;
  if (ymax >= pix_max) ymax = pix_max - 1;
  V_printf(verbflag > 2, "\t\t  ", "BB clip: (%.3f, %.3f) X (%.3f, %.3f)\n", 
       xmin, xmax, ymin, ymax);
  // integer boundaries
  pInfo->xmin = floor(xmin); 
  pInfo->xmax = ceil(xmax);
  pInfo->ymin = floor(ymin); 
  pInfo->ymax = ceil(ymax);
  V_printf(verbflag > 2, "\t\t  ", "BB end : (%d, %d) X (%d, %d)\n", 
       pInfo->xmin, pInfo->xmax, pInfo->ymin, pInfo->ymax);
  
  return 0;
}

// clean namespace
#undef UPDATE_RANGE


/*
 * get distance of points x and y from center to set up WCS
 *   x, y start from 0
 */
static
void
distance2center(wcs_t wcs, double x, double y, double *xDist, double *yDist)
{
  // standard local variables for WCS
  double crvalx, crvaly, crpix1, crpix2, cdelt, crota2, sina, cosa;
  WCS2LOCALS(wcs); // initialize the above vars from wcs
  
  // center location
  double pix_x = x + 1.0;
  double pix_y = y + 1.0;
    
  // standard macros rely upon some of above "WCS2LOCALS" vars
  *xDist = WX(pix_x,pix_y);
  *yDist = WY(pix_x,pix_y);
    
  // printf("%f,%f\n", xDist, yDist);
}

---