proj/globalhs/apps/undistortmdi.c

/* this module is based on jv2helio, but instead interpolates its input
 * to the desired resolution, optionally correcting for distortion and 
 * errors in the p-angle and carrington inclination.
 */


#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <mkl_lapack.h>  //needed for dsecnd()
#include "jsoc_main.h"
#include "fstats.h"
#include "drms_dsdsapi.h"
#include "projection.h"

#define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
#define QUAL_NODATA (0x80000000)

#define PI      (M_PI)
#define RADSINDEG   (PI/180)
#define ARCSECSINRAD    (3600*180/PI)
#define DAYSINYEAR  (365.2425)
#define SECSINDAY   (86400)
#define TAU_A           (499.004783806) // light travel time in seconds, = 1 AU/c
//#define TAU_A     (499.004782)    // this value used in old v2helio
#define TCARR       (25.38)     // days
#define RTRUE       (6.96000000e8)  // meters
#define AU      (149597870691)  // meters/au
//#define AU        (1.49597870e11) // this value used in old v2helio
#define MAXLEN      (256)
#define NO_DATASET  (-1)
#define NO_IMAGE    (-1)
#define kMAX_SKIPERRMSG  1024
#define kMAXROWS        65536

#define kNOTSPECIFIED   "not specified"

/* Must do decls */
double dsecnd();

typedef enum
{
   V2HStatus_Success,
   V2HStatus_MissingParameter,
   V2HStatus_IllegalTimeFormat,
   V2HStatus_TimeConvFailed,
   V2HStatus_Unimplemented,
   V2HStatus_IllegalOrientation
} V2HStatus_t;

char *module_name = "undistortmdi";
char *cvsinfo_undistortmdi = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/undistortmdi.c,v 1.4 2013/05/03 21:05:08 tplarson Exp $";

