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

version 1.20, 2013/11/02 19:53:05 version 1.24, 2014/02/18 23:18:07
Line 309  int computeBtotalderivative(float *bt, i
Line 309  int computeBtotalderivative(float *bt, i
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
             {             {
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                if isnan(derx_bt[j * nx + i]) continue;                 if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue;
                if isnan(dery_bt[j * nx + i]) continue;  
                if isnan(bt[j * nx + i])      continue;                if isnan(bt[j * nx + i])      continue;
                  if isnan(bt[(j+1) * nx + i])  continue;
                  if isnan(bt[(j-1) * nx + i])  continue;
                  if isnan(bt[j * nx + i-1])    continue;
                  if isnan(bt[j * nx + i+1])    continue;
                if isnan(bt_err[j * nx + i])  continue;                if isnan(bt_err[j * nx + i])  continue;
                  if isnan(derx_bt[j * nx + i]) continue;
                  if isnan(dery_bt[j * nx + i]) 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 */
                err += (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+                err += (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+
                       (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;                       (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;
Line 396  int computeBhderivative(float *bh, float
Line 401  int computeBhderivative(float *bh, float
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue;
                  if isnan(bh[j * nx + i])      continue;
                  if isnan(bh[(j+1) * nx + i])  continue;
                  if isnan(bh[(j-1) * nx + i])  continue;
                  if isnan(bh[j * nx + i-1])    continue;
                  if isnan(bh[j * nx + i+1])    continue;
                  if isnan(bh_err[j * nx + i])  continue;
                if isnan(derx_bh[j * nx + i]) continue;                if isnan(derx_bh[j * nx + i]) continue;
                if isnan(dery_bh[j * nx + i]) continue;                if isnan(dery_bh[j * nx + i]) 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 */
Line 485  int computeBzderivative(float *bz, float
Line 497  int computeBzderivative(float *bz, float
           {           {
             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 ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
                if isnan(bz[j * nx + i]) continue;                if isnan(bz[j * nx + i]) continue;
                //if isnan(bz_err[j * nx + i]) continue;                 if isnan(bz[(j+1) * nx + i])  continue;
                  if isnan(bz[(j-1) * nx + i])  continue;
                  if isnan(bz[j * nx + i-1])    continue;
                  if isnan(bz[j * nx + i+1])    continue;
                  if isnan(bz_err[j * nx + i])  continue;
                if isnan(derx_bz[j * nx + i]) continue;                if isnan(derx_bz[j * nx + i]) continue;
                if isnan(dery_bz[j * nx + i]) continue;                if isnan(dery_bz[j * nx + i]) 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 */
                err += (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i])) / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ))+                 err += (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i])) /
                       (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)])) / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;                          (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +
                         (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)])) /
                           (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;
                count_mask++;                count_mask++;
             }             }
           }           }
Line 797  int computeHelicity(float *jz_err, float
Line 815  int computeHelicity(float *jz_err, float
         *total_us_ih_err_ptr  = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // error in the quantity TOTUSJH         *total_us_ih_err_ptr  = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // error in the quantity TOTUSJH
         *total_abs_ih_err_ptr = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // error in the quantity ABSNJZH         *total_abs_ih_err_ptr = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // error in the quantity ABSNJZH
  
         printf("MEANJZH=%f\n",*mean_ih_ptr);          //printf("MEANJZH=%f\n",*mean_ih_ptr);
         printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);          //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
  
         printf("TOTUSJH=%f\n",*total_us_ih_ptr);          //printf("TOTUSJH=%f\n",*total_us_ih_ptr);
         printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);          //printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
  
         printf("ABSNJZH=%f\n",*total_abs_ih_ptr);          //printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
         printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);          //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
  
         return 0;         return 0;
 } }
