(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.20 and 1.21

version 1.20, 2014/02/13 04:41:06 version 1.21, 2014/02/18 23:35:31
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 60 
Line 61 
  *              Fixed disambig cutout for patch: 0 for even, 7 for odd  *              Fixed disambig cutout for patch: 0 for even, 7 for odd
  *  *
  *  *
  *      Example:   *      Example Calls:
   B (full disk disambiguation)   *      [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"   *      > 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"
  BHARP (patch disambiguation)   *      [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"   *      > 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"
  *  
  *      For comparison:  
  *      bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \  
  "out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"  
  *  *
  *  *
  */  */
Line 178  struct swIndex {
Line 175  struct swIndex {
     float meanpot_err;     float meanpot_err;
     float totpot_err;     float totpot_err;
     float meanshear_angle_err;     float meanshear_angle_err;
           float Rparam;
 }; };
  
 // Mapping method // Mapping method
Line 1886  void computeSWIndex(struct swIndex *swKe
Line 1884  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
  
           //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_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);         DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
         float *bz_err = (float *) bz_errArray->data;            // bz_err         float *bz_err = (float *) bz_errArray->data;            // bz_err
Line 1911  void computeSWIndex(struct swIndex *swKe
Line 1914  void computeSWIndex(struct swIndex *swKe
  
         // convert cdelt1_orig from degrees to arcsec         // convert cdelt1_orig from degrees to arcsec
         float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);         float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
           int nx1 = nx*cdelt1/2;
           int ny1 = ny*cdelt1/2;
  
         // Temp arrays         // Temp arrays
         float *bh      = (float *) (malloc(nxny * sizeof(float)));         float *bh      = (float *) (malloc(nxny * sizeof(float)));
Line 1934  void computeSWIndex(struct swIndex *swKe
Line 1939  void computeSWIndex(struct swIndex *swKe
         float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));         float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
         float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));         float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
         float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));         float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
         //spaceweather quantities computed  
  
           // malloc some arrays for the R parameter calculation
           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(nx1*ny1*sizeof(float));
  
           //spaceweather quantities computed
         if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),  &(swKeys_ptr->mean_vf_err),         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->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
     {     {
Line 2033  void computeSWIndex(struct swIndex *swKe
Line 2046  void computeSWIndex(struct swIndex *swKe
         swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;         swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;
         }         }
  
           if (computeR(bz_err, los , dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,
                        p1p, p1n, p1, pmap, nx1, ny1))
           {
                   swKeys_ptr->Rparam = DRMS_MISSING_FLOAT;                // If fail, fill in NaN
           }
   
   
         // Clean up the arrays         // Clean up the arrays
  
         drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong         drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong
Line 2050  void computeSWIndex(struct swIndex *swKe
Line 2070  void computeSWIndex(struct swIndex *swKe
         free(bt_err); free(bh_err);  free(jz_err);         free(bt_err); free(bh_err);  free(jz_err);
     free(jz_err_squared); free(jz_rms_err);     free(jz_err_squared); free(jz_rms_err);
     free(jz_err_squared_smooth);     free(jz_err_squared_smooth);
   
           free(rim);
           free(p1p0);
           free(p1n0);
           free(p1p);
           free(p1n);
           free(p1);
           free(pmap);
   
 } }
  
 /* /*
Line 2091  void setSWIndex(DRMS_Record_t *outRec, s
Line 2120  void setSWIndex(DRMS_Record_t *outRec, s
     drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);     drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);
     drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);     drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);
     drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);     drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);
           drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);
 }; };
  
 /* /*
Line 2176  void setKeys(DRMS_Record_t *outRec, DRMS
Line 2206  void setKeys(DRMS_Record_t *outRec, DRMS
  
         }         }
  
         char timebuf[1024];         TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
         float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */         tnow = (double)time(NULL);
         double val;         tnow += UNIX_epoch;
   
         val = drms_getkey_double(bharpRec, "DATE", &status);         val = drms_getkey_time(bharpRec, "DATE", &status);
         drms_setkey_double(outRec, "DATE_B", val);         drms_setkey_time(outRec, "DATE_B", val);
         sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);         drms_setkey_time(outRec, "DATE", tnow);
         drms_setkey_string(outRec, "DATE", timebuf);  
  
         // set cvs commit version into keyword HEADER         // set cvs commit version into keyword HEADER
         char *cvsinfo = strdup("$Id$");         char *cvsinfo = strdup("$Id$");

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