(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.1 and 1.5

version 1.1, 2012/08/23 08:38:33 version 1.5, 2012/10/23 19:35:37
Line 19 
Line 19 
  *      Version:  *      Version:
  *              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
  *  *
  *      Notes:  *      Notes:
  *              v0.0  *              v0.0
Line 31 
Line 32 
  *              Fixed char I/O thanks to Art  *              Fixed char I/O thanks to Art
  *              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
  *  *
  *      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 564  int createCeaRecord(DRMS_Record_t *mharp
Line 566  int createCeaRecord(DRMS_Record_t *mharp
  
         // Get ephemeris         // Get ephemeris
  
         if (getEphemeris(mharpRec, &(mInfo.ephem))) return 1;          if (getEphemeris(mharpRec, &(mInfo.ephem))) {
                   SHOW("CEA: get ephemeris error\n");
                   return 1;
           }
  
         // Find position         // Find position
  
         if (findPosition(mharpRec, &mInfo)) return 1;          if (findPosition(mharpRec, &mInfo)) {
                   SHOW("CEA: find position error\n");
                   return 1;
           }
  
         // 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
Line 584  int createCeaRecord(DRMS_Record_t *mharp
Line 592  int createCeaRecord(DRMS_Record_t *mharp
  
         // Mapping single segment: Mharp, etc.         // Mapping single segment: Mharp, etc.
  
         if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) return 1;          if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
         if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) return 1;                  SHOW("CEA: mapping magnetogram error\n");
                   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")) return 1;          if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
                   SHOW("CEA: mapping dopplergram error\n");
                   return 1;
           }
         printf("Dopplergram mapping done.\n");         printf("Dopplergram mapping done.\n");
  
         if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) return 1;          if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
                   SHOW("CEA: mapping continuum error\n");
                   return 1;
           }
         printf("Intensitygram mapping done.\n");         printf("Intensitygram mapping done.\n");
  
         //if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) return 1;          if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
         //printf("Conf disambig mapping done.\n");                  SHOW("CEA: mapping conf_disambig error\n");
                   return 1;
         if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "disambig")) return 1;          }
         printf("disambig mapping done.\n");          printf("Conf disambig mapping done.\n");
  
         // Mapping vector B         // Mapping vector B
  
         if (mapVectorB(sharpRec, bharpRec, &mInfo)) return 1;          if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
                   SHOW("CEA: mapping vector B error\n");
                   return 1;
           }
         printf("Vector B mapping done.\n");         printf("Vector B mapping done.\n");
  
         // Mapping vector B errors         // Mapping vector B errors
  
         if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) return 1;          if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) {
                   SHOW("CEA: mapping vector B uncertainty error\n");
                   return 1;
           }
         printf("Vector B error done.\n");         printf("Vector B error done.\n");
  
         // Keywords & Links         // Keywords & Links
Line 621  int createCeaRecord(DRMS_Record_t *mharp
Line 647  int createCeaRecord(DRMS_Record_t *mharp
         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);            // Set all other keywords
           drms_copykey(sharpRec, mharpRec, "QUALITY");            // copied from los records
  
         // Space weather         // Space weather
  
Line 1457  int createCutRecord(DRMS_Record_t *mharp
Line 1484  int createCutRecord(DRMS_Record_t *mharp
                         break;                         break;
                 }                 }
         }         }
         if (iHarpSeg != nMharpSegs) return 1;           // if failed          if (iHarpSeg != nMharpSegs) {
                   SHOW("Cutout: segment number unmatch\n");
                   return 1;               // if failed
           }
         printf("Magnetogram cutout done.\n");         printf("Magnetogram cutout done.\n");
  
         // Cutout Doppler         // Cutout Doppler
Line 1614  void computeSWIndex(struct swIndex *swKe
Line 1644  void computeSWIndex(struct swIndex *swKe
         // 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
         //DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "bitmap");          DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
         //DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);          DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
         //int *mask = (int *) maskArray->data;          // get the previously made mask array          int *bitmask = (int *) bitmaskArray->data;              // get the previously made mask array
  
         //Use conf_disambig map as a threshold on spaceweather quantities         //Use conf_disambig map as a threshold on spaceweather quantities
         //DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");          DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
         DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "disambig");  
         DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);         DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
         int *mask = (int *) maskArray->data;            // get the previously made mask array         int *mask = (int *) maskArray->data;            // get the previously made mask array
  
Line 1676  void computeSWIndex(struct swIndex *swKe
Line 1705  void computeSWIndex(struct swIndex *swKe
         // Compute         // Compute
  
         if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),         if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
                                            mask, cdelt1, rsun_ref, rsun_obs)){                                             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;
         }         }
Line 1684  void computeSWIndex(struct swIndex *swKe
Line 1713  void computeSWIndex(struct swIndex *swKe
         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);          computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
  
         if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask))          if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask, bitmask))
                 swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
  
         computeB_total(bx, by, bz, bt, dims, mask);          computeB_total(bx, by, bz, bt, dims, mask, bitmask);
  
         if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt))          if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt))
                 swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
  
         if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh))          if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, bitmask, derx_bh, dery_bh))
                 swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;                 swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
  
         if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, 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
  
  
  
         if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask,          if(computeJz(bx, by, dims, jz, &(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;
Line 1710  void computeSWIndex(struct swIndex *swKe
Line 1739  void computeSWIndex(struct swIndex *swKe
  
                         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, cdelt1, rsun_ref, rsun_obs))          if (computeAlpha(bz, dims, jz, &(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, &(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, 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;
                 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;
         }         }
  
         if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),         if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),
                                                                  mask, cdelt1, rsun_ref, rsun_obs))                                                                   mask, bitmask, cdelt1, rsun_ref, rsun_obs))
                 swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;                 swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
  
  
         if (computeFreeEnergy(bx, by, bpx, bpy, dims,         if (computeFreeEnergy(bx, by, bpx, bpy, dims,
                                                   &(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),                                                   &(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
                                                   mask, 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;
         }         }
Line 1736  void computeSWIndex(struct swIndex *swKe
Line 1765  void computeSWIndex(struct swIndex *swKe
         if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,         if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,
                                                   &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),                                                   &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),
                                                   &(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),                                                   &(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
                                                   mask)) {                                                    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_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
Line 1797  void setKeys(DRMS_Record_t *outRec, DRMS
Line 1826  void setKeys(DRMS_Record_t *outRec, DRMS
    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);
   
      char timebuf[1024];
      float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
      double val;
      int status = DRMS_SUCCESS;
   
      val = drms_getkey_double(inRec, "DATE",&status);
      drms_setkey_double(outRec, "DATE_B", val);
      sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
      drms_setkey_string(outRec, "DATE", timebuf);
   
      // set cvs commit version into keyword HEADER
      char *cvsinfo = strdup("$Header$");
   //   status = drms_setkey_string(outRec, "HEADER", cvsinfo);
           status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);
   
 }; };
  
 // //

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
Karen Tian
 Powered by
ViewCVS 0.9.4