(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.31

version 1.20, 2014/02/13 04:41:06 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 26 
Line 27 
  *              v0.6    Aug 12 2013  *              v0.6    Aug 12 2013
  *      v0.7    Jan 02 2014  *      v0.7    Jan 02 2014
  *              v0.8    Feb 12 2014  *              v0.8    Feb 12 2014
    *                              v0.9 Mar 04 2014
  *  *
  *      Notes:  *      Notes:
  *              v0.0  *              v0.0
Line 58 
Line 60 
  *              Added disambig to azimuth during error propagation  *              Added disambig to azimuth during error propagation
  *              Changed usage for disambig: bit 2 (radial acute) for full disk, bit 0 for patch  *              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  *              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:   *      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 179  struct swIndex {
     float meanpot_err;     float meanpot_err;
     float totpot_err;     float totpot_err;
     float meanshear_angle_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 346  char *CutBunits[] = {"Mx/cm^2", " ", "cm
Line 355  char *CutBunits[] = {"Mx/cm^2", " ", "cm
     " ", " ", " ",     " ", " ", " ",
     " ", " ", " ",     " ", " ", " ",
     " ", " "};     " ", " "};
 char *CEABunits[] = {"Mx/cm^2", " ", "cm/s", "DN/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"};      "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", " "};      // Mar 4 2014 XS
  
 /* ========================================================================================================== */ /* ========================================================================================================== */
  
Line 355  char *module_name = "sharp";
Line 364  char *module_name = "sharp";
 int seed; int seed;
  
 int fullDisk;       // full disk mode 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[] =
 { {
Line 366  ModuleArgs_t module_args[] =
Line 376  ModuleArgs_t module_args[] =
         {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,    "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}
 }; };
  
Line 395  int DoIt(void)
Line 406  int DoIt(void)
         sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");         sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
  
     seed = params_get_int(&cmdparams, "seed");     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 */
  
Line 409  int DoIt(void)
Line 421  int DoIt(void)
         SHOW("Harp mode\n");         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     // 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.");
Line 1029  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 1064  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 1419  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 1595  int readVectorBErr(DRMS_Record_t *inRec,
Line 1617  int readVectorBErr(DRMS_Record_t *inRec,
  
         // Add disambig, Feb 12 2014         // Add disambig, Feb 12 2014
  
         DRMS_Segment_t *inSeg = drms_segment_lookup(inRec, "disambig");          DRMS_Segment_t *inSeg;
         DRMS_Array_t *inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);      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;         if (status) return 1;
         char *ambig = (char *)inArray_ambig->data;         char *ambig = (char *)inArray_ambig->data;
  
Line 1626  int readVectorBErr(DRMS_Record_t *inRec,
Line 1653  int readVectorBErr(DRMS_Record_t *inRec,
         }         }
   }   }
  
       }
   
         // 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 1650  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]);
         drms_free_array(inArray_ambig);         // Feb 12 2014          if (amb4err) drms_free_array(inArray_ambig);            // Feb 12; Mar 04 2014
  
         return 0;         return 0;
  
Line 1886  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
  
          //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 1912  void computeSWIndex(struct swIndex *swKe
Line 1946  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.);
  
           //if (nx1 > floor((nx-1)/scale + 1) )
           //      DIE("X-dimension of output array in fsample() is too large.");
           //if (ny1 > floor((ny-1)/scale + 1) )
           //      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)));
Line 1934  void computeSWIndex(struct swIndex *swKe
Line 1973  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  
  
           // 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),         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 1959  void computeSWIndex(struct swIndex *swKe
Line 2019  void computeSWIndex(struct swIndex *swKe
  
         computeB_total(bx_err, by_err, bz_err, bt_err, 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, bt_err, &(swKeys_ptr->mean_derivative_btotal_err)))          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;                 swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
     }     }
  
         if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh), &(swKeys_ptr->mean_derivative_bh_err), 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;         swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT;
         }         }
  
         if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err), 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;         swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT;
Line 1991  void computeSWIndex(struct swIndex *swKe
Line 2054  void computeSWIndex(struct swIndex *swKe
         swKeys_ptr->us_i_err           = DRMS_MISSING_FLOAT;         swKeys_ptr->us_i_err           = DRMS_MISSING_FLOAT;
         }         }
  
         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))          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;         swKeys_ptr->mean_alpha_err     = DRMS_MISSING_FLOAT;
     }     }
  
         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),          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_err), &(swKeys_ptr->total_abs_ih_err), 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;
Line 2025  void computeSWIndex(struct swIndex *swKe
Line 2090  void computeSWIndex(struct swIndex *swKe
         swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;         swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;
         }         }
  
   
         if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,         if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,
                                                   &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),                                                   &(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_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, scale, p1pad, nxp, nyp, pmapn))
           {
                   swKeys_ptr->Rparam = DRMS_MISSING_FLOAT;                // If fail, fill in NaN
           }
   
   
           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         // Clean up the arrays
  
         drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong         drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong
Line 2040  void computeSWIndex(struct swIndex *swKe
Line 2128  void computeSWIndex(struct swIndex *swKe
         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 2050  void computeSWIndex(struct swIndex *swKe
Line 2142  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 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 the arrays that are related to the lorentz calculation
           free(fx); free(fy); free(fz);
 } }
  
 /* /*
Line 2091  void setSWIndex(DRMS_Record_t *outRec, s
Line 2197  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);
       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);
 }; };
  
 /* /*
Line 2176  void setKeys(DRMS_Record_t *outRec, DRMS
Line 2290  void setKeys(DRMS_Record_t *outRec, DRMS
  
         }         }
  
         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);
                   drms_setkey_int(outRec, "AMBWEAK", 2);
         val = drms_getkey_double(bharpRec, "DATE", &status);          } else {
         drms_setkey_double(outRec, "DATE_B", val);                  drms_setkey_int(outRec, "AMBPATCH", 1);
         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$");
Line 2212  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.20  
changed lines
  Added in v.1.31

Karen Tian
Powered by
ViewCVS 0.9.4