|
version 1.20, 2013/11/02 19:53:05
|
version 1.21, 2013/11/11 23:18:40
|
| 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; | return 0; |
| } | } |
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps