(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.36 and 1.38

version 1.36, 2015/10/30 18:53:02 version 1.38, 2020/06/30 22:38:17
Line 1 
Line 1 
   
 /*=========================================== /*===========================================
  
  The following 14 functions calculate the following spaceweather indices:  The following 14 functions calculate the following spaceweather indices:
Line 196  int computeGamma(float *bz_err, float *b
Line 197  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 249  int computeB_total(float *bx_err, float
Line 247  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_term1, float *err_term2)  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 nx = dims[0];     int nx = dims[0];
Line 269  int computeBtotalderivative(float *bt, i
Line 267  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_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)])) ;             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)])) ;
         }         }
     }     }
  
Line 279  int computeBtotalderivative(float *bt, i
Line 277  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_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])) ;             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])) ;
         }         }
     }     }
  
Line 326  int computeBtotalderivative(float *bt, i
Line 324  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_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 += 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_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]  )) ;                     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]  )) ;
             count_mask++;             count_mask++;
         }         }
     }     }
Line 344  int computeBtotalderivative(float *bt, i
Line 342  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_term1, float *err_term2)  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 nx = dims[0];     int nx = dims[0];
Line 364  int computeBhderivative(float *bh, float
Line 362  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_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)]));             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)]));
         }         }
     }     }
  
Line 374  int computeBhderivative(float *bh, float
Line 372  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_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]));            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]));
        }        }
     }     }
  
Line 421  int computeBhderivative(float *bh, float
Line 419  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_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 += 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_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]  )) ;                     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]  )) ;
             count_mask++;             count_mask++;
         }         }
     }     }
Line 438  int computeBhderivative(float *bh, float
Line 436  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_term1, float *err_term2)  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 nx = dims[0];     int nx = dims[0];
Line 458  int computeBzderivative(float *bz, float
Line 456  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_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)]));             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)]));
         }         }
     }     }
  
Line 468  int computeBzderivative(float *bz, float
Line 466  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_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]));             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]));
         }         }
     }     }
  
Line 515  int computeBzderivative(float *bz, float
Line 513  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_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 += 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_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]  )) ;                     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]  )) ;
             count_mask++;             count_mask++;
         }         }
     }     }
Line 992  int computeShearAngle(float *bx_err, flo
Line 990  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 1241  int computeLorentz(float *bx, float *by
Line 1234  int computeLorentz(float *bx, float *by
  
 } }
  
   /*===========================================*/
   
   /* Example function 17: Compute total unsigned flux in units of G/cm^2 on the LOS field */
   
   //  To compute the unsigned flux, we simply calculate
   //  flux = surface integral [(vector LOS) dot (normal vector)],
   //       = surface integral [(magnitude LOS)*(magnitude normal)*(cos theta)].
   //  However, since the field is radial, we will assume cos theta = 1.
   //  Therefore the pixels only need to be corrected for the projection.
   
   //  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
   //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
   //  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
   //  =Gauss*cm^2
   
   int computeAbsFlux_los(float *los, int *dims, float *absFlux,
                          float *mean_vf_ptr, float *count_mask_ptr,
                          int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
   
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       *absFlux = 0.0;
       *mean_vf_ptr = 0.0;
   
   
       if (nx <= 0 || ny <= 0) return 1;
   
           for (i = 0; i < nx; i++)
           {
              for (j = 0; j < ny; j++)
              {
               if ( bitmask[j * nx + i] < 30 ) continue;
               if isnan(los[j * nx + i]) continue;
               sum += (fabs(los[j * nx + i]));
               count_mask++;
              }
           }
   
       *mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
       *count_mask_ptr  = count_mask;
   
       return 0;
   }
   
   /*===========================================*/
   /* Example function 18:  Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */
   
   int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los)
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       *mean_derivative_los_ptr = 0.0;
   
       if (nx <= 0 || ny <= 0) return 1;
   
       /* brute force method of calculating the derivative (no consideration for edges) */
       for (i = 1; i <= nx-2; i++)
       {
           for (j = 0; j <= ny-1; j++)
           {
              derx_los[j * nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])*0.5;
           }
       }
   
       /* brute force method of calculating the derivative (no consideration for edges) */
       for (i = 0; i <= nx-1; i++)
       {
           for (j = 1; j <= ny-2; j++)
           {
              dery_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
           }
       }
   
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
       ignore the edges for the error terms as those arrays have been initialized to zero.
       this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
       i=0;
       for (j = 0; j <= ny-1; j++)
       {
           derx_los[j * nx + i] = ( (-3*los[j * nx + i]) + (4*los[j * nx + (i+1)]) - (los[j * nx + (i+2)]) )*0.5;
       }
   
       i=nx-1;
       for (j = 0; j <= ny-1; j++)
       {
           derx_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[j * nx + (i-1)]) - (-los[j * nx + (i-2)]) )*0.5;
       }
   
       j=0;
       for (i = 0; i <= nx-1; i++)
       {
           dery_los[j * nx + i] = ( (-3*los[j*nx + i]) + (4*los[(j+1) * nx + i]) - (los[(j+2) * nx + i]) )*0.5;
       }
   
       j=ny-1;
       for (i = 0; i <= nx-1; i++)
       {
           dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
       }
   
   
       for (i = 0; i <= nx-1; i++)
       {
           for (j = 0; j <= ny-1; j++)
           {
               if ( bitmask[j * nx + i] < 30 ) continue;
               if ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) continue;
               if isnan(los[j * nx + i])      continue;
               if isnan(los[(j+1) * nx + i])  continue;
               if isnan(los[(j-1) * nx + i])  continue;
               if isnan(los[j * nx + i-1])    continue;
               if isnan(los[j * nx + i+1])    continue;
               if isnan(derx_los[j * nx + i]) continue;
               if isnan(dery_los[j * nx + i]) continue;
               sum += sqrt( derx_los[j * nx + i]*derx_los[j * nx + i]  + dery_los[j * nx + i]*dery_los[j * nx + i] ); /* Units of Gauss */
               count_mask++;
           }
       }
   
       *mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
       //printf("mean_derivative_los_ptr=%f\n",*mean_derivative_los_ptr);
   
           return 0;
   }
   
 /*==================KEIJI'S CODE =========================*/ /*==================KEIJI'S CODE =========================*/
  
 // #include <omp.h> // #include <omp.h>
Line 1322  void greenpot(float *bx, float *by, floa
Line 1451  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<2;j2++){for (i2=i2s;i2<2;i2++)              for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;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 1348  void greenpot(float *bx, float *by, floa
Line 1478  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 1365  char *sw_functions_version() // Returns
Line 1491  char *sw_functions_version() // Returns
     return strdup("$Id");     return strdup("$Id");
 } }
  
   
 /* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/

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