(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.32 and 1.39

version 1.32, 2015/01/23 19:21:09 version 1.39, 2021/05/24 22:17:37
Line 110 
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)
  
   // Maximum variation of LONDTMAX-LONDTMIN
   #define MAXLONDIFF      (1.2e-4)
   
 // Some other things // Some other things
 #ifndef MIN #ifndef MIN
 #define MIN(a,b) (((a)<(b)) ? (a) : (b)) #define MIN(a,b) (((a)<(b)) ? (a) : (b))
Line 144 
Line 147 
 // Space weather keywords // Space weather keywords
 struct swIndex { struct swIndex {
     float mean_vf;     float mean_vf;
       float mean_vf_los;
     float count_mask;     float count_mask;
       float count_mask_los;
     float absFlux;     float absFlux;
       float absFlux_los;
     float mean_hf;     float mean_hf;
     float mean_gamma;     float mean_gamma;
     float mean_derivative_btotal;     float mean_derivative_btotal;
     float mean_derivative_bh;     float mean_derivative_bh;
     float mean_derivative_bz;     float mean_derivative_bz;
       float mean_derivative_los;
     float mean_jz;     float mean_jz;
     float us_i;     float us_i;
     float mean_alpha;     float mean_alpha;
Line 1086  int findPosition(DRMS_Record_t *inRec, s
Line 1093  int findPosition(DRMS_Record_t *inRec, s
  
         /* Size */         /* Size */
     // Rounded to 1.d3 precision first. Jun 16 2014 XS     // Rounded to 1.d3 precision first. Jun 16 2014 XS
       // The previous fix does not work. LONDTMAX-LONDTMIN varies from frame to frame
       // Need to find out the maximum possible difference, MAXLONDIFF (1.2e-4)
       // Now, ncol = (maxlon-minlon)/xscale, if the decimal part is outside 0.5 \pm (MAXLONDIFF/xscale)
       // proceed as it is. else, all use floor on ncol
   
           float dpix = (MAXLONDIFF / mInfo->xscale) * 1.5;                // "danger zone"
           float ncol = (maxlon - minlon) / mInfo->xscale;
           float d_ncol = fabs(ncol - floor(ncol) - 0.5);                  // distance to 0.5
           if (d_ncol < dpix) {
                   mInfo->ncol = floor(ncol);
           } else {
                   mInfo->ncol = round(ncol);
           }
  
         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);  
           printf("xcol=%f, ncol=%d, nrow=%d\n", ncol, mInfo->ncol, mInfo->nrow);
  
         return 0;         return 0;
  
Line 1961  void computeSWIndex(struct swIndex *swKe
Line 1982  void computeSWIndex(struct swIndex *swKe
         float *bpz     = (float *) (malloc(nxny * sizeof(float)));         float *bpz     = (float *) (malloc(nxny * sizeof(float)));
         float *derx    = (float *) (malloc(nxny * sizeof(float)));         float *derx    = (float *) (malloc(nxny * sizeof(float)));
         float *dery    = (float *) (malloc(nxny * sizeof(float)));         float *dery    = (float *) (malloc(nxny * sizeof(float)));
           float *derx_los = (float *) (malloc(nxny * sizeof(float)));
           float *dery_los = (float *) (malloc(nxny * sizeof(float)));
         float *derx_bt = (float *) (malloc(nxny * sizeof(float)));         float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
         float *dery_bt = (float *) (malloc(nxny * sizeof(float)));         float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
         float *derx_bh = (float *) (malloc(nxny * sizeof(float)));         float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
Line 1973  void computeSWIndex(struct swIndex *swKe
Line 1996  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)));
         float *err_term1  = (float *) (malloc(nxny * sizeof(float)));          float *err_term1   = (float *) (calloc(nxny, sizeof(float)));
         float *err_term2  = (float *) (malloc(nxny * sizeof(float)));          float *err_term2   = (float *) (calloc(nxny, sizeof(float)));
           float *err_termA   = (float *) (calloc(nxny, sizeof(float)));
           float *err_termB   = (float *) (calloc(nxny, sizeof(float)));
           float *err_termAt  = (float *) (calloc(nxny, sizeof(float)));
           float *err_termBt  = (float *) (calloc(nxny, sizeof(float)));
           float *err_termAh  = (float *) (calloc(nxny, sizeof(float)));
           float *err_termBh  = (float *) (calloc(nxny, sizeof(float)));
  
         // define some values for the R calculation         // define some values for the R calculation
         int scale = round(2.0/cdelt1);         int scale = round(2.0/cdelt1);
Line 2008  void computeSWIndex(struct swIndex *swKe
Line 2037  void computeSWIndex(struct swIndex *swKe
                 swKeys_ptr->count_mask  = DRMS_MISSING_INT;                 swKeys_ptr->count_mask  = DRMS_MISSING_INT;
         }         }
  
           if (computeAbsFlux_los(los, dims, &(swKeys_ptr->absFlux_los), &(swKeys_ptr->mean_vf_los),
                              &(swKeys_ptr->count_mask_los), bitmask, cdelt1, rsun_ref, rsun_obs))
           {
                   swKeys_ptr->absFlux_los = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
                   swKeys_ptr->mean_vf_los = DRMS_MISSING_FLOAT;
                   swKeys_ptr->count_mask_los = 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);
  
Line 2022  void computeSWIndex(struct swIndex *swKe
Line 2059  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,         if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt,
                                     dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err)))                                      dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err), err_termAt, err_termBt))
         {         {
                 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),         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_err), mask, bitmask, derx_bh, dery_bh, err_termAh, err_termBh))
         {         {
                 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),         if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err),
                                 mask, bitmask, derx_bz, dery_bz))                                  mask, bitmask, derx_bz, dery_bz, err_termA, err_termB))
         {         {
                 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;
         }         }
  
           if (computeLOSderivative(los, dims, &(swKeys_ptr->mean_derivative_los), bitmask, derx_los, dery_los))
           {
                   swKeys_ptr->mean_derivative_los = DRMS_MISSING_FLOAT; // If fail, fill in NaN
           }
   
         computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,         computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
                   derx, dery, err_term1, err_term2);                   derx, dery, err_term1, err_term2);
  
