(file) Return to sharp.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

Diff for /JSOC/proj/sharp/apps/sharp.c between version 1.11 and 1.31

version 1.11, 2013/01/24 00:05:58 version 1.31, 2014/06/16 04:45:18
Line 1 
Line 1 
 /* /*
  *  sharp.c  *  sharp.c
  *  *
  *      This module creates the pipeline space weather harps   *      This module creates the pipeline Space Weather Active Region Patches (SHARPs).
  *      It is a hard-coded strip-down version of bmap.c   *      It is a hard-coded strip-down version of bmap.c.
  *      It takes the Mharp and Bharp series and crete the following quantities   *      It takes the Mharp and Bharp series and create the following quantities:
    *
  *  Series 1: Sharp_CEA  *  Series 1: Sharp_CEA
  *        CEA remapped magnetogram, bitmap, continuum, doppler (same size in map coordinate, need manual spec?)  *        CEA remapped magnetogram, bitmap, continuum, doppler (same size in map coordinate, need manual spec?)
  *        CEA remapped vector field (Br, Bt, Bp) (same as above)  *        CEA remapped vector field (Br, Bt, Bp) (same as above)
  *    Space weather indices based on vector cutouts (step 2)  *    Space weather indices based on vector cutouts (step 2)
    *
  *  Series 2: Sharp_cutout:  *  Series 2: Sharp_cutout:
  *        cutouts of magnetogram, bitmap, continuum, doppler (HARP defined, various sizes in CCD pixels)  *        cutouts of magnetogram, bitmap, continuum, doppler (HARP defined, various sizes in CCD pixels)
  *        cutouts of all vector data segments (same as above)  *        cutouts of all vector data segments (same as above)
  *      Series 3: Other remaps  
  *  *
  *      Author:  *      Author:
  *              Xudong Sun; Monica Bobra  *              Xudong Sun; Monica Bobra
Line 22 
Line 23 
  *              v0.2    Sep 04 2012  *              v0.2    Sep 04 2012
  *              v0.3    Dec 18 2012  *              v0.3    Dec 18 2012
  *              v0.4    Jan 02 2013  *              v0.4    Jan 02 2013
  *              v0.5    Jam 23 2012   *              v0.5 Jan 23 2013
    *              v0.6 Aug 12 2013
    *              v0.7 Jan 02 2014
    *              v0.8 Feb 12 2014
    *                              v0.9 Mar 04 2014
  *  *
  *      Notes:  *      Notes:
  *              v0.0  *              v0.0
Line 42 
Line 47 
  *              Took out convert_inplace(). Was causing all the images to be int  *              Took out convert_inplace(). Was causing all the images to be int
  *              v0.5  *              v0.5
  *              Corrected ephemeris keywords, added argument mInfo for setKeys()  *              Corrected ephemeris keywords, added argument mInfo for setKeys()
  *   *              v0.6
  *      Example:   *              Changes in remapping of bitmap and conf_disambig, now near neighbor without anti-aliasing
  *      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \   *              v0.7
  "bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \   *              Added full disk as "b"
  "dop=hmi.V_720s[2012.02.20_10:00]" \   *              Global flag fullDisk is set if "b" is set
  "cont=hmi.Ic_720s[2012.02.20_10:00]" \   *              Utilize BharpRS and BharpRec all around
  "sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"   *              Pass mharpRec to set_keys() too in case of full disk
  *      For comparison:   *              Fixed Bunit (removed from copy_me_keys(), added loops for Bunits in set_keys() here)
  *      bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \   *              Error for CEA still does account for disambiguation yet
  "out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"   *              v0.8
    *              Added disambig to azimuth during error propagation
    *              Changed usage for disambig: bit 2 (radial acute) for full disk, bit 0 for patch
    *              Fixed disambig cutout for patch: 0 for even, 7 for odd
    *              v0.9
    *              Fixed unit
    *              Check whether in PATCH of FD mode, so the error propagation uses disambiguated azimuth or not
    *
    *
    *      Example Calls:
    *      [I]   B (full disk disambiguation)
    *      > sharp "mharp=hmi.Mharp_720s[1832][2012.07.12_15:24]" "b=hmi_test.B_720s[2012.07.12_15:24]" "dop=hmi.V_720s[2012.07.12_15:24]" "cont=hmi.Ic_720s[2012.07.12_15:24]" "sharp_cea=su_xudong.sharp_cea_720s" "sharp_cut=su_xudong.sharp_720s"
    *      [II]  BHARP (patch disambiguation)
    *      > sharp "mharp=hmi.Mharp_720s[1832][2012.07.12_15:24]" "bharp=hmi.Bharp_720s[1832][2012.07.12_15:24]" "dop=hmi.V_720s[2012.07.12_15:24]" "cont=hmi.Ic_720s[2012.07.12_15:24]" "sharp_cea=su_xudong.sharp_cea_720s" "sharp_cut=su_xudong.sharp_720s"
  *  *
  *  *
  */  */
Line 73 
Line 91 
 #include "errorprop.c" #include "errorprop.c"
 #include "sw_functions.c" #include "sw_functions.c"
  
   //#include <mkl.h> // Comment out mkl.h, which can only run on solar3
 #include <mkl_blas.h> #include <mkl_blas.h>
 #include <mkl_service.h> #include <mkl_service.h>
 #include <mkl_lapack.h> #include <mkl_lapack.h>
Line 91 
Line 110 
 // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
 #define NYQVIST         (0.015) #define NYQVIST         (0.015)
  
   // Some other things
 #ifndef MIN #ifndef MIN
 #define MIN(a,b) (((a)<(b)) ? (a) : (b)) #define MIN(a,b) (((a)<(b)) ? (a) : (b))
 #endif #endif