ModuleArgs_t module_args[] = 
{
   {ARG_STRING,  "in", "", "input data records"},
   {ARG_STRING,  "out", "", "output data series"},
   {ARG_STRING,  "segin", kNOTSPECIFIED, "name of input segment if not using segment 0"},
   {ARG_STRING,  "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
   {ARG_STRING,  "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
   {ARG_STRING,  "srclink",  "SRCDATA", "name of link to source data"},
   {ARG_INT,     "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
   {ARG_INT,     "VERB", "1", "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", ""},  /* debug messages */
   {ARG_INT,     "OUTCOLS", "4096", "number of output columns"},
   {ARG_INT,     "OUTROWS", "4096", "number of output rows"},
   {ARG_FLOAT,   "MAPRMAX",  "1.1", "maximum image radius", ""},
   {ARG_INT,     "RMAXFLAG", "1", "set to nonzero to set output to DRMS_MISSING outside MAPRMAX, otherwise uses 0.0 outside MAPRMAX"},  /* Non 0 sets data outside RMAX MISSING */   
   {ARG_INT,     "INTERPO",  "1", "option for interpolation: 0=bilinear; 1=cubic convolution", ""}, /* 2 methods - see soi_fun.h */
   {ARG_FLOAT,   "OUTSCALE", "1.0", "bscale to use for output", ""},    /* scale for output */
   {ARG_FLOAT,   "OUTBIAS",  "0.0", "bzero to use for output", ""},    /* bias for scaled output */
   {ARG_INT,     "DISTORT",  "0", "option for distortion correction: 0=none; 1=full disk(fd) data; 2=vector-weighted(vw) data", ""}, /* 0 for none, 1 for FD, 2 for vw */
   {ARG_FLOAT,   "CUBIC",    "7.06E-9", "cubic distortion in fd units", ""}, /* Cubic distortion in FD units */
   {ARG_FLOAT,   "TILTALPHA","2.59", "tilt of CCD, degrees", ""}, /* TILT of CCD in degrees */
   {ARG_FLOAT,   "TILTBETA", "56.0", "direction of CCD tilt, degrees", ""}, /* Direction of TILT in degrees */
   {ARG_FLOAT,   "TILTFEFF", "12972.629", "effective focal length", ""}, /* Effective focal length */
   {ARG_INT,     "DATASIGN",  "0", "value to multiply data; set to 0 to take DATASIGN from keyword, or 1.0 if not found", ""},  /* Non 0 forces datasign to value*/
   {ARG_INT,     "MAXMISSVALS", "0", "if >0, this becomes threshold on MISSVALS from keyword", ""},  /* max. allowed MISSING pixels */
   {ARG_FLOAT,   "SOLAR_P",  "999.0", "P-angle; if unset, taken from keywords", ""},    /* can't use D_MISSING here */
   {ARG_FLOAT,   "PSIGN",    "1.0", "sign of SOLAR_P", ""}, /* Sign of P. For MWO data. */
   {ARG_FLOAT,   "PERR",     "0.0", "fixed P-angle error, likely -0.22", ""},   /* Fixed P-angle error. Maybe -0.22. */
   {ARG_FLOAT,   "IERR",     "0.0", "error in Carrington inclination, likely -0.10", ""},   /* Error in Carrington inclination. Maybe -0.10. */
   {ARG_TIME,    "REF_TB0",  "2001.06.06_06:57:22_TAI", "reference time for computing correction to P and B angles, roughly when B0=0", ""},
   {ARG_END}
};

#include "imageinterp.c"
#include "saveparm.c"
#include "timing.c"
#include "set_history.c"

static inline void ParameterDef(int status, 
                char *pname, 
                double defaultp, 
                double *p, 
                int iRec,
                                int verbflag)
{
   /* logs warning and sets parameter to default value */
   if (status != 0)
   {
      *p = defaultp;
      if (verbflag)
         fprintf(stderr, "WARNING: default value %g used for %s, iRec = %d, status = %d\n", defaultp, pname, iRec, status);
   }
}


#define PARAMETER_ERROR(PNAME)                                     \
  if (status != DRMS_SUCCESS)                                      \
  {                                                                \
    CreateBlankRecord(inrec, outrec, quality);                     \
    fprintf(stderr,                                                \
        "SKIP: error getting keyword %s: iRec = %d, status = %d, T_REC = %s, recnum = %lld\n",  \
        PNAME,                                                 \
            iRec,                                                  \
        status,                                                \
            trecstr,                                               \
            inrec->recnum);                                        \
    if (inarr)  drms_free_array(inarr);                            \
    if (outarr) drms_free_array(outarr);                           \
    if (orientation) free(orientation);                            \
    continue;                                                      \
  }


/* Segment will be empty, but there will be a record! */
static void CreateBlankRecord(DRMS_Record_t *inrec, DRMS_Record_t *outrec, int quality)
{
  /* create 'blank' data */
// might insert 'quality = quality | MASK' here.
   quality = quality | QUAL_NODATA;
   drms_copykey(outrec, inrec, "T_REC");
   drms_setkey_int(outrec, "QUALITY", quality);
   drms_close_record(outrec, DRMS_INSERT_RECORD);
}


int DoIt(void)
{

   int newstat = 0;
   int status = DRMS_SUCCESS;
   V2HStatus_t vstat = V2HStatus_Success;

   const char *orientationdef = "SESW    ";
   char *orientation = NULL;
   int paramsign;
   int longitude_shift, velocity_correction, interpolation, apodization;
   int mag_correction;
   int mag_offset;
   int row;
   int mapped_lmax, sinb_divisions, mapcols, maprows, nmax, nmin;
   int length[2];
   int carrStretch = 0;
   float diffrotA = 0.0;
   float diffrotB = 0.0;
   float diffrotC = 0.0;
   double tobs, tmearth, tref, trefb0;
   double smajor, sminor, sangle;
   double xscale, yscale, imagescale;
   int xpixels, ypixels, pixels;
   double obs_vr, obs_vw, obs_vn;
   double b0, bmax, bmin, sinBdelta;
   double longmax, longmin, longmax_adjusted, longmin_adjusted, longinterval;
   double p0, p, rmax;
   double ierr, perr, psign;
   double x0, y0;
   double obsdist, longshift, obsl0, refl0, mapl0, longrate, rtrue, rsun, S;
   double rsunDef;
   int obsCR;
   int apel;
   double apinner, apwidth, apx, apy;
   double scale, bias;
   double colsperdeg;

   int quality, NaN_beyond_rmax, nRecs;
   double satrot, instrot;
   double dsignout, vsign;
   int distsave;
   double cubsave, tiltasave, tiltbsave, tiltfsave;
   TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
   char trecstr[100];

   LIBPROJECTION_Dist_t distP;
   DRMS_Segment_t *segin = NULL;
   DRMS_Segment_t *segout = NULL;
   DRMS_Array_t *inarr = NULL;
   DRMS_Array_t *outarr = NULL;
   DRMS_Record_t *inrec = NULL;
   DRMS_Record_t *outrec = NULL;
   DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
   DRMS_RecLifetime_t lifetime;

   long long histrecnum=-1;

   int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
   int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);

   struct timeval tv0, tv1, tv;
   double ct0, ct1, ct;
   double wt0, wt1, wt;
   double ut0, ut1, ut;
   double st0, st1, st;
   double tt0, tt1, tt;

   wt0=getwalltime();
   getcputime(&ut0, &st0);

   gettimeofday(&tv0, NULL);
   ct0=dsecnd();
   tt0=1000*((double)clock())/CLOCKS_PER_SEC;

   char *inRecQuery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
   char *outSeries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);

   int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
   int maxmissvals = cmdparams_save_int(&cmdparams, "MAXMISSVALS", &newstat);
   int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
   if (permflag)
      lifetime = DRMS_PERMANENT;
   else
      lifetime = DRMS_TRANSIENT;


   interpolation = cmdparams_save_int(&cmdparams, "INTERPO", &newstat);
   paramsign = cmdparams_save_int(&cmdparams, "DATASIGN", &newstat);

   distsave = cmdparams_save_int(&cmdparams, "DISTORT", &newstat);
   cubsave = cmdparams_save_double(&cmdparams, "CUBIC", &newstat);
   tiltasave = cmdparams_save_double(&cmdparams, "TILTALPHA", &newstat);
   tiltbsave = cmdparams_save_double(&cmdparams, "TILTBETA", &newstat);
   tiltfsave = cmdparams_save_double(&cmdparams, "TILTFEFF", &newstat);

   scale = cmdparams_save_double(&cmdparams, "OUTSCALE", &newstat);
   bias = cmdparams_save_double(&cmdparams, "OUTBIAS", &newstat);
   p0 = cmdparams_save_double(&cmdparams, "SOLAR_P", &newstat);
   psign = cmdparams_save_double(&cmdparams, "PSIGN", &newstat);
   perr = cmdparams_save_double(&cmdparams, "PERR", &newstat);
   ierr = cmdparams_save_double(&cmdparams, "IERR", &newstat);
   trefb0 = cmdparams_save_time(&cmdparams, "REF_TB0", &newstat);
   rmax = cmdparams_save_double(&cmdparams, "MAPRMAX", &newstat);
   NaN_beyond_rmax = cmdparams_save_int(&cmdparams, "RMAXFLAG", &newstat);

   SetDistort(distsave, cubsave, tiltasave, tiltbsave, tiltfsave, &distP);

   mapcols = cmdparams_save_int(&cmdparams, "OUTCOLS", &newstat);
   maprows = cmdparams_save_int(&cmdparams, "OUTROWS", &newstat);

   length[0] = mapcols;
   length[1] = maprows;

   char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
   char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
   int seginflag = strcmp(kNOTSPECIFIED, segnamein);
   int segoutflag = strcmp(kNOTSPECIFIED, segnameout);

   char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
   char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);

   if (newstat) 
   {
     fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
     cpsave_decode_error(newstat);
     return 1;
   }  
   else if (savestrlen != strlen(savestr)) 
   {
     fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
     return 1;
   }

   // set up ancillary dataseries for processing metadata
   DRMS_Record_t *tempoutrec = drms_create_record(drms_env, 
                           outSeries,
                           DRMS_TRANSIENT, 
                           &status);

   if (status != DRMS_SUCCESS) 
   {
     fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outSeries, status);
     return 1;
   }

   DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
   DRMS_Link_t *srclink = hcon_lookup_lower(&tempoutrec->links, srclinkname);

   if (histlink != NULL) 
   {
     histrecnum=set_history(histlink);
     if (histrecnum < 0)
     {
       drms_close_record(tempoutrec, DRMS_FREE_RECORD);
       return 1;
     }
   }
   else
   {
     fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
   }

  int itest;
