(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.35 and 1.36

version 1.35, 2015/03/02 21:41:31 version 1.36, 2015/10/30 18:53:02
Line 196  int computeGamma(float *bz_err, float *b
Line 196  int computeGamma(float *bz_err, float *b
     *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);     *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);
     //printf("MEANGAM=%f\n",*mean_gamma_ptr);     //printf("MEANGAM=%f\n",*mean_gamma_ptr);
     //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);     //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
       printf("sum,count_mask=%f,%d\n",sum,count_mask);
       printf("bh[200,200],bh_err[200,200]=%f,%f\n",bh[200,200],bh_err[200,200]);
   
     return 0;     return 0;
 } }
  
Line 246  int computeB_total(float *bx_err, float
Line 249  int computeB_total(float *bx_err, float
 /*===========================================*/ /*===========================================*/
 /* 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, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr, float *err_termAt, float *err_termBt)  int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr, float *err_term1, float *err_term2)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 266  int computeBtotalderivative(float *bt, i
Line 269  int computeBtotalderivative(float *bt, i
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
         {         {
            derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;            derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
            err_termAt[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)])) ;             err_term1[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)])) ;
         }         }
     }     }
  
Line 276  int computeBtotalderivative(float *bt, i
Line 279  int computeBtotalderivative(float *bt, i
         for (j = 1; j <= ny-2; j++)         for (j = 1; j <= ny-2; j++)
         {         {
            dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;            dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
            err_termBt[j * nx + i] = (((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])) ;             err_term2[j * nx + i] = (((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])) ;
         }         }
     }     }
  
Line 323  int computeBtotalderivative(float *bt, i
Line 326  int computeBtotalderivative(float *bt, i
             if isnan(derx_bt[j * nx + i]) continue;             if isnan(derx_bt[j * nx + i]) continue;
             if isnan(dery_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 += err_termBt[j * 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 += err_term2[j * 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_termAt[j * 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_term1[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;
             count_mask++;             count_mask++;
         }         }
     }     }
Line 341  int computeBtotalderivative(float *bt, i
Line 344  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, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh, float *err_termAh, float *err_termBh)  int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh, float *err_term1, float *err_term2)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 361  int computeBhderivative(float *bh, float
Line 364  int computeBhderivative(float *bh, float
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
         {         {
            derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;            derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
            err_termAh[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));             err_term1[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));
         }         }
     }     }
  
Line 371  int computeBhderivative(float *bh, float
Line 374  int computeBhderivative(float *bh, float
        for (j = 1; j <= ny-2; j++)        for (j = 1; j <= ny-2; j++)
        {        {
           dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;           dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
           err_termBh[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));            err_term2[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));
        }        }
     }     }
  
Line 418  int computeBhderivative(float *bh, float
Line 421  int computeBhderivative(float *bh, float
             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 */
             err += err_termBh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+              err += err_term2[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+
                    err_termAh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;                     err_term1[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;
             count_mask++;             count_mask++;
         }         }
     }     }