Line 118 
Line 138 
 #define INTERP                  0 #define INTERP                  0
 #define dpath    "/home/jsoc/cvs/Development/JSOC" #define dpath    "/home/jsoc/cvs/Development/JSOC"
  
   
 /* ========================================================================================================== */ /* ========================================================================================================== */
  
 // Space weather keywords // Space weather keywords
 struct swIndex { struct swIndex {
         float mean_vf;         float mean_vf;
       float count_mask;
         float absFlux;         float absFlux;
         float mean_hf;         float mean_hf;
         float mean_gamma;         float mean_gamma;
Line 142  struct swIndex {
Line 164  struct swIndex {
         float meanshear_angle;         float meanshear_angle;
         float area_w_shear_gt_45h;         float area_w_shear_gt_45h;
         float meanshear_angleh;         float meanshear_angleh;
       float mean_derivative_btotal_err;
       float mean_vf_err;
       float mean_gamma_err;
       float mean_derivative_bh_err;
       float mean_derivative_bz_err;
       float mean_jz_err;
       float us_i_err;
       float mean_alpha_err;
       float mean_ih_err;
       float total_us_ih_err;
       float total_abs_ih_err;
       float totaljz_err;
       float meanpot_err;
       float totpot_err;
       float meanshear_angle_err;
       float Rparam;
       float totfx;
       float totfy;
       float totfz;
       float totbsq;
       float epsx;
       float epsy;
       float epsz;
 }; };
  
 // Mapping method // Mapping method
Line 162  enum projection {
Line 207  enum projection {
 char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG", char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
         "SIN", "ZEA"};         "SIN", "ZEA"};
  
 // Ephemeris  // Ephemeris information
 struct ephemeris { struct ephemeris {
         double disk_lonc, disk_latc;         double disk_lonc, disk_latc;
         double disk_xc, disk_yc;         double disk_xc, disk_yc;
Line 180  struct mapInfo {
Line 225  struct mapInfo {
         float *xi_out, *zeta_out;       // coordinate on full disk image to sample at         float *xi_out, *zeta_out;       // coordinate on full disk image to sample at
 }; };
  
   
 /* ========================================================================================================== */ /* ========================================================================================================== */
  
 /* Get all input data series */ /* Get all input data series */
Line 196  int getInputRS_aux(DRMS_RecordSet_t **in
Line 240  int getInputRS_aux(DRMS_RecordSet_t **in
 /* Find record from record set with given T_rec */ /* Find record from record set with given T_rec */
 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec); int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
  
 // ===================  
   
 /* Create CEA record */ /* Create CEA record */
 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
                                         DRMS_Record_t *dopRec, DRMS_Record_t *contRec,                                         DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
Line 223  int getEphemeris(DRMS_Record_t *inRec, s
Line 265  int getEphemeris(DRMS_Record_t *inRec, s
 void findCoord(struct mapInfo *mInfo); void findCoord(struct mapInfo *mInfo);
  
 /* Mapping function */ /* Mapping function */
 int performSampling(float *outData, float *inData, struct mapInfo *mInfo);  int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
  
 /* Performing local vector transformation */ /* Performing local vector transformation */
 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo); void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
Line 260  void computeSWIndex(struct swIndex *swKe
Line 302  void computeSWIndex(struct swIndex *swKe
 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr); void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
  
 /* Set all keywords, no error checking for now */ /* Set all keywords, no error checking for now */
 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo);  // Changed Dec 30 XS
   void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
  
 // =================== // ===================
  
Line 300  char *CutSegs[] = {"magnetogram", "bitma
Line 343  char *CutSegs[] = {"magnetogram", "bitma
         "field_inclination_err", "field_az_err", "inclin_azimuth_err",         "field_inclination_err", "field_az_err", "inclin_azimuth_err",
         "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",         "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
         "disambig", "conf_disambig"};         "disambig", "conf_disambig"};
 char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",  char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
         BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};          BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA, "conf_disambig"};
 /* ========================================================================================================== */  // For Bunits, added Dec 30 XS
   char *CutBunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
       "degree", "degree", "Mx/cm^2", "cm/s", "mA", " ",
       "length units", "DN/s", "DN/s", " ", " ",
       " ",
       " ", " ",
       "degree", "degree", "Mx/cm^2", "cm/s", " ",
       " ", " ", " ",
       " ", " ", " ",
       " ", " "};
   char *CEABunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
       "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", " "};      // Mar 4 2014 XS
  
   /* ========================================================================================================== */
  
 char *module_name = "sharp"; char *module_name = "sharp";
 char *version_id = "2012 Dec 18";  /* Version number */  int seed;
   
   int fullDisk;       // full disk mode
   int amb4err;      // Use azimuth disambiguation for error propagation, default is 0 for patch and 1 for FD
  
 ModuleArgs_t module_args[] = ModuleArgs_t module_args[] =
 { {
         {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},         {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
         {ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},         {ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},
       {ARG_STRING, "b", kNotSpecified, "Input B series, if set, overrides bharp."},
         {ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},         {ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},
         {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},         {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
         {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},         {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
         {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},         {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
       {ARG_INT,    "seed", "987654", "Seed for the random number generator."},
       {ARG_INT,   "f_amb4err", "0", "Force using disambiguation in error propagation"},     // Mar 4 2014 XS
         {ARG_END}         {ARG_END}
 }; };
  
 int DoIt(void) int DoIt(void)
 { {
       int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
       int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
  
         int status = DRMS_SUCCESS;         int status = DRMS_SUCCESS;
         int nrecs, irec;         int nrecs, irec;
  
         char *mharpQuery, *bharpQuery;          char *mharpQuery, *bharpQuery, *bQuery;
         char *dopQuery, *contQuery;         char *dopQuery, *contQuery;
         char *sharpCeaQuery, *sharpCutQuery;         char *sharpCeaQuery, *sharpCutQuery;
  
Line 337  int DoIt(void)
Line 399  int DoIt(void)
  
         mharpQuery = (char *) params_get_str(&cmdparams, "mharp");         mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
         bharpQuery = (char *) params_get_str(&cmdparams, "bharp");         bharpQuery = (char *) params_get_str(&cmdparams, "bharp");
       bQuery = (char *) params_get_str(&cmdparams, "b");
         dopQuery = (char *) params_get_str(&cmdparams, "dop");         dopQuery = (char *) params_get_str(&cmdparams, "dop");
         contQuery = (char *) params_get_str(&cmdparams, "cont");         contQuery = (char *) params_get_str(&cmdparams, "cont");
         sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");         sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
         sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");         sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
  
       seed = params_get_int(&cmdparams, "seed");
       int f_amb4err = params_get_int(&cmdparams, "f_amb4err");
   
         /* Get input data, check everything */         /* Get input data, check everything */
  
       // Full disk mode if "b" is set
       if (strcmp(bQuery, kNotSpecified)) {
           fullDisk = 1;
           bharpQuery = bQuery;
           //        SHOW(bharpQuery); SHOW("\n");
           SHOW("Full disk mode\n");
       } else {
           fullDisk = 0;
           SHOW("Harp mode\n");
       }
   
       // Mar 4 2014
       if (f_amb4err == 0) {         // no forcing, 0 for patch and 1 for FD
           amb4err = fullDisk ? 1 : 0;
       } else {
           amb4err = 1;
       }
       printf("amb4err=%d\n", amb4err);
   
       // Bharp point to B if full disk
         if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))         if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
                 DIE("Input harp data error.");                 DIE("Input harp data error.");
         nrecs = mharpRS->n;         nrecs = mharpRS->n;
Line 363  int DoIt(void)
Line 449  int DoIt(void)
                 /* Records in work */                 /* Records in work */
  
                 DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;                 DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;
   
                 mharpRec = mharpRS->records[irec];                 mharpRec = mharpRS->records[irec];
                 bharpRec = bharpRS->records[irec];  
  
                 TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);                 TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
  
                   if (!fullDisk) {
               bharpRec = bharpRS->records[irec];
           } else {
               if (getInputRec_aux(&bharpRec, bharpRS, trec)) {     // Bharp point to full disk B
                   printf("Fetching B failed, image #%d skipped.\n", irec);
                   continue;
               }
           }
   
                 struct swIndex swKeys;                 struct swIndex swKeys;
  
                 DRMS_Record_t *dopRec = NULL, *contRec = NULL;                 DRMS_Record_t *dopRec = NULL, *contRec = NULL;
   
                 if (getInputRec_aux(&dopRec, dopRS, trec)) {                 if (getInputRec_aux(&dopRec, dopRS, trec)) {
                         printf("Fetching Doppler failed, image #%d skipped.\n", irec);                         printf("Fetching Doppler failed, image #%d skipped.\n", irec);
                         continue;                         continue;
Line 449  int getInputRS(DRMS_RecordSet_t **mharpR
Line 545  int getInputRS(DRMS_RecordSet_t **mharpR
         *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);         *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
     if (status || (*mharpRS_ptr)->n == 0) return 1;     if (status || (*mharpRS_ptr)->n == 0) return 1;
  
           if (fullDisk) {
           if (getInputRS_aux(bharpRS_ptr, bharpQuery, *mharpRS_ptr)) return 1;
       } else {
         *bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);         *bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
     if (status || (*bharpRS_ptr)->n == 0) return 1;     if (status || (*bharpRS_ptr)->n == 0) return 1;
   
         if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;         if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
       }
  
         return 0;         return 0;
  
Line 578  int createCeaRecord(DRMS_Record_t *mharp
Line 677  int createCeaRecord(DRMS_Record_t *mharp
         mInfo.proj = (enum projection) cyleqa;          // projection method         mInfo.proj = (enum projection) cyleqa;          // projection method
         mInfo.xscale = XSCALE;         mInfo.xscale = XSCALE;
         mInfo.yscale = YSCALE;         mInfo.yscale = YSCALE;
         mInfo.nbin = NBIN;  
       int ncol0, nrow0;           // oversampled map size
  
         // Get ephemeris         // Get ephemeris
  
Line 594  int createCeaRecord(DRMS_Record_t *mharp
Line 694  int createCeaRecord(DRMS_Record_t *mharp
                 return 1;                 return 1;
         }         }
  
           // ========================================
           // Do this for all bitmaps, Aug 12 2013 XS
           // ========================================
   
       mInfo.nbin = 1;                     // for bitmaps. suppress anti-aliasing
           ncol0 = mInfo.ncol;
           nrow0 = mInfo.nrow;
   
           mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
           mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
   
           findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
   
           if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
                   SHOW("CEA: mapping bitmap error\n");
                   return 1;
           }
           printf("Bitmap mapping done.\n");
   
       if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
                   SHOW("CEA: mapping conf_disambig error\n");
                   return 1;
           }
           printf("Conf disambig mapping done.\n");
   
       free(mInfo.xi_out);
           free(mInfo.zeta_out);
   
           // ========================================
           // Do this again for floats, Aug 12 2013 XS
           // ========================================
         // Create xi_out, zeta_out array in mInfo:         // Create xi_out, zeta_out array in mInfo:
         // Coordinates to sample in original full disk image         // Coordinates to sample in original full disk image
  
         int ncol0, nrow0;               // oversampled map size          mInfo.nbin = NBIN;
         ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN         ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
         nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;         nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
  
Line 612  int createCeaRecord(DRMS_Record_t *mharp
Line 743  int createCeaRecord(DRMS_Record_t *mharp
                 SHOW("CEA: mapping magnetogram error\n");                 SHOW("CEA: mapping magnetogram error\n");
                 return 1;                 return 1;
         }         }
         if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {  
                 SHOW("CEA: mapping bitmap error\n");  
                 return 1;  
         }  
         printf("Magnetogram mapping done.\n");         printf("Magnetogram mapping done.\n");
  
         if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {         if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
Line 630  int createCeaRecord(DRMS_Record_t *mharp
Line 757  int createCeaRecord(DRMS_Record_t *mharp
         }         }
         printf("Intensitygram mapping done.\n");         printf("Intensitygram mapping done.\n");
  
         if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {  
                 SHOW("CEA: mapping conf_disambig error\n");  
                 return 1;  
         }  
         printf("Conf disambig mapping done.\n");  
   
         // Mapping vector B         // Mapping vector B
  
         if (mapVectorB(sharpRec, bharpRec, &mInfo)) {         if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
Line 657  int createCeaRecord(DRMS_Record_t *mharp
Line 778  int createCeaRecord(DRMS_Record_t *mharp
         drms_copykey(sharpRec, mharpRec, "T_REC");         drms_copykey(sharpRec, mharpRec, "T_REC");
         drms_copykey(sharpRec, mharpRec, "HARPNUM");         drms_copykey(sharpRec, mharpRec, "HARPNUM");
  
         DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");      if (fullDisk) {
         if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);          DRMS_Link_t *bLink = hcon_lookup_lower(&sharpRec->links, "B");
           if (bLink) drms_link_set("B", sharpRec, bharpRec);
       } else {
         DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");         DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
         if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);         if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
       }
           DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
           if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
  
         setKeys(sharpRec, bharpRec, &mInfo);            // Set all other keywords      setKeys(sharpRec, mharpRec, bharpRec, &mInfo);            // Set all other keywords
         drms_copykey(sharpRec, mharpRec, "QUALITY");            // copied from los records         drms_copykey(sharpRec, mharpRec, "QUALITY");            // copied from los records
  
         // Space weather         // Space weather
Line 704  int mapScaler(DRMS_Record_t *sharpRec, D
Line 830  int mapScaler(DRMS_Record_t *sharpRec, D
         int status = 0;         int status = 0;
         int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;         int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
         int dims[2] = {nx, ny};         int dims[2] = {nx, ny};
           int interpOpt = INTERP;         // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
  
         // Input full disk array         // Input full disk array
  
Line 715  int mapScaler(DRMS_Record_t *sharpRec, D
Line 842  int mapScaler(DRMS_Record_t *sharpRec, D
         inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);         inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
         if (!inArray) return 1;         if (!inArray) return 1;
  
       if (!strcmp(segName, "conf_disambig") || !strcmp(segName, "bitmap")) {
           // Moved out so it works for FD conf_disambig as well
           // Jan 2 2014 XS
           interpOpt = 3;          // Aug 12 XS, near neighbor
       }
   
         float *inData;         float *inData;
         int xsz = inArray->axis[0], ysz = inArray->axis[1];         int xsz = inArray->axis[0], ysz = inArray->axis[1];
         if ((xsz != FOURK) || (ysz != FOURK)) {         // for bitmap, make tmp full disk         if ((xsz != FOURK) || (ysz != FOURK)) {         // for bitmap, make tmp full disk
Line 737  int mapScaler(DRMS_Record_t *sharpRec, D
Line 870  int mapScaler(DRMS_Record_t *sharpRec, D
         // Mapping         // Mapping
  
         float *map = (float *) (malloc(nxny * sizeof(float)));         float *map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(map, inData, mInfo))          if (performSampling(map, inData, mInfo, interpOpt))             // Add interpOpt for different types, Aug 12 XS
         {if (inArray) drms_free_array(inArray); free(map); return 1;}         {if (inArray) drms_free_array(inArray); free(map); return 1;}
  
         // Write out         // Write out
Line 800  int mapVectorB(DRMS_Record_t *sharpRec,
Line 933  int mapVectorB(DRMS_Record_t *sharpRec,
         float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation         float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation
  
         bx_map = (float *) (malloc(nxny * sizeof(float)));         bx_map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(bx_map, bx_img, mInfo))          if (performSampling(bx_map, bx_img, mInfo, INTERP))
         {free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}         {free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
  
         by_map = (float *) (malloc(nxny * sizeof(float)));         by_map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(by_map, by_img, mInfo))          if (performSampling(by_map, by_img, mInfo, INTERP))
         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
  
         bz_map = (float *) (malloc(nxny * sizeof(float)));         bz_map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(bz_map, bz_img, mInfo))          if (performSampling(bz_map, bz_img, mInfo, INTERP))
         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
  
         free(bx_img); free(by_img); free(bz_img);         free(bx_img); free(by_img); free(bz_img);
Line 916  int findPosition(DRMS_Record_t *inRec, s
Line 1049  int findPosition(DRMS_Record_t *inRec, s
         float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);         float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
  
         /* Center coord */         /* Center coord */
       // Changed into double Jun 16 2014 XS
  
         float minlon = drms_getkey_float(inRec, "LONDTMIN", &status); if (status) return 1;             // Stonyhurst lon          double minlon = drms_getkey_double(inRec, "LONDTMIN", &status); if (status) return 1;           // Stonyhurst lon
         float maxlon = drms_getkey_float(inRec, "LONDTMAX", &status); if (status) return 1;          double maxlon = drms_getkey_double(inRec, "LONDTMAX", &status); if (status) return 1;
         float minlat = drms_getkey_float(inRec, "LATDTMIN", &status); if (status) return 1;          double minlat = drms_getkey_double(inRec, "LATDTMIN", &status); if (status) return 1;
         float maxlat = drms_getkey_float(inRec, "LATDTMAX", &status); if (status) return 1;          double maxlat = drms_getkey_double(inRec, "LATDTMAX", &status); if (status) return 1;
  
         // A bug fixer for HARP (per M. Turmon)         // A bug fixer for HARP (per M. Turmon)
         // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong         // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
Line 928  int findPosition(DRMS_Record_t *inRec, s
Line 1062  int findPosition(DRMS_Record_t *inRec, s
         // We compute minlon & minlat then by         // We compute minlon & minlat then by
         // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT         // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
  
         float psize = drms_getkey_float(inRec, "SIZE", &status);      //  float psize = drms_getkey_float(inRec, "SIZE", &status);
         if (psize != psize) {      //  if (psize != psize) {
                 TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1;  
       if (minlon != minlon || maxlon != maxlon) {         // check lons instead of SIZE
                   TIME t0 = drms_getkey_time(inRec, "T_FRST1", &status); if (status) return 1;                    // changed from T_FRST to T_FRST1, T_FRST may not exist
                 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;                 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
                 char firstRecQuery[100], t0_str[100];                 char firstRecQuery[100], t0_str[100];
                 sprint_time(t0_str, t0, "TAI", 0);                 sprint_time(t0_str, t0, "TAI", 0);
Line 949  int findPosition(DRMS_Record_t *inRec, s
Line 1085  int findPosition(DRMS_Record_t *inRec, s
         mInfo->yc = (maxlat + minlat) / 2.;         mInfo->yc = (maxlat + minlat) / 2.;
  
         /* Size */         /* Size */
       // Rounded to 1.d3 precision first. Jun 16 2014 XS
  
         mInfo->ncol = round((maxlon - minlon) / mInfo->xscale);          mInfo->ncol = round(round((maxlon - minlon) * 1.e3) / 1.e3 / mInfo->xscale);
         mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);          mInfo->nrow = round(round((maxlat - minlat) * 1.e3) / 1.e3 / mInfo->yscale);
  
         return 0;         return 0;
  
Line 982  int getEphemeris(DRMS_Record_t *inRec, s
Line 1119  int getEphemeris(DRMS_Record_t *inRec, s
         float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);         float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
         float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);         float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
         float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy         float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
         printf("cdelt=%f\n",cdelt);  
         ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0         ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0
         ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;         ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
  
Line 1086  void findCoord(struct mapInfo *mInfo)
Line 1222  void findCoord(struct mapInfo *mInfo)
  *  *
  */  */
  
 int performSampling(float *outData, float *inData, struct mapInfo *mInfo)  int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
 { {
  
         int status = 0;         int status = 0;
           int ind_map;
  
         int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN         int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN
         int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;         int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
  
         float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));          // Changed Aug 12 2013, XS, for bitmaps
           float *outData0;
           if (interpOpt == 3 && mInfo->nbin == 1) {
           outData0 = outData;
           } else {
           outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
           }
  
         float *xi_out = mInfo->xi_out;         float *xi_out = mInfo->xi_out;
         float *zeta_out = mInfo->zeta_out;         float *zeta_out = mInfo->zeta_out;
Line 1102  int performSampling(float *outData, floa
Line 1245  int performSampling(float *outData, floa
         // Interpolation         // Interpolation
  
         struct fint_struct pars;         struct fint_struct pars;
         int interpOpt = INTERP;         // Use Wiener by default, 6 order, 1 constraint          // Aug 12 2013, passed in as argument now
  
         switch (interpOpt) {         switch (interpOpt) {
                 case 0:                 // Wiener, 6 order, 1 constraint                 case 0:                 // Wiener, 6 order, 1 constraint
Line 1114  int performSampling(float *outData, floa
Line 1257  int performSampling(float *outData, floa
                 case 2:                 // Bilinear                 case 2:                 // Bilinear
                         init_finterpolate_linear(&pars, 1.);                         init_finterpolate_linear(&pars, 1.);
                         break;                         break;
                   case 3:                 // Near neighbor
               break;
                 default:                 default:
                         return 1;                         return 1;
         }         }
  
           printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
           if (interpOpt == 3) {                   // Aug 6 2013, Xudong
                   for (int row0 = 0; row0 < nrow0; row0++) {
               for (int col0 = 0; col0 < ncol0; col0++) {
                   ind_map = row0 * ncol0 + col0;
                   outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
               }
           }
           } else {
         finterpolate(&pars, inData, xi_out, zeta_out, outData0,         finterpolate(&pars, inData, xi_out, zeta_out, outData0,
                                  FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);                                  FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
           }
  
         // Rebinning, smoothing         // Rebinning, smoothing
  
           if (interpOpt == 3 && mInfo->nbin == 1) {
           return 0;
           } else {
         frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian         frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
         free(outData0);         // Dec 18 2012         free(outData0);         // Dec 18 2012
           }
  
         //         //
  
Line 1282  int getBErr(float *bx_err, float *by_err
Line 1441  int getBErr(float *bx_err, float *by_err
                         if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,                         if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
                                                         disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {                                                         disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
                                 bx_err[ind_map] = DRMS_MISSING_FLOAT;                                 bx_err[ind_map] = DRMS_MISSING_FLOAT;
                                 bx_err[ind_map] = DRMS_MISSING_FLOAT;                                  by_err[ind_map] = DRMS_MISSING_FLOAT;
                                 bx_err[ind_map] = DRMS_MISSING_FLOAT;                                  bz_err[ind_map] = DRMS_MISSING_FLOAT;       // Mar 7
                                 continue;                                 continue;
                         }                         }
  
Line 1363  int readVectorB(DRMS_Record_t *inRec, fl
Line 1522  int readVectorB(DRMS_Record_t *inRec, fl
         int llx, lly;           // lower-left corner         int llx, lly;           // lower-left corner
         int bmx, bmy;           // bitmap size         int bmx, bmy;           // bitmap size
  
       if (fullDisk) {
           llx = lly = 0;
           bmx = bmy = FOURK;
       } else {
         llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;         llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
         lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;         lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
   
         bmx = inArray_ambig->axis[0];         bmx = inArray_ambig->axis[0];
         bmy = inArray_ambig->axis[1];         bmy = inArray_ambig->axis[1];
       }
  
         int kx, ky, kOff;         int kx, ky, kOff;
         int ix = 0, jy = 0, yOff = 0, iData = 0;         int ix = 0, jy = 0, yOff = 0, iData = 0;
         int xDim = FOURK, yDim = FOURK;         int xDim = FOURK, yDim = FOURK;
           int amb = 0;
  
         for (jy = 0; jy < yDim; jy++)         for (jy = 0; jy < yDim; jy++)
         {         {
Line 1394  int readVectorB(DRMS_Record_t *inRec, fl
Line 1558  int readVectorB(DRMS_Record_t *inRec, fl
                                 continue;                                 continue;
                         } else {                         } else {
                                 kOff = ky * bmx + kx;                                 kOff = ky * bmx + kx;
                                 if (ambig[kOff] % 2) {          // 180                  //                              if (ambig[kOff] % 2) {          // 180
                                   // Feb 12 2014, use bit #2 for full disk, lowest bit for patch
                                   if (fullDisk) { amb = (ambig[kOff] / 4) % 2; } else { amb = ambig[kOff] % 2; }
                                   if (amb) {                              // Feb 12 2014, use bit #2
                                         bx_img[iData] *= -1.; by_img[iData] *= -1.;                                         bx_img[iData] *= -1.; by_img[iData] *= -1.;
                                 }                                 }
                         }                         }
Line 1434  int readVectorBErr(DRMS_Record_t *inRec,
Line 1601  int readVectorBErr(DRMS_Record_t *inRec,
         DRMS_Array_t *inArrays[9];         DRMS_Array_t *inArrays[9];
  
         // Read full disk images         // Read full disk images
       // Do we need disambig? Dec 30 XS
  
         for (int iSeg = 0; iSeg < 9; iSeg++) {         for (int iSeg = 0; iSeg < 9; iSeg++) {
  
Line 1447  int readVectorBErr(DRMS_Record_t *inRec,
Line 1615  int readVectorBErr(DRMS_Record_t *inRec,
         float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];         float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];
         float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];         float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];
  
           // Add disambig, Feb 12 2014
   
           DRMS_Segment_t *inSeg;
       DRMS_Array_t *inArray_ambig;
   
       if (amb4err) {              // Mar 4 2014
   
           inSeg = drms_segment_lookup(inRec, "disambig");
           inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
           if (status) return 1;
           char *ambig = (char *)inArray_ambig->data;
   
           int llx, lly;           // lower-left corner
           int bmx, bmy;           // bitmap size
   
           if (fullDisk) {
               llx = lly = 0;
               bmx = bmy = FOURK;
           } else {
               llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
               lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
               bmx = inArray_ambig->axis[0];
               bmy = inArray_ambig->axis[1];
           }
   
           int idx, idx_a;
           int amb;
   
           for (int j = 0; j < bmy; j++) {
               for (int i = 0; i < bmx; i++) {
                   idx_a = j * bmx + i;
                   idx = (j + lly) * FOURK + (i + llx);
                   // Feb 12 2014, use bit #2 for full disk, lowest bit for patch
                   if (fullDisk) { amb = (ambig[idx_a] / 4) % 2; } else { amb = ambig[idx_a] % 2; }
                   if (amb) { bA0[idx] += 180.; }
               }
           }
   
       }
   
         // Convert errors to variances, correlation coefficients to covariances         // Convert errors to variances, correlation coefficients to covariances
  
         for (int i = 0; i < FOURK2; i++) {         for (int i = 0; i < FOURK2; i++) {
Line 1455  int readVectorBErr(DRMS_Record_t *inRec,
Line 1663  int readVectorBErr(DRMS_Record_t *inRec,
                 if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;                 if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;
  
                 bT[i] = bT0[i];                 bT[i] = bT0[i];
                 bI[i] = bI0[i];                  bI[i] = bI0[i];         // in deg, coverted in errorprop
                 bA[i] = bA0[i];                 bA[i] = bA0[i];
  
                 errbT[i] = errbT0[i] * errbT0[i];                 errbT[i] = errbT0[i] * errbT0[i];
Line 1471  int readVectorBErr(DRMS_Record_t *inRec,
Line 1679  int readVectorBErr(DRMS_Record_t *inRec,
         //         //
  
         for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);         for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
           if (amb4err) drms_free_array(inArray_ambig);            // Feb 12; Mar 04 2014
  
         return 0;         return 0;
  
Line 1546  int createCutRecord(DRMS_Record_t *mharp
Line 1755  int createCutRecord(DRMS_Record_t *mharp
         if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);         if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
  
         setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices         setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
         setKeys(sharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout          setKeys(sharpRec, mharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout
  
         // Stats         // Stats
  
Line 1603  int writeCutout(DRMS_Record_t *outRec, D
Line 1812  int writeCutout(DRMS_Record_t *outRec, D
                 return 1;                 return 1;
         }         }
  
           // Feb 12 2014, fool-proof, for patch, change everything to 0 or 7!!!
           // This is a fix for disambiguation before Aug 2013
   
           if (!strcmp(SegName, "disambig") && !fullDisk) {
                   double *disamb = (double *) (cutoutArray->data);
                   for (int i = 0; i < nxny; i++) {
                           if (((int)disamb[i]) % 2) { disamb[i] = 7; } else { disamb[i] = 0; }
                   }
           }
   
         /* Adding disambiguation resolution to cutout azimuth? */         /* Adding disambiguation resolution to cutout azimuth? */
  
 #if DISAMB_AZI #if DISAMB_AZI
           int amb;
         if (!strcmp(SegName, "azimuth")) {         if (!strcmp(SegName, "azimuth")) {
                 DRMS_Segment_t *disambSeg = NULL;                 DRMS_Segment_t *disambSeg = NULL;
                 disambSeg = drms_segment_lookup(inRec, "disambig");                 disambSeg = drms_segment_lookup(inRec, "disambig");
                 if (!disambSeg) {drms_free_array(cutoutArray); return 1;}                 if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
                 DRMS_Array_t *disambArray;                 DRMS_Array_t *disambArray;
           if (fullDisk) { // Jan 2 2014 XS
               disambArray = drms_segment_readslice(disambSeg, DRMS_TYPE_CHAR, ll, ur, &status);
               if (status) return 1;
           } else {
                 if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {                 if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
                         disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);                         disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
                         if (status) {drms_free_array(cutoutArray); return 1;}                         if (status) {drms_free_array(cutoutArray); return 1;}
Line 1618  int writeCutout(DRMS_Record_t *outRec, D
Line 1842  int writeCutout(DRMS_Record_t *outRec, D
                   drms_free_array(cutoutArray);                   drms_free_array(cutoutArray);
                         return 1;                         return 1;
                 }                 }
           }
                 double *azimuth = (double *) cutoutArray->data;                 double *azimuth = (double *) cutoutArray->data;
                 char *disamb = (char *) disambArray->data;                 char *disamb = (char *) disambArray->data;
                 for (int n = 0; n < nxny; n++) {                 for (int n = 0; n < nxny; n++) {
                         if (disamb[n]) azimuth[n] += 180.;              //                  if (disamb[n] % 2) azimuth[n] += 180.;      // Nov 12 2013 Fixed!!!
                           // Feb 12 2014, use bit #2 for full disk, lowest bit for patch
                           if (fullDisk) { amb = (disamb[n] / 4) % 2; } else { amb = disamb[n] % 2; }
                           if (amb) azimuth[n] += 180.;
                 }                 }
                 drms_free_array(disambArray);                 drms_free_array(disambArray);
         }         }
Line 1661  void computeSWIndex(struct swIndex *swKe
Line 1889  void computeSWIndex(struct swIndex *swKe
         int nx = mInfo->ncol, ny = mInfo->nrow;         int nx = mInfo->ncol, ny = mInfo->nrow;
         int nxny = nx * ny;         int nxny = nx * ny;
         int dims[2] = {nx, ny};         int dims[2] = {nx, ny};
   
         // Get bx, by, bz, mask         // Get bx, by, bz, mask
  
         // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities         // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
Line 1686  void computeSWIndex(struct swIndex *swKe
Line 1915  void computeSWIndex(struct swIndex *swKe
         DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);         DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
         float *bz = (float *) bzArray->data;            // bz         float *bz = (float *) bzArray->data;            // bz
  
         // Get emphemeris         //Use magnetogram map to compute R
          DRMS_Segment_t *losSeg = drms_segment_lookup(inRec, "magnetogram");
          DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
          float *los = (float *) losArray->data;          // los
   
           DRMS_Segment_t *bz_errSeg = drms_segment_lookup(inRec, BR_ERR_SEG_CEA);
           DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
           float *bz_err = (float *) bz_errArray->data;            // bz_err
   
           DRMS_Segment_t *by_errSeg = drms_segment_lookup(inRec, BT_ERR_SEG_CEA);
           DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
           float *by_err = (float *) by_errArray->data;            // by_err
           //for (int i = 0; i < nxny; i++) by_err[i] *= -1;
   
           DRMS_Segment_t *bx_errSeg = drms_segment_lookup(inRec, BP_ERR_SEG_CEA);
           DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
           float *bx_err = (float *) bx_errArray->data;            // bx_err
  
         //float cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);          // Get emphemeris
         float  cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);          float  cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
         float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);         float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
         double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);         double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);
         double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);         double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);
Line 1698  void computeSWIndex(struct swIndex *swKe
Line 1943  void computeSWIndex(struct swIndex *swKe
         float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);         float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);
         float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);         float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);
  
         //float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)          // convert cdelt1_orig from degrees to arcsec
           float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
  
         printf("cdelt1=%f\n",cdelt1);          //if (nx1 > floor((nx-1)/scale + 1) )
         printf("rsun_ref=%f\n",rsun_ref);          //      DIE("X-dimension of output array in fsample() is too large.");
         printf("rsun_obs=%f\n",rsun_obs);          //if (ny1 > floor((ny-1)/scale + 1) )
         printf("dsun_obs=%f\n",dsun_obs);          //      DIE("Y-dimension of output array in fsample() is too large.");
  
         // Temp arrays         // Temp arrays
   
         float *bh = (float *) (malloc(nxny * sizeof(float)));         float *bh = (float *) (malloc(nxny * sizeof(float)));
         float *bt = (float *) (malloc(nxny * sizeof(float)));         float *bt = (float *) (malloc(nxny * sizeof(float)));
         float *jz = (float *) (malloc(nxny * sizeof(float)));         float *jz = (float *) (malloc(nxny * sizeof(float)));
Line 1722  void computeSWIndex(struct swIndex *swKe
Line 1967  void computeSWIndex(struct swIndex *swKe
         float *dery_bh = (float *) (malloc(nxny * sizeof(float)));         float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
         float *derx_bz = (float *) (malloc(nxny * sizeof(float)));         float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
         float *dery_bz = (float *) (malloc(nxny * sizeof(float)));         float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
           float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
         // The Computations          float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
           float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
         if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),          float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
                                            mask, bitmask, cdelt1, rsun_ref, rsun_obs)){          float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
           float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
   
           // define some values for the R calculation
           int scale = round(2.0/cdelt1);
           int nx1 = nx/scale;
           int ny1 = ny/scale;
           int nxp = nx1+40;
           int nyp = ny1+40;
           float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1n0    = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1p     = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1n     = (float *)malloc(nx1*ny1*sizeof(float));
           float *p1      = (float *)malloc(nx1*ny1*sizeof(float));
           float *pmap    = (float *)malloc(nxp*nyp*sizeof(float));
           float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
           float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
   
           // define some arrays for the lorentz force calculation
           float *fx = (float *) (malloc(nxny * sizeof(float)));
           float *fy = (float *) (malloc(nxny * sizeof(float)));
           float *fz = (float *) (malloc(nxny * sizeof(float)));
   
   
           //spaceweather quantities computed
           if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),  &(swKeys_ptr->mean_vf_err),
                              &(swKeys_ptr->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
           {
                 swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;               // If fail, fill in NaN                 swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
                 swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
                   swKeys_ptr->mean_vf_err = DRMS_MISSING_FLOAT;
                   swKeys_ptr->count_mask  = DRMS_MISSING_INT;
         }         }
  
         for (int i = 0; i < nxny; i++) bpz[i] = bz[i];         for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
         greenpot(bpx, bpy, bpz, nx, ny);         greenpot(bpx, bpy, bpz, nx, ny);
  
         computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);          computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
  
         if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask, bitmask))          if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), &(swKeys_ptr->mean_gamma_err),mask, bitmask))
           {
                 swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
                   swKeys_ptr->mean_gamma_err =  DRMS_MISSING_FLOAT;
           }
  
         computeB_total(bx, by, bz, bt, dims, mask, bitmask);          computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);
  
         if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt))          if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt,
                                       dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err)))
           {
                 swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
                   swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
           }
  
         if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, bitmask, derx_bh, dery_bh))          if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh),
                                   &(swKeys_ptr->mean_derivative_bh_err), mask, bitmask, derx_bh, dery_bh))
           {
                 swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
                   swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT;
           }
  
         if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, bitmask, derx_bz, dery_bz))          if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err),
                                   mask, bitmask, derx_bz, dery_bz))
           {
                 swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN                 swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
                   swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT;
           }
  
         computeJz(bx, by, dims, jz, mask, bitmask, cdelt1, rsun_ref, rsun_obs, derx, dery);          computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
                 derx, dery);
  
         struct fresize_struct convolution_array;  
         init_fresize_gaussian(&convolution_array, 4.0f, 12, 1);  
         fresize(&convolution_array, jz, jz_smooth, nx, ny, nx, nx, ny, nx, 0, 0, 0.0f);  
  
         if(computeJzsmooth(bx, by, dims, jz_smooth, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask,          if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &(swKeys_ptr->mean_jz),
                                  cdelt1, rsun_ref, rsun_obs, derx, dery))                         &(swKeys_ptr->mean_jz_err), &(swKeys_ptr->us_i), &(swKeys_ptr->us_i_err), mask, bitmask, cdelt1,
                          rsun_ref, rsun_obs, derx, dery))
         {         {
                 swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
                 swKeys_ptr->us_i = DRMS_MISSING_FLOAT;                 swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
                   swKeys_ptr->mean_jz_err        = DRMS_MISSING_FLOAT;
                   swKeys_ptr->us_i_err           = DRMS_MISSING_FLOAT;
         }         }
         printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);  
   
  
         if (computeAlpha(bz, dims, jz_smooth, &(swKeys_ptr->mean_alpha), mask, bitmask, cdelt1, rsun_ref, rsun_obs))          if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &(swKeys_ptr->mean_alpha), &(swKeys_ptr->mean_alpha_err),
                            mask, bitmask, cdelt1, rsun_ref, rsun_obs))
           {
                 swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
                   swKeys_ptr->mean_alpha_err     = DRMS_MISSING_FLOAT;
           }
  
         if (computeHelicity(bz, dims, jz_smooth, &(swKeys_ptr->mean_ih),          if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &(swKeys_ptr->mean_ih), &(swKeys_ptr->mean_ih_err),
                                                 &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),                                                 &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
                                                 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {                              &(swKeys_ptr->total_us_ih_err), &(swKeys_ptr->total_abs_ih_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
           {
                 swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
                 swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;                 swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
                 swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;                 swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
                   swKeys_ptr->mean_ih_err        = DRMS_MISSING_FLOAT;
                   swKeys_ptr->total_us_ih_err    = DRMS_MISSING_FLOAT;
                   swKeys_ptr->total_abs_ih_err   = DRMS_MISSING_FLOAT;
         }         }
  
         if (computeSumAbsPerPolarity(bz, jz_smooth, dims, &(swKeys_ptr->totaljz),          if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &(swKeys_ptr->totaljz), &(swKeys_ptr->totaljz_err),
                                                                  mask, bitmask, cdelt1, rsun_ref, rsun_obs))                                                                  mask, bitmask, cdelt1, rsun_ref, rsun_obs))
           {
                 swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;                 swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
                   swKeys_ptr->totaljz_err        = DRMS_MISSING_FLOAT;
           }
  
           if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims,
         if (computeFreeEnergy(bx, by, bpx, bpy, dims,                                &(swKeys_ptr->meanpot), &(swKeys_ptr->meanpot_err), &(swKeys_ptr->totpot), &(swKeys_ptr->totpot_err),
                                                   &(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),                                mask, bitmask, cdelt1, rsun_ref, rsun_obs))
                                                   mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {          {
                 swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN                 swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN
                 swKeys_ptr->totpot = DRMS_MISSING_FLOAT;                 swKeys_ptr->totpot = DRMS_MISSING_FLOAT;
                   swKeys_ptr->meanpot_err        = DRMS_MISSING_FLOAT;
                   swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;
         }         }
  
         if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,  
                                                   &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),          if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,
                                                   &(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),                                &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),
                                                   mask, bitmask)) {                                mask, bitmask))
           {
                 swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN                 swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
                 swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;                 swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
                 swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN                  swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;
                 swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;          }
   
           if (computeR(bz_err, los , dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,
                        p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
           {
                   swKeys_ptr->Rparam = DRMS_MISSING_FLOAT;                // If fail, fill in NaN
         }         }
  
         // Clean up  
           if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &(swKeys_ptr->totfx), &(swKeys_ptr->totfy), &(swKeys_ptr->totfz), &(swKeys_ptr->totbsq),
              &(swKeys_ptr->epsx), &(swKeys_ptr->epsy), &(swKeys_ptr->epsz), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
           {
                   swKeys_ptr->totfx             = DRMS_MISSING_FLOAT;
                   swKeys_ptr->totfy             = DRMS_MISSING_FLOAT;
                   swKeys_ptr->totfz             = DRMS_MISSING_FLOAT;
                   swKeys_ptr->totbsq            = DRMS_MISSING_FLOAT;
                   swKeys_ptr->epsx              = DRMS_MISSING_FLOAT;
                   swKeys_ptr->epsy              = DRMS_MISSING_FLOAT;
                   swKeys_ptr->epsz              = DRMS_MISSING_FLOAT;
   
           }
   
   
           // Clean up the arrays
  
         drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong         drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong
         drms_free_array(maskArray);         drms_free_array(maskArray);
         drms_free_array(bxArray);         drms_free_array(bxArray);
         drms_free_array(byArray);         drms_free_array(byArray);
         drms_free_array(bzArray);         drms_free_array(bzArray);
           drms_free_array(losArray);              // Mar 7
           drms_free_array(bx_errArray);
           drms_free_array(by_errArray);
           drms_free_array(bz_errArray);
  
         free(bh); free(bt); free(jz); free(jz_smooth);         free(bh); free(bt); free(jz); free(jz_smooth);
         free(bpx); free(bpy); free(bpz);         free(bpx); free(bpy); free(bpz);
Line 1812  void computeSWIndex(struct swIndex *swKe
Line 2139  void computeSWIndex(struct swIndex *swKe
         free(derx_bt); free(dery_bt);         free(derx_bt); free(dery_bt);
         free(derx_bz); free(dery_bz);         free(derx_bz); free(dery_bz);
         free(derx_bh); free(dery_bh);         free(derx_bh); free(dery_bh);
           free(bt_err); free(bh_err);  free(jz_err);
           free(jz_err_squared); free(jz_rms_err);
           free(jz_err_squared_smooth);
   
           // free the arrays that are related to the r calculation
           free(rim);
           free(p1p0);
           free(p1n0);
           free(p1p);
           free(p1n);
           free(p1);
           free(pmap);
           free(p1pad);
           free(pmapn);
  
         free_fresize(&convolution_array);          // free the arrays that are related to the lorentz calculation
           free(fx); free(fy); free(fz);
 } }
  
   
 /* /*
  * Set space weather indices, no error checking for now  * Set space weather indices, no error checking for now
  *  *
Line 1840  void setSWIndex(DRMS_Record_t *outRec, s
Line 2181  void setSWIndex(DRMS_Record_t *outRec, s
         drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);         drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);
         drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);         drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
         drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);         drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
       drms_setkey_float(outRec, "CMASK",   swKeys_ptr->count_mask);
       drms_setkey_float(outRec, "ERRBT",   swKeys_ptr->mean_derivative_btotal_err);
       drms_setkey_float(outRec, "ERRVF",   swKeys_ptr->mean_vf_err);
       drms_setkey_float(outRec, "ERRGAM",  swKeys_ptr->mean_gamma_err);
       drms_setkey_float(outRec, "ERRBH",   swKeys_ptr->mean_derivative_bh_err);
       drms_setkey_float(outRec, "ERRBZ",   swKeys_ptr->mean_derivative_bz_err);
       drms_setkey_float(outRec, "ERRJZ",   swKeys_ptr->mean_jz_err);
       drms_setkey_float(outRec, "ERRUSI",  swKeys_ptr->us_i_err);
       drms_setkey_float(outRec, "ERRALP",  swKeys_ptr->mean_alpha_err);
       drms_setkey_float(outRec, "ERRMIH",  swKeys_ptr->mean_ih_err);
       drms_setkey_float(outRec, "ERRTUI",  swKeys_ptr->total_us_ih_err);
       drms_setkey_float(outRec, "ERRTAI",  swKeys_ptr->total_abs_ih_err);
       drms_setkey_float(outRec, "ERRJHT",  swKeys_ptr->totaljz_err);
       drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);
       drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);
       drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);
       drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);
       drms_setkey_float(outRec, "TOTFX",   swKeys_ptr->totfx);
       drms_setkey_float(outRec, "TOTFY",   swKeys_ptr->totfy);
       drms_setkey_float(outRec, "TOTFZ",   swKeys_ptr->totfz);
       drms_setkey_float(outRec, "TOTBSQ",  swKeys_ptr->totbsq);
       drms_setkey_float(outRec, "EPSX",    swKeys_ptr->epsx);
       drms_setkey_float(outRec, "EPSY",    swKeys_ptr->epsy);
       drms_setkey_float(outRec, "EPSZ",    swKeys_ptr->epsz);
 }; };
  
   
 /* /*
  * Set all keywords, no error checking for now  * Set all keywords, no error checking for now
  *  *
  */  */
  
 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo)  void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
 { {
         copy_me_keys(inRec, outRec);  
         copy_patch_keys(inRec, outRec);          copy_me_keys(bharpRec, outRec);
         copy_geo_keys(inRec, outRec);          copy_patch_keys(mharpRec, outRec);      // Dec 30
         copy_ambig_keys(inRec, outRec);          copy_geo_keys(mharpRec, outRec);        // Dec 30
           copy_ambig_keys(bharpRec, outRec);
  
         int status = 0;         int status = 0;
  
         // Change a few geometry keywords for CEA records          // Change a few geometry keywords for CEA & cutout records
         if (mInfo != NULL) {          if (mInfo != NULL) {        // CEA
  
           drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);           drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
                 drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);                 drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
Line 1877  void setKeys(DRMS_Record_t *outRec, DRMS
Line 2242  void setKeys(DRMS_Record_t *outRec, DRMS
                 drms_setkey_string(outRec, "CTYPE2", key);                 drms_setkey_string(outRec, "CTYPE2", key);
                 drms_setkey_float(outRec, "CROTA2", 0.0);                 drms_setkey_float(outRec, "CROTA2", 0.0);
  
         } else {          // Jan 2 2014 XS
           int nSeg = ARRLENGTH(CEASegs);
           for (int iSeg = 0; iSeg < nSeg; iSeg++) {
               DRMS_Segment_t *outSeg = NULL;
               outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);
               if (!outSeg) continue;
               // Set Bunit
               char bunit_xxx[20];
               sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
               //printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);
               drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);
           }
  
           float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status);          } else {        // Cutout
           float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status);  
           float x_ll = drms_getkey_float(inRec, "CRPIX1", &status);          float disk_xc, disk_yc;
           float y_ll = drms_getkey_float(inRec, "CRPIX2", &status);          if (fullDisk) {
               disk_xc = drms_getkey_float(bharpRec, "CRPIX1", &status);
               disk_yc = drms_getkey_float(bharpRec, "CRPIX2", &status);
           } else {
               disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);
               disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);
           }
           float x_ll = drms_getkey_float(mharpRec, "CRPIX1", &status);
           float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &status);
           // Defined as disk center's pixel address wrt lower-left of cutout           // Defined as disk center's pixel address wrt lower-left of cutout
           drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);           drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
                 drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);                 drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