// these must be present in the output dataseries and variable, not links or constants
/*
   char *outchecklist[] = {"T_REC", "QUALITY", "DATASIGN",
                           "CRPIX1", "CRVAL1", "CROTA1", "CDELT1",
                           "CRPIX2", "CRVAL2", "CROTA2", "CDELT2" };

   for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
   {
      DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
      if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1)
      {
         fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
         return 1;
      }
   }
*/
   drms_close_record(tempoutrec, DRMS_FREE_RECORD);


   DRMS_RecordSet_t *inRecSet = drms_open_records(drms_env, inRecQuery, &status);
   if (status != DRMS_SUCCESS || inRecSet == NULL)
   {
      fprintf(stderr, "ERROR: problem opening input recordset: status = %d\n", status);
      return 1;
   }

   nRecs = inRecSet->n;
   if (nRecs == 0)
   {
      printf("WARNING: input recordset contains no records\nmodule %s successful completion\n", cmdparams.argv[0]);
      return 1;
   }

   if (verbflag) 
      printf("input recordset opened, nRecs = %d\n",nRecs);


// go ahead and check for the presence of these input keywords as well
   char *inchecklist[] = {"T_REC", "QUALITY", "T_OBS", "CRLT_OBS", "CRLN_OBS", 
                          "CDELT1", "CDELT2"};

   DRMS_Record_t *tempinrec = inRecSet->records[0];
   for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
   {
      DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&tempinrec->keywords, inchecklist[itest]);
      if (!inkeytest)
      {
         fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
         return 1;
      }
   }


   int iRec;
   int error=0; // only set error before a break
   int nsuccess=0;
   for (iRec = 0; iRec < nRecs; iRec++)
   {
      if (verbflag > 1) 
      {
         wt1=getwalltime();
         getcputime(&ut1, &st1);
         gettimeofday(&tv1, NULL);
         ct1=dsecnd(); //((double)clock())/CLOCKS_PER_SEC;
         tt1=1000*((double)clock())/CLOCKS_PER_SEC;
         printf("processing record %d\n", iRec);
      }
      inrec = inRecSet->records[iRec];

      // create an output record
      outrec = drms_create_record(drms_env, 
                                  outSeries,
                                  lifetime, 
                                  &status);

      if (status != DRMS_SUCCESS || outrec==NULL)
      {
         fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status);
         error=2;
         break;
      }

      if (histlink)
         drms_setlink_static(outrec, histlinkname,  histrecnum);

      if (srclink)
         drms_setlink_static(outrec, srclinkname,  inrec->recnum);


      if (segoutflag)
         segout = drms_segment_lookup(outrec, segnameout);
      else
         segout = drms_segment_lookupnum(outrec, 0);

      // create an output array
      outarr = drms_array_create(usetype, 2, length, NULL, &status);

      if (!segout || !outarr || status != DRMS_SUCCESS)
      {
         fprintf(stderr, "ERROR: problem with output segment or array: iRec = %d, status = %d\n", iRec, status);
         if (outarr)
            drms_free_array(outarr);
         error=1;
         break;
      }

      trec = drms_getkey_time(inrec, "T_REC", &status);
      if (status != DRMS_SUCCESS)
      {
         fprintf(stderr,
        "SKIP: error getting prime keyword %s: iRec = %d, status = %d, recnum = %lld\n",
        "T_REC",
            iRec,
        status,
            inrec->recnum);
         drms_free_array(outarr);
         continue;
      }
      sprint_time(trecstr, trec, "TAI", 0); 

      quality = drms_getkey_int(inrec, "QUALITY", &status);
      PARAMETER_ERROR("QUALITY")
      // insert tests on quality here.
      // if we encounter a keyword error causing a continue within the loop, we should set a bit in quality for this.
      // in other words, CreateBlankRecord should always be preceded by setting quality, or could move this inside CreateBlankRecord.
      if (quality & QUAL_NODATA)
      {
         CreateBlankRecord(inrec, outrec, quality);
         fprintf(stderr,"SKIP: record rejected based on quality = %d: iRec = %d, T_REC = %s, recnum = %lld\n", quality, iRec, trecstr, inrec->recnum);
         drms_free_array(outarr);
         continue;
      }

      if (seginflag)
         segin = drms_segment_lookup(inrec, segnamein);
      else
         segin = drms_segment_lookupnum(inrec, 0);

      if (segin)
         inarr = drms_segment_read(segin, usetype, &status);

      if (!segin || !inarr || status != DRMS_SUCCESS)
      {
//         CreateBlankRecord(inrec, outrec, quality);
         fprintf(stderr, "ERROR: problem with input segment or array: iRec = %d, status = %d, T_REC = %s, recnum = %lld \n", iRec, status, trecstr, inrec->recnum);
         drms_free_array(outarr);
         if (inarr)
            drms_free_array(inarr);
         error=1;
         break;
      }

      if (maxmissvals > 0) 
      {
         int missvals = drms_getkey_int(inrec, "MISSVALS", &status);
         if (status == DRMS_ERROR_UNKNOWNKEYWORD)
         {
            fprintf(stderr, "ERROR: required keyword %s missing, iRec = %d, T_REC = %s, recnum = %lld\n", "MISSVALS", iRec, trecstr, inrec->recnum);
            drms_free_array(inarr);
            drms_free_array(outarr);
            free(orientation);
            error=1;
            break;
         }
         PARAMETER_ERROR("MISSVALS")
         if (missvals > maxmissvals) 
         {
            CreateBlankRecord(inrec, outrec, quality);
            fprintf(stderr, "SKIP: %d pixels MISSING, max allowed is %d: iRec = %d, status = %d, T_REC = %s, recnum = %lld\n", missvals, maxmissvals, iRec, status, trecstr, inrec->recnum);
            drms_free_array(inarr);
            drms_free_array(outarr);
            free(orientation);
            continue;
         }
      }

      tobs = drms_getkey_time(inrec, "T_OBS", &status);
      PARAMETER_ERROR("T_OBS")

      // MDI keyword was OBS_B0
      b0 = drms_getkey_double(inrec, "CRLT_OBS", &status);
      PARAMETER_ERROR("CRLT_OBS")

      // MDI keyword was OBS_L0
      obsl0 = drms_getkey_double(inrec, "CRLN_OBS", &status);
      PARAMETER_ERROR("CRLN_OBS")

      if (p0 == 999.0) 
      {
         // MDI keyword was SOLAR_P = -(SAT_ROT + INST_ROT)
/*
         satrot = drms_getkey_double(inrec, "SAT_ROT", &status);
         PARAMETER_ERROR("SAT_ROT")
         instrot = drms_getkey_double(inrec, "INST_ROT", &status);
         PARAMETER_ERROR("INST_ROT")
         p=-(satrot+instrot);
*/
         double crota = drms_getkey_double(inrec, "CROTA2", &status);
         PARAMETER_ERROR("CROTA2")
         p=-crota;
      } 
      else 
      {
         p = p0;
      }

      // fix for 1988 MWO
      // p = 180.0 - p ;

      p = psign * p ;

      // 991839442. corresponds to hour 73878 minute 57 second 22
      // or 73878.956 or day 3078.2898, roughly when B0 is 0
      // b0=b0 * 0.986207; One way of correcting. 
      // The following is pretty good 
      // b0=b0-0.1*sin((tobs-991839442.)/31557600.*2*PI);
      // b0 = b0 + ierr * sin((tobs - 991839442.) / 31557600. * 2 * PI);
      b0 = b0 + ierr * sin((tobs - trefb0) / 31557600. * 2 * PI);

      // p=p-0.2; 
      // p=p-0.2+0.1*cos((tobs-991839442.)/31557600.*2*PI); 
      // p = p + perr - ierr * cos((tobs - 991839442.) / 31557600. * 2 * PI);
      p = p + perr - ierr * cos((tobs - trefb0) / 31557600. * 2 * PI);

      if (paramsign != 0)
      {
         vsign = paramsign;
      }
      else
      {
         vsign = drms_getkey_double(inrec, "DATASIGN", &status);
         ParameterDef(status, "DATASIGN", 1.0, &vsign, iRec, 1);
      }

      xscale = drms_getkey_double(inrec, "CDELT1", &status);
      PARAMETER_ERROR("CDELT1")
      yscale = drms_getkey_double(inrec, "CDELT2", &status);
      PARAMETER_ERROR("CDELT2")

      if (xscale != yscale)
      {
        fprintf(stderr, "ERROR: %s != %s not supported, iRec = %d, T_REC = %s, recnum = %lld \n", "CDELT1", "CDELT2", iRec, trecstr, inrec->recnum);
        drms_free_array(inarr);
        error++;
        continue;
      }
      imagescale=xscale;

      // MDI keyword was OBS_DIST, in AU
      obsdist = drms_getkey_double(inrec, "DSUN_OBS", &status);
      obsdist /= AU;  // note that an incorrect value of 1.49597892e11 has sometimes been used to convert between OBS_DIST and DSUN_OBS when porting data from DSDS to DRMS 
      ParameterDef(status, "DSUN_OBS", 1.0, &obsdist, iRec, 1);
      S = rtrue / obsdist; // radians - approx. arcsin(rtrue/obsdist)


      rsun = drms_getkey_double(inrec, "R_SUN", &status);
      if (status != DRMS_SUCCESS)
      {
        rsun = drms_getkey_double(inrec, "RSUN_OBS", &status);
        if (status != DRMS_SUCCESS)
          rsun = ARCSECSINRAD * S / sqrt(1.0 - S * S) / imagescale;    
        else
          rsun /= imagescale;
      }
                   
      xpixels = inarr->axis[0];
      ypixels = inarr->axis[1];

      // MDI keyword was X0
      x0 = drms_getkey_double(inrec, "CRPIX1", &status);
      ParameterDef(status, "CRPIX1", xpixels / 2, &x0, iRec, 1);
      x0 -= 1.0;

      // MDI keyword was Y0
      y0 = drms_getkey_double(inrec, "CRPIX2", &status);
      ParameterDef(status, "CRPIX2", ypixels / 2, &y0, iRec, 1);
      y0 -= 1.0;

      if (status = imageinterp((float *)inarr->data, 
                              (float *)outarr->data,
                              xpixels, 
                              ypixels, 
                              x0, 
                              y0, 
                              p * RADSINDEG, 
                              rsun, 
                              rmax, 
                              NaN_beyond_rmax, 
                              interpolation, 
                              mapcols, 
                              maprows, 
                              vsign, 
                              &distP)) 
      {
         CreateBlankRecord(inrec, outrec, quality);
         fprintf(stderr, "failure in imageinterp: iRec = %d, status = %d, T_REC = %s, recnum = %lld\n", iRec, status, trecstr, inrec->recnum);
         drms_free_array(inarr);
         drms_free_array(outarr);
         free(orientation);
         continue; // go to next image 
      }

      drms_setkey_int(outrec, "TOTVALS", maprows*mapcols);
      set_statistics(segout, outarr, 1);

      // Write to the output record
      outarr->bzero=bias;
      outarr->bscale=scale;
      segout->axis[0]=outarr->axis[0]; // required for vardim output
      segout->axis[1]=outarr->axis[1];
      drms_segment_write(segout, outarr, 0);