Line 2122  void computeSWIndex(struct swIndex *swKe
Line 2164  void computeSWIndex(struct swIndex *swKe
  
         }         }
  
   
         // 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 2141  void computeSWIndex(struct swIndex *swKe
Line 2182  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(derx_los); free(dery_los);
         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 numerical derivatives
           free(err_term2);
           free(err_term1);
           free(err_termB);
           free(err_termA);
           free(err_termBt);
           free(err_termAt);
           free(err_termBh);
           free(err_termAh);
   
         // free the arrays that are related to the r calculation         // free the arrays that are related to the r calculation
         free(rim);         free(rim);
         free(p1p0);         free(p1p0);
Line 2168  void computeSWIndex(struct swIndex *swKe
Line 2220  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)
 { {
     drms_setkey_float(outRec, "USFLUX",  swKeys_ptr->mean_vf);     drms_setkey_float(outRec, "USFLUX",  swKeys_ptr->mean_vf);
       drms_setkey_float(outRec, "USFLUXL", swKeys_ptr->mean_vf_los);
     drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);     drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
     drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);     drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
     drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);     drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
     drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);     drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
       drms_setkey_float(outRec, "MEANGBL", swKeys_ptr->mean_derivative_los);
     drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);     drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
     drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);     drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
     drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);     drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
Line 2184  void setSWIndex(DRMS_Record_t *outRec, s
Line 2238  void setSWIndex(DRMS_Record_t *outRec, s
     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, "CMASK",   swKeys_ptr->count_mask);
       drms_setkey_float(outRec, "CMASKL",  swKeys_ptr->count_mask_los);
     drms_setkey_float(outRec, "ERRBT",   swKeys_ptr->mean_derivative_btotal_err);     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, "ERRVF",   swKeys_ptr->mean_vf_err);
     drms_setkey_float(outRec, "ERRGAM",  swKeys_ptr->mean_gamma_err);     drms_setkey_float(outRec, "ERRGAM",  swKeys_ptr->mean_gamma_err);


Legend:
Removed from v.1.32  
changed lines
  Added in v.1.39

Karen Tian
Powered by
ViewCVS 0.9.4