Line 435  int computeBhderivative(float *bh, float
Line 438  int computeBhderivative(float *bh, float
 /*===========================================*/ /*===========================================*/
 /* 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, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz, float *err_termA, float *err_termB)  int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz, float *err_term1, float *err_term2)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 455  int computeBzderivative(float *bz, float
Line 458  int computeBzderivative(float *bz, float
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
         {         {
            derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;            derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
            err_termA[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)]));             err_term1[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)]));
         }         }
     }     }
  
Line 465  int computeBzderivative(float *bz, float
Line 468  int computeBzderivative(float *bz, float
         for (j = 1; j <= ny-2; j++)         for (j = 1; j <= ny-2; j++)
         {         {
            dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;            dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
            err_termB[j * nx + i] = (((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]));             err_term2[j * nx + i] = (((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]));
         }         }
     }     }
  
Line 512  int computeBzderivative(float *bz, float
Line 515  int computeBzderivative(float *bz, float
             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 += err_termB[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]  )) +              err += err_term2[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]  )) +
                    err_termA[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]  )) ;                     err_term1[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]  )) ;
             count_mask++;             count_mask++;
         }         }
     }     }
Line 830  int computeHelicity(float *jz_err, float
Line 833  int computeHelicity(float *jz_err, float
 /* Example function 12:  Sum of Absolute Value per polarity  */ /* Example function 12:  Sum of Absolute Value per polarity  */
  
 //  The Sum of the Absolute Value per polarity is defined as the following: //  The Sum of the Absolute Value per polarity is defined as the following:
 //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per arcsecond.  //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per square arcsecond.
 //  The units of jz are in G/pix. In this case, we would have the following: //  The units of jz are in G/pix. In this case, we would have the following:
 //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF), //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
 //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
Line 989  int computeShearAngle(float *bx_err, flo
Line 992  int computeShearAngle(float *bx_err, flo
             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("i,j=%d,%d\n",i,j);
             //printf("dotproduct=%f\n",dotproduct);             //printf("dotproduct=%f\n",dotproduct);
             //printf("magnitude_potential=%f\n",magnitude_potential);             //printf("magnitude_potential=%f\n",magnitude_potential);
             //printf("magnitude_vector=%f\n",magnitude_vector);             //printf("magnitude_vector=%f\n",magnitude_vector);
               //printf("dotproduct/(magnitude_potential*magnitude_vector)=%f\n",dotproduct/(magnitude_potential*magnitude_vector));
               //printf("acos(dotproduct/(magnitude_potential*magnitude_vector))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector)));
               //printf("acos(dotproduct/(magnitude_potential*magnitude_vector)*(180./PI))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI));
  
             shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);             shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
             sumsum                  += shear_angle;             sumsum                  += shear_angle;
Line 1314  void greenpot(float *bx, float *by, floa
Line 1322  void greenpot(float *bx, float *by, floa
             i2e = inx + iwindow;             i2e = inx + iwindow;
             if (i2s <   0){i2s =   0;}             if (i2s <   0){i2s =   0;}
             if (i2e > nnx){i2e = nnx;}             if (i2e > nnx){i2e = nnx;}
               //            for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
             for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)              for (j2=j2s;j2<2;j2++){for (i2=i2s;i2<2;i2++)
             {             {
                 float val1 = bztmp[nnx*j2 + i2];                 float val1 = bztmp[nnx*j2 + i2];
                 float rr, r1, r2;  
                 //        r1 = (float)(i2 - inx);  
                 //        r2 = (float)(j2 - iny);  
                 //        rr = r1*r1 + r2*r2;  
                 //        if (rr < rwindow)  
                 //        {  
                 int   di, dj;                 int   di, dj;
                 di = abs(i2 - inx);                 di = abs(i2 - inx);
                 dj = abs(j2 - iny);                 dj = abs(j2 - iny);
                 sum = sum + val1 * rdist[nnx * dj + di] * dz;                 sum = sum + val1 * rdist[nnx * dj + di] * dz;
                 //        }                  printf("sum,val1,rdist[nnx * dj + di]=%f,%f,%f\n",sum,val1,rdist[nnx * dj + di]);
             } }             } }
             pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.             pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
         }         }
     } } // end of OpenMP parallelism     } } // end of OpenMP parallelism
  
       printf("bztmp[0]=%f\n",bztmp[0]);
       printf("bztmp[91746]=%f\n",bztmp[91746]);
       printf("pfpot[0]=%f\n",pfpot[0]);
       printf("pfpot[91746]=%f\n",pfpot[91746]);
   
     // #pragma omp parallel for private(inx)     // #pragma omp parallel for private(inx)
     for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)     for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
     {     {
Line 1341  void greenpot(float *bx, float *by, floa
Line 1348  void greenpot(float *bx, float *by, floa
         by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;         by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
     } } // end of OpenMP parallelism     } } // end of OpenMP parallelism
  
       printf("bx[0]=%f\n",bx[0]);
       printf("by[0]=%f\n",by[0]);
       printf("bx[91746]=%f\n",bx[91746]);
       printf("by[91746]=%f\n",by[91746]);
     free(rdist);     free(rdist);
     free(pfpot);     free(pfpot);
     free(bztmp);     free(bztmp);
Line 1354  char *sw_functions_version() // Returns
Line 1365  char *sw_functions_version() // Returns
     return strdup("$Id");     return strdup("$Id");
 } }
  
   
 /* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/

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