Line 926  int computeFreeEnergy(float *bx_err, flo
Line 944  int computeFreeEnergy(float *bx_err, flo
  
 int computeShearAngle(float *bx_err, float *by_err, float *bz_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, int computeShearAngle(float *bx_err, float *by_err, float *bz_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
                       float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)                       float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
   
   
 { {
         int nx = dims[0];         int nx = dims[0];
         int ny = dims[1];         int ny = dims[1];
         int i = 0;         int i = 0;
         int j = 0;         int j = 0;
         int count_mask = 0;          float count_mask = 0;
           float count = 0;
         double dotproduct = 0.0;         double dotproduct = 0.0;
         double magnitude_potential = 0.0;         double magnitude_potential = 0.0;
         double magnitude_vector = 0.0;         double magnitude_vector = 0.0;
         double shear_angle = 0.0;         double shear_angle = 0.0;
         double denominator = 0.0;         double denominator = 0.0;
         double err = 0.0;  
         double sum = 0.0;  
         double count = 0.0;  
         double term1 = 0.0;         double term1 = 0.0;
         double term2 = 0.0;         double term2 = 0.0;
         double term3 = 0.0;         double term3 = 0.0;
           double sumsum = 0.0;
           double err = 0.0;
           double part1 = 0.0;
           double part2 = 0.0;
           double part3 = 0.0;
         *area_w_shear_gt_45ptr = 0.0;         *area_w_shear_gt_45ptr = 0.0;
         *meanshear_angleptr = 0.0;         *meanshear_angleptr = 0.0;
  
Line 959  int computeShearAngle(float *bx_err, flo
Line 982  int computeShearAngle(float *bx_err, flo
                  if isnan(bz[j * nx + i]) continue;                  if isnan(bz[j * nx + i]) continue;
                  if isnan(bx[j * nx + i]) continue;                  if isnan(bx[j * nx + i]) continue;
                  if isnan(by[j * nx + i]) continue;                  if isnan(by[j * nx + i]) continue;
                    if isnan(bx_err[j * nx + i]) continue;
                    if isnan(by_err[j * nx + i]) continue;
                    if isnan(bz_err[j * nx + i]) continue;
   
                  /* For mean 3D shear angle, percentage with shear greater than 45*/                  /* For mean 3D shear angle, percentage with shear greater than 45*/
                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);
                  magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));                  magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));
                  magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );                  magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );
                    //printf("dotproduct=%f\n",dotproduct);
                    //printf("magnitude_potential=%f\n",magnitude_potential);
                    //printf("magnitude_vector=%f\n",magnitude_vector);
   
                  shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);                  shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
                  sum                   += shear_angle ;                   sumsum                  += shear_angle;
                    //printf("shear_angle=%f\n",shear_angle);
                  count ++;                  count ++;
   
                  if (shear_angle > 45) count_mask ++;                  if (shear_angle > 45) count_mask ++;
  
                  /* For the error analysis*/                   // For the error analysis
                  // terms 1,2, and 3 are not squared  
                  term1 = -by[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bx[j * nx + i]*by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*(bx[j * nx + i]*bpz[j * nx + i] - bz[j * nx + i]*bpx[j * nx + i]);                   term1 = bx[j * nx + i]*by[j * nx + i]*bpy[j * nx + i] - by[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bz[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bz[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i];
                  term2 =  bx[j * nx + i]*bx[j * nx + i]*bpy[j * nx + i] - bx[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bz[j * nx + i]*(bz[j * nx + i]*bpy[j * nx + i] - by[j * nx + i]*bpz[j * nx + i]);                   term2 = bx[j * nx + i]*bx[j * nx + i]*bpy[j * nx + i] - bx[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bx[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i] - bz[j * nx + i]*by[j * nx + i]*bpz[j * nx + i];
                  term3 =  bx[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bx[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*(by[j * nx + i]*bpz[j * nx + i] - bz[j * nx + i]*bpy[j * nx + i]);                   term3 = bx[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bx[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*by[j * nx + i]*bpz[j * nx + i] - by[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i];
   
                    part1 = bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i];
                    part2 = bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] + bpz[j * nx + i]*bpz[j * nx + i];
                    part3 = bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i];
  
                  // denominator is squared                  // denominator is squared
                  denominator = (bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]) *                   denominator = part1*part1*part1*part2*(1.0-((part3*part3)/(part1*part2)));
                                (bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]) *  
                                (bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]) *                   err = (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
                                (bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] * bpz[j * nx + i]*bpz[j * nx + i]) *                         (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
                                ( 1-( ((bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i])  *                         (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) ;
                                       (bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i])) /  
                                      ((bx[j * nx + i]*bx[j * nx + i]  + by[j * nx + i]*by[j * nx + i]  + bz[j * nx + i]*bz[j * nx + i]) * (bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] + bpz[j * nx + i]*bpz[j * nx + i]))  ));  
   
                  err += ((term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator)) +  
                         ((term2*term2*by_err[j * nx + i]*by_err[j * nx + i])/(denominator)) +  
                         ((term3*term3*bz_err[j * nx + i]*bz_err[j * nx + i])/(denominator)) ;  
               }               }
           }           }
   
         /* For mean 3D shear angle, area with shear greater than 45*/         /* For mean 3D shear angle, area with shear greater than 45*/
         *meanshear_angleptr = (sum)/(count);                 /* Units are degrees */              *meanshear_angleptr = (sumsum)/(count);                 /* Units are degrees */
         *meanshear_angle_err_ptr = ((sqrt(err))/(count))*(180./PI);          *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);
  
         /* The area here is a fractional area -- the % of the total area. This has no error associated with it. */         /* The area here is a fractional area -- the % of the total area. This has no error associated with it. */
         *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);         *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);
  
         printf("MEANSHR=%f\n",*meanshear_angleptr);          //printf("MEANSHR=%f\n",*meanshear_angleptr);
         printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);          //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
           //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);
   
           return 0;
   }
   
   /*===========================================*/
   int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1,
                float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
                float *pmap, int nx1, int ny1)
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int index;
       double sum = 0.0;
       double err = 0.0;
       *Rparam = 0.0;
       struct fresize_struct fresboxcar, fresgauss;
       int scale = round(2.0/cdelt1);
       float sigma = 10.0/2.3548;
   
       init_fresize_boxcar(&fresboxcar,1,1);
   
       // set up convolution kernel
       init_fresize_gaussian(&fresgauss,sigma,20,1);
   
       // make sure convolution kernel is smaller than or equal to array size
       if ((nx || ny) < 41.) return -1;
   
       fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
       for (i = 0; i < nx1; i++)
       {
         for (j = 0; j < ny1; j++)
         {
           index = j * nx1 + i;
           if (rim[index] > 150)
             p1p0[index]=1.0;
           else
             p1p0[index]=0.0;
           if (rim[index] < -150)
             p1n0[index]=1.0;
           else
             p1n0[index]=0.0;
         }
       }
   
       fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
       fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
   
       for (i = 0; i < nx1; i++)
       {
         for (j = 0; j < ny1; j++)
         {
           index = j * nx1 + i;
           if (p1p[index] > 0 && p1n[index] > 0)
             p1[index]=1.0;
           else
             p1[index]=0.0;
         }
       }
   
       fresize(&fresgauss, p1, pmap, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
   
       for (i = 0; i < nx1; i++)
       {
         for (j = 0; j < ny1; j++)
         {
           index = j * nx1 + i;
           sum += pmap[index]*abs(rim[index]);
         }
       }
   
       if (sum < 1.0)
         *Rparam = 0.0;
       else
         *Rparam = log10(sum);
   
       free_fresize(&fresboxcar);
       free_fresize(&fresgauss);
  
         return 0;         return 0;
 } }


Legend:
Removed from v.1.20  
changed lines
  Added in v.1.24

Karen Tian
Powered by
ViewCVS 0.9.4