Return to sw_functions.c CVS log 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
//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("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]  )) ;
}         }
}     }
 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]  )) ;
}         }
}     }
 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]  )) ;
}         }
}     }
 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