(file) Return to sw_functions.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

Diff for /JSOC/proj/sharp/apps/sw_functions.c between version 1.2 and 1.3

version 1.2, 2012/08/27 19:55:49 version 1.3, 2012/10/23 18:42:56
Line 18 
Line 18 
     TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter     TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
     MEANSHR Mean shear angle (measured using Btotal) in degrees     MEANSHR Mean shear angle (measured using Btotal) in degrees
  
    The indices are calculated on the pixels in which the disambig bitmap equals 5 or 7:     The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
     5: pixels for which the radial acute disambiguation solution was chosen     pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
     7: pixels for which the radial acute and NRWA disambiguation agree     coordinate bitmaps are interpolated.
   
      In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
      and the pixels that equal 33 or 44 in bitmap. Here are the definitions of the pixel values:
   
      For conf_disambig:
      50 : not all solutions agree (weak field method applied)
      60 : not all solutions agree (weak field + annealed)
      90 : all solutions agree (strong field + annealed)
       0 : not disambiguated
   
      For bitmap:
      1  : weak field outside smooth bounding curve
      2  : strong field outside smooth bounding curve
      33 : weak field inside smooth bounding curve
      34 : strong field inside smooth bounding curve
  
    Written by Monica Bobra 15 August 2012    Written by Monica Bobra 15 August 2012
    Potential Field code (appended) written by Keiji Hayashi    Potential Field code (appended) written by Keiji Hayashi
