(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.5 and 1.10

version 1.5, 2012/10/23 19:35:37 version 1.10, 2013/01/23 21:48:26
Line 20 
Line 20 
  *              v0.0    Jul 02 2012  *              v0.0    Jul 02 2012
  *              v0.1    Jul 23 2012  *              v0.1    Jul 23 2012
  *              v0.2    Sep 04 2012  *              v0.2    Sep 04 2012
    *              v0.3    Dec 18 2012
    *              v0.4    Jan 02 2013
    *              v0.5    Jam 23 2012
  *  *
  *      Notes:  *      Notes:
  *              v0.0  *              v0.0
Line 33 
Line 36 
  *              SW indices fixed  *              SW indices fixed
  *              Added doppler and continuum  *              Added doppler and continuum
  *              Added other keywords: HEADER (populated by cvs build version), DATE_B  *              Added other keywords: HEADER (populated by cvs build version), DATE_B
    *              v0.3
    *              Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"
    *              v0.4
    *              Took out convert_inplace(). Was causing all the images to be int
    *              v0.5
    *              Corrected ephemeris keywords, added argument mInfo for setKeys()
  *  *
  *      Example:  *      Example:
  *      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \  *      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
Line 149  enum projection {
Line 158  enum projection {
         lambert         lambert
 }; };
  
   // WSC code
   char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
           "SIN", "ZEA"};
   
 // Ephemeris // Ephemeris
 struct ephemeris { struct ephemeris {
         double disk_lonc, disk_latc;         double disk_lonc, disk_latc;
Line 247  void computeSWIndex(struct swIndex *swKe
Line 260  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);  void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo);
  
 // =================== // ===================
  
Line 294  char *CEASegs[] = {"magnetogram", "bitma
Line 307  char *CEASegs[] = {"magnetogram", "bitma
  
  
 char *module_name = "sharp"; char *module_name = "sharp";
 char *version_id = "2012 Jul 02";  /* Version number */  char *version_id = "2012 Dec 18";  /* Version number */
  
 ModuleArgs_t module_args[] = ModuleArgs_t module_args[] =
 { {
Line 405  int DoIt(void)
Line 418  int DoIt(void)
  
         } // irec         } // irec
  
   
         drms_close_records(mharpRS, DRMS_FREE_RECORD);         drms_close_records(mharpRS, DRMS_FREE_RECORD);
         drms_close_records(bharpRS, DRMS_FREE_RECORD);         drms_close_records(bharpRS, DRMS_FREE_RECORD);
           drms_close_records(dopRS, DRMS_FREE_RECORD);                            // Dec 18 2012
           drms_close_records(contRS, DRMS_FREE_RECORD);                           // Dec 18 2012
  
         return 0;         return 0;
  
Line 646  int createCeaRecord(DRMS_Record_t *mharp
Line 662  int createCeaRecord(DRMS_Record_t *mharp
         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);
  
         setKeys(sharpRec, bharpRec);            // Set all other keywords          setKeys(sharpRec, 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 730  int mapScaler(DRMS_Record_t *sharpRec, D
Line 746  int mapScaler(DRMS_Record_t *sharpRec, D
         outSeg = drms_segment_lookup(sharpRec, segName);         outSeg = drms_segment_lookup(sharpRec, segName);
         if (!outSeg) return 1;         if (!outSeg) return 1;
  
         DRMS_Type_t arrayType = outSeg->info->type;  //      DRMS_Type_t arrayType = outSeg->info->type;
         DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);         DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
         if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}         if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
  
         // convert to needed data type         // convert to needed data type
  
         drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);  //      drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);         // Jan 02 2013
  
         outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];         outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
         outArray->parent_segment = outSeg;  //      outArray->parent_segment = outSeg;
         outArray->israw = 0;            // always compressed         outArray->israw = 0;            // always compressed
         outArray->bzero = outSeg->bzero;         outArray->bzero = outSeg->bzero;
         outArray->bscale = outSeg->bscale;         outArray->bscale = outSeg->bscale;
Line 748  int mapScaler(DRMS_Record_t *sharpRec, D
Line 764  int mapScaler(DRMS_Record_t *sharpRec, D
         if (status) return 0;         if (status) return 0;
  
         if (inArray) drms_free_array(inArray);         if (inArray) drms_free_array(inArray);
           if ((xsz != FOURK) || (ysz != FOURK)) free(inData);                     // Dec 18 2012
         if (outArray) drms_free_array(outArray);         if (outArray) drms_free_array(outArray);
         return 0;         return 0;
  
Line 814  int mapVectorB(DRMS_Record_t *sharpRec,
Line 831  int mapVectorB(DRMS_Record_t *sharpRec,
                 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);                 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
                 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);                 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
                 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];                 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
                 outArray->parent_segment = outSeg;  //              outArray->parent_segment = outSeg;
                 outArray->israw = 0;                 outArray->israw = 0;
                 outArray->bzero = outSeg->bzero;                 outArray->bzero = outSeg->bzero;
                 outArray->bscale = outSeg->bscale;                 outArray->bscale = outSeg->bscale;
Line 866  int mapVectorBErr(DRMS_Record_t *sharpRe
Line 883  int mapVectorBErr(DRMS_Record_t *sharpRe
                 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);                 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
                 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);                 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
                 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];                 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
                 outArray->parent_segment = outSeg;  //              outArray->parent_segment = outSeg;
                 outArray->israw = 0;                 outArray->israw = 0;
                 outArray->bzero = outSeg->bzero;                 outArray->bzero = outSeg->bzero;
                 outArray->bscale = outSeg->bscale;                 outArray->bscale = outSeg->bscale;
Line 1107  int performSampling(float *outData, floa
Line 1124  int performSampling(float *outData, floa
         // Rebinning, smoothing         // Rebinning, smoothing
  
         frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian         frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
           free(outData0);         // Dec 18 2012
  
         //         //
  
Line 1528  int createCutRecord(DRMS_Record_t *mharp
Line 1546  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);            // Set all other keywords          setKeys(sharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout
  
         // Stats         // Stats
  
Line 1589  int writeCutout(DRMS_Record_t *outRec, D
Line 1607  int writeCutout(DRMS_Record_t *outRec, D
  
 #if DISAMB_AZI #if DISAMB_AZI
         if (!strcmp(SegName, "azimuth")) {         if (!strcmp(SegName, "azimuth")) {
                 DRMS_Segment_t *disambSeg = drms_segment_lookup(inRec, "disambig");                  DRMS_Segment_t *disambSeg = NULL;
                   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 (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;}
                 } else {                 } else {
                     drms_free_array(cutoutArray);
                         return 1;                         return 1;
                 }                 }
                 double *azimuth = (double *) cutoutArray->data;                 double *azimuth = (double *) cutoutArray->data;
Line 1611  int writeCutout(DRMS_Record_t *outRec, D
Line 1631  int writeCutout(DRMS_Record_t *outRec, D
  
         outSeg = drms_segment_lookup(outRec, SegName);         outSeg = drms_segment_lookup(outRec, SegName);
         if (!outSeg) return 1;         if (!outSeg) return 1;
         drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);  //      drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);      // Jan 02 2013
         outSeg->axis[0] = cutoutArray->axis[0];         outSeg->axis[0] = cutoutArray->axis[0];
         outSeg->axis[1] = cutoutArray->axis[1];         outSeg->axis[1] = cutoutArray->axis[1];
         cutoutArray->parent_segment = outSeg;  //      cutoutArray->parent_segment = outSeg;
         cutoutArray->israw = 0;         // always compressed         cutoutArray->israw = 0;         // always compressed
     cutoutArray->bzero = outSeg->bzero;     cutoutArray->bzero = outSeg->bzero;
     cutoutArray->bscale = outSeg->bscale;               // Same as inArray's     cutoutArray->bscale = outSeg->bscale;               // Same as inArray's
Line 1690  void computeSWIndex(struct swIndex *swKe
Line 1710  void computeSWIndex(struct swIndex *swKe
         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)));
           float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
         float *bpx = (float *) (malloc(nxny * sizeof(float)));         float *bpx = (float *) (malloc(nxny * sizeof(float)));
         float *bpy = (float *) (malloc(nxny * sizeof(float)));         float *bpy = (float *) (malloc(nxny * sizeof(float)));
         float *bpz = (float *) (malloc(nxny * sizeof(float)));         float *bpz = (float *) (malloc(nxny * sizeof(float)));
Line 1702  void computeSWIndex(struct swIndex *swKe
Line 1723  void computeSWIndex(struct swIndex *swKe
         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)));
  
         // Compute          // The Computations
  
         if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),         if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
                                            mask, bitmask, cdelt1, rsun_ref, rsun_obs)){                                            mask, bitmask, cdelt1, rsun_ref, rsun_obs)){
Line 1729  void computeSWIndex(struct swIndex *swKe
Line 1750  void computeSWIndex(struct swIndex *swKe
         if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, bitmask, derx_bz, dery_bz))         if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), 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
  
           computeJz(bx, by, dims, jz, 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(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask,          if(computeJzsmooth(bx, by, dims, jz_smooth, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask,
                      cdelt1, rsun_ref, rsun_obs, derx, dery)) {                                   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;
         }         }
   
                         printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);                         printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
  
         if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, bitmask, cdelt1, rsun_ref, rsun_obs))  
           if (computeAlpha(bz, dims, jz_smooth, &(swKeys_ptr->mean_alpha), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
                 swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
  
         if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih),          if (computeHelicity(bz, dims, jz_smooth, &(swKeys_ptr->mean_ih),
                                                 &(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)) {                                                 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
                 swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
Line 1750  void computeSWIndex(struct swIndex *swKe
Line 1776  void computeSWIndex(struct swIndex *swKe
                 swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;                 swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
         }         }
  
         if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),          if (computeSumAbsPerPolarity(bz, jz_smooth, dims, &(swKeys_ptr->totaljz),
                                                                  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;
  
Line 1774  void computeSWIndex(struct swIndex *swKe
Line 1800  void computeSWIndex(struct swIndex *swKe
  
         // Clean up         // Clean up
  
           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);
  
         free(bh); free(bt); free(jz);          free(bh); free(bt); free(jz); free(jz_smooth);
         free(bpx); free(bpy); free(bpz);         free(bpx); free(bpy); free(bpz);
         free(derx); free(dery);         free(derx); free(dery);
         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_fresize(&convolution_array);
 } }
  
  
Line 1820  void setSWIndex(DRMS_Record_t *outRec, s
Line 1848  void setSWIndex(DRMS_Record_t *outRec, s
  *  *
  */  */
  
 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)  void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo)
 { {
    copy_me_keys(inRec, outRec);    copy_me_keys(inRec, outRec);
    copy_patch_keys(inRec, outRec);    copy_patch_keys(inRec, outRec);
    copy_geo_keys(inRec, outRec);    copy_geo_keys(inRec, outRec);
    copy_ambig_keys(inRec, outRec);    copy_ambig_keys(inRec, outRec);
  
           int status = 0;
   
           // Change a few geometry keywords for CEA records
           if (mInfo != NULL) {
   
             drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
                   drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
   
                   drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
                   drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
                   drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
                   drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
                   drms_setkey_string(outRec, "CUNIT1", "degree");
                   drms_setkey_string(outRec, "CUNIT2", "degree");
   
                   char key[64];
                   snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
                   drms_setkey_string(outRec, "CTYPE1", key);
                   snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
                   drms_setkey_string(outRec, "CTYPE2", key);
                   drms_setkey_float(outRec, "CROTA2", 0.0);
   
           } else {
   
             float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status);
             float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status);
             float x_ll = drms_getkey_float(inRec, "CRPIX1", &status);
             float y_ll = drms_getkey_float(inRec, "CRPIX2", &status);
             // 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, "CRPIX2", disk_yc - y_ll + 1.);
                   // Always 0.
                   drms_setkey_float(outRec, "CRVAL1", 0);
                   drms_setkey_float(outRec, "CRVAL2", 0);
   
           }
   
    char timebuf[1024];    char timebuf[1024];
    float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */    float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
    double val;    double val;
    int status = DRMS_SUCCESS;          status = 0;
  
    val = drms_getkey_double(inRec, "DATE",&status);    val = drms_getkey_double(inRec, "DATE",&status);
    drms_setkey_double(outRec, "DATE_B", val);    drms_setkey_double(outRec, "DATE_B", val);


Legend:
Removed from v.1.5  
changed lines
  Added in v.1.10

Karen Tian
Powered by
ViewCVS 0.9.4