//      drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
      drms_copykey(outrec, inrec, "T_REC");
      drms_setkey_int(outrec, "QUALITY", quality);

      double xratio = (double)mapcols/xpixels;
      double yratio = (double)maprows/ypixels;

      drms_setkey_float(outrec, "CRPIX1", (drms_getkey_float(inrec, "CRPIX1", &status)-0.5)*xratio+0.5);
      drms_setkey_float(outrec, "CRPIX2", (drms_getkey_float(inrec, "CRPIX2", &status)-0.5)*yratio+0.5);
      drms_setkey_float(outrec, "CDELT1", drms_getkey_float(inrec, "CDELT1", &status)/xratio);
      drms_setkey_float(outrec, "CDELT2", drms_getkey_float(inrec, "CDELT2", &status)/yratio);
      drms_setkey_float(outrec, "CROTA2", 0.0);
      drms_setkey_float(outrec, "X0", (x0 + 0.5)*xratio - 0.5);
      drms_setkey_float(outrec, "Y0", (y0 + 0.5)*yratio - 0.5);
      drms_setkey_float(outrec, "IM_SCALE",  imagescale/xratio);
      drms_setkey_float(outrec, "R_SUN",  rsun*xratio);
      drms_setkey_float(outrec, "RSUN_OBS", drms_getkey_float(inrec, "RSUN_OBS", &status)*xratio);

      drms_copykey(outrec, inrec, "CRVAL1");
      drms_copykey(outrec, inrec, "CRVAL2");
      drms_copykey(outrec, inrec, "T_OBS");
      drms_copykey(outrec, inrec, "CADENCE");
      drms_copykey(outrec, inrec, "EXPTIME");
      drms_copykey(outrec, inrec, "CRLT_OBS");
      drms_copykey(outrec, inrec, "CRLN_OBS");
      drms_copykey(outrec, inrec, "CAR_ROT");
      drms_copykey(outrec, inrec, "DSUN_OBS");
      drms_copykey(outrec, inrec, "R_SUN_REF");
      drms_copykey(outrec, inrec, "OBS_VR");
      drms_copykey(outrec, inrec, "OBS_VW");
      drms_copykey(outrec, inrec, "OBS_VN");

      dsignout=vsign*drms_getkey_double(inrec, "DATASIGN", &status);
      if (status != DRMS_SUCCESS) dsignout=vsign;
      dsignout/=fabs(dsignout);
      drms_setkey_int(outrec, "DATASIGN", (int)dsignout);

      tnow = (double)time(NULL);
      tnow += UNIX_epoch;
      drms_setkey_time(outrec, "DATE", tnow);

      drms_close_record(outrec, DRMS_INSERT_RECORD);

      if (verbflag > 1) 
      {
         wt=getwalltime();
         getcputime(&ut, &st);
         gettimeofday(&tv, NULL);
         ct=dsecnd(); //((double)clock())/CLOCKS_PER_SEC;
         tt=1000*((double)clock())/CLOCKS_PER_SEC;
         fprintf(stdout, 
                "record %d done, %f ms wall time, %f ms cpu time\n", 
                iRec, 
                (float)((tv.tv_sec * 1000000.0 + tv.tv_usec -
                        (tv1.tv_sec * 1000000.0 + tv1.tv_usec)) / 1000.0),
                (ct-ct1)*1000);
printf("test: %f ms wall time, %f ms cpu time\n", wt-wt1, (ut+st)-(ut1+st1));
printf("clock: %f ms\n",tt-tt1);
      }

      drms_free_array(inarr);
      drms_free_array(outarr);
      free(orientation);
      nsuccess++;

   } // end loop

   if (inRecSet)
      drms_close_records(inRecSet, DRMS_FREE_RECORD);

   if (error == 1)
      drms_close_record(outrec, DRMS_FREE_RECORD);

   wt=getwalltime();
   getcputime(&ut, &st);

   gettimeofday(&tv, NULL);
   ct=dsecnd(); //((double)clock())/CLOCKS_PER_SEC;
   tt=1000*((double)clock())/CLOCKS_PER_SEC;
   if (verbflag) 
   {
      printf("number of records processed = %d\n", nsuccess);
      fprintf(stdout, "total time spent: %f ms wall time, %f ms cpu time\n", 
         (float)((tv.tv_sec * 1000000.0 + tv.tv_usec - (tv0.tv_sec * 1000000.0 + tv0.tv_usec)) / 1000.0),
             (ct-ct0)*1000);
printf("test: %f ms wall time, %f ms cpu time\n", wt-wt0, (ut+st)-(ut0+st0));
printf("clock: %f ms\n",tt-tt0);
      if (!error)
         printf("module %s successful completion\n", cmdparams.argv[0]);
   }

   return 0;

} // end DoIt

---