Line 52 
Line 67 
 //  7: pixels for which the radial acute and NRWA disambiguation agree //  7: pixels for which the radial acute and NRWA disambiguation agree
  
 int computeAbsFlux(float *bz, int *dims, float *absFlux, int computeAbsFlux(float *bz, int *dims, float *absFlux,
                                    float *mean_vf_ptr, int *mask,                                     float *mean_vf_ptr, int *mask,  int *bitmask,
                                    float cdelt1, double rsun_ref, double rsun_obs)                                    float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
Line 70  int computeAbsFlux(float *bz, int *dims,
Line 85  int computeAbsFlux(float *bz, int *dims,
         {         {
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                   if ( mask[j * nx + i] < 70 ) continue;                    if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                   sum += (fabs(bz[j * nx + i]));                   sum += (fabs(bz[j * nx + i]));
                   //sum += (fabs(bz[j * nx + i]))*inverseMu[j * nx + i]; // use this with mu function  
                   count_mask++;                   count_mask++;
                 }                 }
         }         }
Line 88  int computeAbsFlux(float *bz, int *dims,
Line 102  int computeAbsFlux(float *bz, int *dims,
 // Native units of Bh are Gauss // Native units of Bh are Gauss
  
 int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, int computeBh(float *bx, float *by, float *bz, float *bh, int *dims,
                           float *mean_hf_ptr, int *mask)                            float *mean_hf_ptr, int *mask, int *bitmask)
  
 { {
  
Line 120  int computeBh(float *bx, float *by, floa
Line 134  int computeBh(float *bx, float *by, floa
  
  
 int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims, int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims,
                                  float *mean_gamma_ptr, int *mask)                                   float *mean_gamma_ptr, int *mask, int *bitmask)
 { {
     int nx = dims[0], ny = dims[1];     int nx = dims[0], ny = dims[1];
     int i, j, count_mask=0;     int i, j, count_mask=0;
Line 137  int computeGamma(float *bx, float *by, f
Line 151  int computeGamma(float *bx, float *by, f
               {               {
                 if (bh[j * nx + i] > 100)                 if (bh[j * nx + i] > 100)
                   {                   {
                     if (mask[j * nx + i] < 70 ) continue;                      if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
                     count_mask++;                     count_mask++;
                   }                   }
Line 153  int computeGamma(float *bx, float *by, f
Line 167  int computeGamma(float *bx, float *by, f
 /* Example function 4: Calculate B_Total*/ /* Example function 4: Calculate B_Total*/
 // Native units of B_Total are in gauss // Native units of B_Total are in gauss
  
 int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask)  int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask, int *bitmask)
 { {
  
     int nx = dims[0], ny = dims[1];     int nx = dims[0], ny = dims[1];
Line 174  int computeB_total(float *bx, float *by,
Line 188  int computeB_total(float *bx, float *by,
 /*===========================================*/ /*===========================================*/
 /* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ /* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
  
 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, float *derx_bt, float *dery_bt)  int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt)
 { {
  
     int nx = dims[0], ny = dims[1];     int nx = dims[0], ny = dims[1];
Line 250  int computeBtotalderivative(float *bt, i
Line 264  int computeBtotalderivative(float *bt, i
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bt[j * nx + i]-dery_bt[j * nx + i]) == 0) continue;                // if ( (derx_bt[j * nx + i]-dery_bt[j * nx + i]) == 0) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */
                count_mask++;                count_mask++;
             }             }
Line 265  int computeBtotalderivative(float *bt, i
Line 279  int computeBtotalderivative(float *bt, i
 /*===========================================*/ /*===========================================*/
 /* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ /* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
  
 int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask, float *derx_bh, float *dery_bh)  int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh)
 { {
  
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 340  int computeBhderivative(float *bh, int *
Line 354  int computeBhderivative(float *bh, int *
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;                // if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */
                count_mask++;                count_mask++;
             }             }
Line 353  int computeBhderivative(float *bh, int *
Line 367  int computeBhderivative(float *bh, int *
 /*===========================================*/ /*===========================================*/
 /* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ /* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
  
 int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask, float *derx_bz, float *dery_bz)  int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz)
 { {
  
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 428  int computeBzderivative(float *bz, int *
Line 442  int computeBzderivative(float *bz, int *
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;                // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ); /* Units of Gauss */
                count_mask++;                count_mask++;
             }             }
Line 471  int computeBzderivative(float *bz, int *
Line 485  int computeBzderivative(float *bz, int *
 //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.) //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
  
 int computeJz(float *bx, float *by, int *dims, float *jz, int computeJz(float *bx, float *by, int *dims, float *jz,
                           float *mean_jz_ptr, float *us_i_ptr, int *mask,                            float *mean_jz_ptr, float *us_i_ptr, int *mask, int *bitmask,
                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
  
  
Line 552  int computeJz(float *bx, float *by, int
Line 566  int computeJz(float *bx, float *by, int
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx[j * nx + i]-dery[j * nx + i]) == 0) continue;                // if ( (derx[j * nx + i]-dery[j * nx + i]) == 0) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                curl +=     (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */                curl +=     (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
                us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */                us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */
                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */
Line 582  int computeJz(float *bx, float *by, int
Line 596  int computeJz(float *bx, float *by, int
 //                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6) //                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
 //                               = 1/Mm //                               = 1/Mm
  
 int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask,  int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask, int *bitmask,
                  float cdelt1, double rsun_ref, double rsun_obs)                  float cdelt1, double rsun_ref, double rsun_obs)
 { {
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 597  int computeAlpha(float *bz, int *dims, f
Line 611  int computeAlpha(float *bz, int *dims, f
           {           {
             for (j = 1; j < ny-1; j++)             for (j = 1; j < ny-1; j++)
               {               {
                 if (mask[j * nx + i] < 70 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(jz[j * nx + i]) continue;                 if isnan(jz[j * nx + i]) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
Line 624  int computeAlpha(float *bz, int *dims, f
Line 638  int computeAlpha(float *bz, int *dims, f
  
  
 int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr, int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr,
                                         float *total_abs_ih_ptr, int *mask, float cdelt1, double rsun_ref, double rsun_obs)                                          float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
  
Line 640  int computeHelicity(float *bz, int *dims
Line 654  int computeHelicity(float *bz, int *dims
         {         {
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                 if (mask[j * nx + i] < 70 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(jz[j * nx + i]) continue;                 if isnan(jz[j * nx + i]) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
Line 670  int computeHelicity(float *bz, int *dims
Line 684  int computeHelicity(float *bz, int *dims
 //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
  
 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr, int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,
                                                          int *mask, float cdelt1, double rsun_ref, double rsun_obs)                                                           int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 685  int computeSumAbsPerPolarity(float *bz,
Line 699  int computeSumAbsPerPolarity(float *bz,
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                 if (mask[j * nx + i] < 70 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
               }               }
Line 708  int computeSumAbsPerPolarity(float *bz,
Line 722  int computeSumAbsPerPolarity(float *bz,
 // = erg/cm(1/pix^2) // = erg/cm(1/pix^2)
  
 int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims,
                                           float *meanpotptr, float *totpotptr, int *mask,                                            float *meanpotptr, float *totpotptr, int *mask, int *bitmask,
                                           float cdelt1, double rsun_ref, double rsun_obs)                                           float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
Line 725  int computeFreeEnergy(float *bx, float *
Line 739  int computeFreeEnergy(float *bx, float *
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                  if (mask[j * nx + i] < 70 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  sum += ((    ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) -  ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i]))  )/8.*PI);                  sum += ((    ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) -  ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i]))  )/8.*PI);
                  count_mask++;                  count_mask++;
               }               }
Line 742  int computeFreeEnergy(float *bx, float *
Line 756  int computeFreeEnergy(float *bx, float *
 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
                                           float *meanshear_angleptr, float *area_w_shear_gt_45ptr,                                           float *meanshear_angleptr, float *area_w_shear_gt_45ptr,
                                           float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,                                           float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,
                                           int *mask)                                            int *mask, int *bitmask)
 { {
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
         int i, j;         int i, j;
Line 758  int computeShearAngle(float *bx, float *
Line 772  int computeShearAngle(float *bx, float *
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                  if (mask[j * nx + i] < 70 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if isnan(bpx[j * nx + i]) continue;                  if isnan(bpx[j * nx + i]) continue;
                  if isnan(bpy[j * nx + i]) continue;                  if isnan(bpy[j * nx + i]) continue;
                  if isnan(bpz[j * nx + i]) continue;                  if isnan(bpz[j * nx + i]) continue;

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