Line 1890  void setKeys(DRMS_Record_t *outRec, DRMS
Line 2274  void setKeys(DRMS_Record_t *outRec, DRMS
                 drms_setkey_float(outRec, "CRVAL1", 0);                 drms_setkey_float(outRec, "CRVAL1", 0);
                 drms_setkey_float(outRec, "CRVAL2", 0);                 drms_setkey_float(outRec, "CRVAL2", 0);
  
           // Jan 2 2014 XS
           int nSeg = ARRLENGTH(CutSegs);
           for (int iSeg = 0; iSeg < nSeg; iSeg++) {
               DRMS_Segment_t *outSeg = NULL;
               outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);
               if (!outSeg) continue;
               // Set Bunit
               char bunit_xxx[20];
               sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
               //printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);
               drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);
           }
   
   
         }         }
  
         char timebuf[1024];          // Mar 19 XS
         float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */          if (fullDisk) {
         double val;                  drms_setkey_int(outRec, "AMBPATCH", 0);
         status = 0;                  drms_setkey_int(outRec, "AMBWEAK", 2);
           } else {
         val = drms_getkey_double(inRec, "DATE",&status);                  drms_setkey_int(outRec, "AMBPATCH", 1);
         drms_setkey_double(outRec, "DATE_B", val);          }
         sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);  
         drms_setkey_string(outRec, "DATE", timebuf);      TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
       tnow = (double)time(NULL);
       tnow += UNIX_epoch;
   
       val = drms_getkey_time(bharpRec, "DATE", &status);
       drms_setkey_time(outRec, "DATE_B", val);
       drms_setkey_time(outRec, "DATE", tnow);
  
         // set cvs commit version into keyword HEADER         // set cvs commit version into keyword HEADER
         char *cvsinfo = strdup("$Id$");         char *cvsinfo = strdup("$Id$");
         //   status = drms_setkey_string(outRec, "HEADER", cvsinfo);      char *cvsinfo2 = sw_functions_version();
         status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);      char cvsinfoall[2048];
       strcat(cvsinfoall,cvsinfo);
       strcat(cvsinfoall,"\n");
       strcat(cvsinfoall,cvsinfo2);
       status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
  
 }; };
  
 //  
 //  
  
 /* ############# Nearest neighbour interpolation ############### */ /* ############# Nearest neighbour interpolation ############### */
  
Line 1927  float nnb (float *f, int nx, int ny, dou
Line 2333  float nnb (float *f, int nx, int ny, dou
  
 } }
  
   
 /* ################## Wrapper for Jesper's rebin code ################## */ /* ################## Wrapper for Jesper's rebin code ################## */
  
 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss) void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)


Legend:
Removed from v.1.11  
changed lines
  Added in v.1.31

Karen Tian
Powered by
ViewCVS 0.9.4