JSOC/proj/sharp/apps/sw_functions.c - annotate

Return to sw_functions.c CVS log

Up to [Development] / JSOC / proj / sharp / apps

1 mbobra 1.37
2 xudong 1.1 /*===========================================
3 xudong 1.27 4 The following 14 functions calculate the following spaceweather indices: 5 6 USFLUX Total unsigned flux in Maxwells 7 MEANGAM Mean inclination angle, gamma, in degrees 8 MEANGBT Mean value of the total field gradient, in Gauss/Mm 9 MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm 10 MEANGBH Mean value of the horizontal field gradient, in Gauss/Mm 11 MEANJZD Mean vertical current density, in mA/m2 12 TOTUSJZ Total unsigned vertical current, in Amperes 13 MEANALP Mean twist parameter, alpha, in 1/Mm 14 MEANJZH Mean current helicity in G2/m 15 TOTUSJH Total unsigned current helicity in G2/m 16 ABSNJZH Absolute value of the net current helicity in G2/m 17 SAVNCPP Sum of the Absolute Value of the Net Currents Per Polarity in Amperes 18 MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter 19 TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter 20 MEANSHR Mean shear angle (measured using Btotal) in degrees
21 mbobra 1.29 R_VALUE Karel Schrijver's R parameter
22 xudong 1.27 23 The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and 24 pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD 25 coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data 26 prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI 27 contain nearest-neighbor bitmaps). 28 29 In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig 30 and the pixels that equal 33 or 34 in bitmap. Here are the definitions of the pixel values: 31 32 For conf_disambig: 33 50 : not all solutions agree (weak field method applied) 34 60 : not all solutions agree (weak field + annealed) 35 90 : all solutions agree (strong field + annealed) 36 0 : not disambiguated 37 38 For bitmap: 39 1 : weak field outside smooth bounding curve 40 2 : strong field outside smooth bounding curve 41 33 : weak field inside smooth bounding curve 42 34 : strong field inside smooth bounding curve 43 xudong 1.27 44 Written by Monica Bobra 15 August 2012 45 Potential Field code (appended) written by Keiji Hayashi 46 Error analysis modification 21 October 2013 47 48 ===========================================*/
49 xudong 1.1 #include <math.h>
50 mbobra 1.9 #include <mkl.h>
51 xudong 1.1
52 mbobra 1.9 #define PI (M_PI)
53 xudong 1.1 #define MUNAUGHT (0.0000012566370614) /* magnetic constant / 54 55 /===========================================/ 56 57 / Example function 1: Compute total unsigned flux in units of G/cm^2 */ 58
59 xudong 1.27 // To compute the unsigned flux, we simply calculate
60 xudong 1.1 // flux = surface integral [(vector Bz) dot (normal vector)], 61 // = surface integral [(magnitude Bz)(magnitude normal)(cos theta)]. 62 // However, since the field is radial, we will assume cos theta = 1. 63 // Therefore the pixels only need to be corrected for the projection. 64 65 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel. 66 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 67 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
68 mbobra 1.20 // =Gauss*cm^2
69 xudong 1.1
70 xudong 1.27 int computeAbsFlux(float bz_err, float bz, int dims, float absFlux, 71 float mean_vf_ptr, float mean_vf_err_ptr, float count_mask_ptr, int mask,
72 mbobra 1.9 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
73 xudong 1.1 74 {
75 xudong 1.27
76 mbobra 1.14 int nx = dims[0]; 77 int ny = dims[1];
78 mbobra 1.15 int i = 0; 79 int j = 0; 80 int count_mask = 0; 81 double sum = 0.0; 82 double err = 0.0;
83 mbobra 1.14 absFlux = 0.0; 84 mean_vf_ptr = 0.0;
85 xudong 1.27 86
87 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
88 xudong 1.27 89 for (i = 0; i < nx; i++)
90 xudong 1.1 {
91 mbobra 1.34 for (j = 0; j < ny; j++) 92 {
93 xudong 1.27 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 94 if isnan(bz[j * nx + i]) continue; 95 sum += (fabs(bz[j * nx + i])); 96 err += bz_err[j * nx + i]bz_err[j nx + i]; 97 count_mask++;
98 mbobra 1.34 }
99 xudong 1.1 }
100 xudong 1.27 101 mean_vf_ptr = sumcdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 102 mean_vf_err_ptr = (sqrt(err))fabs(cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0); // error in the unsigned flux 103 *count_mask_ptr = count_mask; 104 return 0;
105 xudong 1.1 } 106 107 /===========================================/
108 mbobra 1.9 /* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
109 xudong 1.1 // Native units of Bh are Gauss 110
111 xudong 1.27 int computeBh(float bx_err, float by_err, float bh_err, float bx, float by, float bz, float bh, int dims,
112 mbobra 1.3 float mean_hf_ptr, int mask, int *bitmask)
113 xudong 1.1 114 {
115 xudong 1.27
116 mbobra 1.14 int nx = dims[0]; 117 int ny = dims[1];
118 mbobra 1.15 int i = 0;
119 xudong 1.27 int j = 0;
120 mbobra 1.15 int count_mask = 0;
121 xudong 1.27 double sum = 0.0;
122 mbobra 1.9 *mean_hf_ptr = 0.0;
123 xudong 1.27
124 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
125 xudong 1.27 126 for (i = 0; i < nx; i++) 127 {
128 mbobra 1.5 for (j = 0; j < ny; j++)
129 xudong 1.27 { 130 if isnan(bx[j * nx + i]) 131 { 132 bh[j * nx + i] = NAN; 133 bh_err[j * nx + i] = NAN; 134 continue; 135 } 136 if isnan(by[j * nx + i]) 137 { 138 bh[j * nx + i] = NAN; 139 bh_err[j * nx + i] = NAN; 140 continue; 141 } 142 bh[j * nx + i] = sqrt( bx[j * nx + i]bx[j nx + i] + by[j * nx + i]by[j nx + i] ); 143 sum += bh[j * nx + i]; 144 bh_err[j * nx + i]=sqrt( bx[j * nx + i]bx[j nx + i]bx_err[j nx + i]bx_err[j nx + i] + by[j * nx + i]by[j nx + i]by_err[j nx + i]by_err[j nx + i])/ bh[j * nx + i]; 145 count_mask++; 146 } 147 } 148
149 xudong 1.1 mean_hf_ptr = sum/(count_mask); // would be divided by nxny if shape of count_mask = shape of magnetogram
150 xudong 1.27
151 xudong 1.1 return 0; 152 } 153 154 /===========================================/ 155 /* Example function 3: Calculate Gamma in units of degrees */ 156 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
157 xudong 1.27 //
158 mbobra 1.20 // Error analysis calculations are done in radians (since derivatives are only true in units of radians), 159 // and multiplied by (180./PI) at the end for consistency in units.
160 xudong 1.1
161 mbobra 1.9 int computeGamma(float bz_err, float bh_err, float bx, float by, float bz, float bh, int dims, 162 float mean_gamma_ptr, float mean_gamma_err_ptr, int mask, int *bitmask)
163 xudong 1.1 {
164 mbobra 1.14 int nx = dims[0]; 165 int ny = dims[1];
166 mbobra 1.15 int i = 0; 167 int j = 0; 168 int count_mask = 0; 169 double sum = 0.0; 170 double err = 0.0; 171 *mean_gamma_ptr = 0.0;
172 xudong 1.27
173 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
174 xudong 1.27 175 for (i = 0; i < nx; i++) 176 { 177 for (j = 0; j < ny; j++) 178 { 179 if (bh[j * nx + i] > 100) 180 { 181 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 182 if isnan(bz[j * nx + i]) continue; 183 if isnan(bz_err[j * nx + i]) continue; 184 if isnan(bh_err[j * nx + i]) continue; 185 if isnan(bh[j * nx + i]) continue; 186 if (bz[j * nx + i] == 0) continue; 187 sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))(180./PI); 188 err += (1/(1+((bh[j nx + i]bh[j nx + i])/(bz[j * nx + i]bz[j nx + i]))))(1/(1+((bh[j nx + i]bh[j nx + i])/(bz[j * nx + i]bz[j nx + i])))) * 189 ( ((bh_err[j * nx + i]bh_err[j nx + i])/(bz[j * nx + i]bz[j nx + i])) + 190 ((bh[j * nx + i]bh[j nx + i]bz_err[j nx + i]bz_err[j nx + i])/(bz[j * nx + i]bz[j nx + i]bz[j nx + i]bz[j nx + i])) ); 191 count_mask++; 192 } 193 } 194 } 195 xudong 1.27 196 mean_gamma_ptr = sum/count_mask; 197 mean_gamma_err_ptr = (sqrt(err)/(count_mask))(180./PI); 198 //printf("MEANGAM=%f\n",mean_gamma_ptr); 199 //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr); 200 return 0;
201 xudong 1.1 } 202 203 /===========================================/ 204 /* Example function 4: Calculate B_Total*/ 205 // Native units of B_Total are in gauss 206
207 mbobra 1.9 int computeB_total(float bx_err, float by_err, float bz_err, float bt_err, float bx, float by, float bz, float bt, int dims, int mask, int *bitmask)
208 xudong 1.1 {
209 xudong 1.27
210 mbobra 1.14 int nx = dims[0]; 211 int ny = dims[1];
212 mbobra 1.15 int i = 0; 213 int j = 0; 214 int count_mask = 0;
215 xudong 1.1 216 if (nx <= 0 \|\| ny <= 0) return 1;
217 xudong 1.27 218 for (i = 0; i < nx; i++) 219 { 220 for (j = 0; j < ny; j++) 221 { 222 if isnan(bx[j * nx + i]) 223 { 224 bt[j * nx + i] = NAN; 225 bt_err[j * nx + i] = NAN; 226 continue; 227 } 228 if isnan(by[j * nx + i]) 229 { 230 bt[j * nx + i] = NAN; 231 bt_err[j * nx + i] = NAN; 232 continue; 233 } 234 if isnan(bz[j * nx + i]) 235 { 236 bt[j * nx + i] = NAN; 237 bt_err[j * nx + i] = NAN; 238 xudong 1.27 continue; 239 } 240 bt[j * nx + i] = 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]); 241 bt_err[j * nx + i]=sqrt(bx[j * nx + i]bx[j nx + i]bx_err[j nx + i]bx_err[j nx + i] + by[j * nx + i]by[j nx + i]by_err[j nx + i]by_err[j nx + i] + bz[j * nx + i]bz[j nx + i]bz_err[j nx + i]bz_err[j nx + i] ) / bt[j * nx + i]; 242 } 243 } 244 return 0;
245 xudong 1.1 } 246 247 /===========================================/ 248 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ 249
250 mbobra 1.37 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)
251 xudong 1.1 {
252 xudong 1.27
253 mbobra 1.14 int nx = dims[0]; 254 int ny = dims[1];
255 mbobra 1.15 int i = 0; 256 int j = 0; 257 int count_mask = 0;
258 xudong 1.27 double sum = 0.0;
259 mbobra 1.15 double err = 0.0;
260 mbobra 1.14 *mean_derivative_btotal_ptr = 0.0;
261 xudong 1.27
262 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
263 xudong 1.27 264 /* brute force method of calculating the derivative (no consideration for edges) */ 265 for (i = 1; i <= nx-2; i++) 266 {
267 mbobra 1.34 for (j = 0; j <= ny-1; j++)
268 xudong 1.27 {
269 mbobra 1.34 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
270 mbobra 1.37 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)])) ;
271 xudong 1.27 } 272 } 273 274 /* brute force method of calculating the derivative (no consideration for edges) */ 275 for (i = 0; i <= nx-1; i++) 276 {
277 mbobra 1.34 for (j = 1; j <= ny-2; j++)
278 xudong 1.27 {
279 mbobra 1.34 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
280 mbobra 1.37 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])) ;
281 xudong 1.27 } 282 } 283
284 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below. 285 ignore the edges for the error terms as those arrays have been initialized to zero. 286 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.*/
287 xudong 1.27 i=0; 288 for (j = 0; j <= ny-1; j++) 289 { 290 derx_bt[j * nx + i] = ( (-3bt[j nx + i]) + (4bt[j nx + (i+1)]) - (bt[j * nx + (i+2)]) )0.5; 291 } 292 293 i=nx-1; 294 for (j = 0; j <= ny-1; j++) 295 { 296 derx_bt[j nx + i] = ( (3bt[j nx + i]) + (-4bt[j nx + (i-1)]) - (-bt[j * nx + (i-2)]) )0.5; 297 } 298 299 j=0; 300 for (i = 0; i <= nx-1; i++) 301 { 302 dery_bt[j nx + i] = ( (-3bt[jnx + i]) + (4bt[(j+1) nx + i]) - (bt[(j+2) * nx + i]) )0.5; 303 } 304 305 j=ny-1; 306 for (i = 0; i <= nx-1; i++) 307 { 308 xudong 1.27 dery_bt[j nx + i] = ( (3bt[j nx + i]) + (-4bt[(j-1) nx + i]) - (-bt[(j-2) * nx + i]) )*0.5; 309 }
310 mbobra 1.34 311 // Calculate the sum only
312 xudong 1.27 for (i = 1; i <= nx-2; i++) 313 { 314 for (j = 1; j <= ny-2; j++) 315 { 316 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 317 if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue; 318 if isnan(bt[j * nx + i]) continue; 319 if isnan(bt[(j+1) * nx + i]) continue; 320 if isnan(bt[(j-1) * nx + i]) continue; 321 if isnan(bt[j * nx + i-1]) continue; 322 if isnan(bt[j * nx + i+1]) continue; 323 if isnan(bt_err[j * nx + i]) continue; 324 if isnan(derx_bt[j * nx + i]) continue; 325 if isnan(dery_bt[j * nx + i]) continue; 326 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 */
327 mbobra 1.37 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] ))+ 328 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] )) ;
329 xudong 1.27 count_mask++; 330 } 331 } 332 333 mean_derivative_btotal_ptr = (sum)/(count_mask); 334 mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); 335 //printf("MEANGBT=%f\n",mean_derivative_btotal_ptr); 336 //printf("MEANGBT_err=%f\n",mean_derivative_btotal_err_ptr); 337 338 return 0;
339 xudong 1.1 } 340 341 342 /===========================================/ 343 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ 344
345 mbobra 1.37 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)
346 xudong 1.1 {
347 xudong 1.27 348 int nx = dims[0]; 349 int ny = dims[1]; 350 int i = 0; 351 int j = 0; 352 int count_mask = 0; 353 double sum= 0.0; 354 double err =0.0; 355 mean_derivative_bh_ptr = 0.0; 356 357 if (nx <= 0 \|\| ny <= 0) return 1; 358 359 / brute force method of calculating the derivative (no consideration for edges) */ 360 for (i = 1; i <= nx-2; i++) 361 {
362 mbobra 1.34 for (j = 0; j <= ny-1; j++)
363 xudong 1.27 {
364 mbobra 1.34 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
365 mbobra 1.37 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)]));
366 xudong 1.27 } 367 } 368 369 /* brute force method of calculating the derivative (no consideration for edges) */ 370 for (i = 0; i <= nx-1; i++) 371 {
372 mbobra 1.34 for (j = 1; j <= ny-2; j++) 373 { 374 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
375 mbobra 1.37 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]));
376 mbobra 1.34 }
377 xudong 1.27 } 378
379 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below. 380 ignore the edges for the error terms as those arrays have been initialized to zero. 381 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.*/
382 xudong 1.27 i=0; 383 for (j = 0; j <= ny-1; j++) 384 { 385 derx_bh[j * nx + i] = ( (-3bh[j nx + i]) + (4bh[j nx + (i+1)]) - (bh[j * nx + (i+2)]) )0.5; 386 } 387 388 i=nx-1; 389 for (j = 0; j <= ny-1; j++) 390 { 391 derx_bh[j nx + i] = ( (3bh[j nx + i]) + (-4bh[j nx + (i-1)]) - (-bh[j * nx + (i-2)]) )0.5; 392 } 393 394 j=0; 395 for (i = 0; i <= nx-1; i++) 396 { 397 dery_bh[j nx + i] = ( (-3bh[jnx + i]) + (4bh[(j+1) nx + i]) - (bh[(j+2) * nx + i]) )0.5; 398 } 399 400 j=ny-1; 401 for (i = 0; i <= nx-1; i++) 402 { 403 xudong 1.27 dery_bh[j nx + i] = ( (3bh[j nx + i]) + (-4bh[(j-1) nx + i]) - (-bh[(j-2) * nx + i]) )0.5; 404 } 405 406 407 for (i = 0; i <= nx-1; i++) 408 { 409 for (j = 0; j <= ny-1; j++) 410 { 411 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 412 if ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue; 413 if isnan(bh[j * nx + i]) continue; 414 if isnan(bh[(j+1) * nx + i]) continue; 415 if isnan(bh[(j-1) * nx + i]) continue; 416 if isnan(bh[j * nx + i-1]) continue; 417 if isnan(bh[j * nx + i+1]) continue; 418 if isnan(bh_err[j * nx + i]) continue; 419 if isnan(derx_bh[j * nx + i]) continue; 420 if isnan(dery_bh[j * nx + i]) continue; 421 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 */
422 mbobra 1.37 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] ))+ 423 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] )) ;
424 xudong 1.27 count_mask++; 425 } 426 } 427 428 mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)(ny-2)) if shape of count_mask = shape of magnetogram 429 mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask) 430 //printf("MEANGBH=%f\n",mean_derivative_bh_ptr); 431 //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr); 432 433 return 0;
434 xudong 1.1 } 435 436 /===========================================/ 437 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ 438
439 mbobra 1.37 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)
440 xudong 1.1 {
441 xudong 1.27 442 int nx = dims[0]; 443 int ny = dims[1]; 444 int i = 0; 445 int j = 0; 446 int count_mask = 0;
447 mbobra 1.34 double sum = 0.0;
448 xudong 1.27 double err = 0.0;
449 mbobra 1.34 *mean_derivative_bz_ptr = 0.0;
450 xudong 1.27
451 mbobra 1.34 if (nx <= 0 \|\| ny <= 0) return 1;
452 xudong 1.27 453 /* brute force method of calculating the derivative (no consideration for edges) */ 454 for (i = 1; i <= nx-2; i++) 455 {
456 mbobra 1.34 for (j = 0; j <= ny-1; j++)
457 xudong 1.27 {
458 mbobra 1.34 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
459 mbobra 1.37 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)]));
460 xudong 1.27 } 461 } 462 463 /* brute force method of calculating the derivative (no consideration for edges) */ 464 for (i = 0; i <= nx-1; i++) 465 {
466 mbobra 1.34 for (j = 1; j <= ny-2; j++)
467 xudong 1.27 {
468 mbobra 1.34 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
469 mbobra 1.37 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]));
470 xudong 1.27 } 471 } 472
473 mbobra 1.34 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
474 mbobra 1.35 ignore the edges for the error terms as those arrays have been initialized to zero.
475 mbobra 1.34 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.*/
476 xudong 1.27 i=0; 477 for (j = 0; j <= ny-1; j++) 478 { 479 derx_bz[j * nx + i] = ( (-3bz[j nx + i]) + (4bz[j nx + (i+1)]) - (bz[j * nx + (i+2)]) )0.5; 480 } 481 482 i=nx-1; 483 for (j = 0; j <= ny-1; j++) 484 { 485 derx_bz[j nx + i] = ( (3bz[j nx + i]) + (-4bz[j nx + (i-1)]) - (-bz[j * nx + (i-2)]) )0.5; 486 } 487 488 j=0; 489 for (i = 0; i <= nx-1; i++) 490 { 491 dery_bz[j nx + i] = ( (-3bz[jnx + i]) + (4bz[(j+1) nx + i]) - (bz[(j+2) * nx + i]) )0.5; 492 } 493 494 j=ny-1; 495 for (i = 0; i <= nx-1; i++) 496 { 497 xudong 1.27 dery_bz[j nx + i] = ( (3bz[j nx + i]) + (-4bz[(j-1) nx + i]) - (-bz[(j-2) * nx + i]) )0.5; 498 } 499 500 501 for (i = 0; i <= nx-1; i++) 502 { 503 for (j = 0; j <= ny-1; j++) 504 { 505 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 506 if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue; 507 if isnan(bz[j * nx + i]) continue; 508 if isnan(bz[(j+1) * nx + i]) continue; 509 if isnan(bz[(j-1) * nx + i]) continue; 510 if isnan(bz[j * nx + i-1]) continue; 511 if isnan(bz[j * nx + i+1]) continue; 512 if isnan(bz_err[j * nx + i]) continue; 513 if isnan(derx_bz[j * nx + i]) continue; 514 if isnan(dery_bz[j * nx + i]) continue;
515 mbobra 1.34 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 */
516 mbobra 1.37 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] )) + 517 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] )) ;
518 xudong 1.27 count_mask++; 519 } 520 } 521
522 mbobra 1.34 mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)(ny-2)) if shape of count_mask = shape of magnetogram
523 xudong 1.27 mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask) 524 //printf("MEANGBZ=%f\n",mean_derivative_bz_ptr); 525 //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr); 526
527 xudong 1.1 return 0; 528 } 529 530 /===========================================/ 531 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */ 532 533 // In discretized space like data pixels, 534 // the current (or curl of B) is calculated as 535 // the integration of the field Bx and By along 536 // the circumference of the data pixel divided by the area of the pixel. 537 // One form of differencing (a word for the differential operator
538 xudong 1.27 // in the discretized space) of the curl is expressed as
539 xudong 1.1 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2 540 // +dy * (By(i+1,j)+By(i,j)) / 2 541 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
542 xudong 1.27 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy) 543 //
544 xudong 1.1 // 545 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003), 546 // one must perform the following unit conversions:
547 mbobra 1.8 // (Gauss)(1/arcsec)(arcsec/meter)(Newton/GaussAmperemeter)(Ampere^2/Newton)(milliAmpere/Ampere), or 548 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4PI10^-7)( 10^3 milliAmpere/Ampere), or
549 xudong 1.27 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
550 xudong 1.1 // where a Tesla is represented as a Newton/Ampere*meter.
551 mbobra 1.4 //
552 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 553 // In that case, we would have the following: 554 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4PI10^7)(10^3), or 555 // jz * (35.0) 556 // 557 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
558 mbobra 1.8 // (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS) 559 // = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
560 xudong 1.1 561
562 xudong 1.27 // Comment out random number generator, which can only run on solar3
563 mbobra 1.33 // int computeJz(float bx_err, float by_err, float bx, float by, int dims, float jz, float jz_err, float jz_err_squared,
564 xudong 1.27 // int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery, float *noisebx,
565 mbobra 1.9 // float noiseby, float noisebz) 566 567 int computeJz(float bx_err, float by_err, float bx, float by, int dims, float jz, float jz_err, float jz_err_squared,
568 mbobra 1.33 int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery, float err_term1, float err_term2)
569 mbobra 1.9
570 xudong 1.1
571 xudong 1.27 { 572 int nx = dims[0]; 573 int ny = dims[1]; 574 int i = 0; 575 int j = 0; 576 int count_mask = 0; 577 578 if (nx <= 0 \|\| ny <= 0) return 1; 579 580 /* Calculate the derivative/ 581 / brute force method of calculating the derivative (no consideration for edges) */ 582 583 for (i = 1; i <= nx-2; i++) 584 {
585 mbobra 1.34 for (j = 0; j <= ny-1; j++)
586 xudong 1.27 {
587 mbobra 1.34 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])0.5; 588 err_term1[j nx + i] = (by_err[j * nx + i+1])(by_err[j nx + i+1]) + (by_err[j * nx + i-1])(by_err[j nx + i-1]);
589 xudong 1.27 } 590 } 591 592 for (i = 0; i <= nx-1; i++) 593 {
594 mbobra 1.34 for (j = 1; j <= ny-2; j++)
595 xudong 1.27 {
596 mbobra 1.34 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])0.5; 597 err_term2[j nx + i] = (bx_err[(j+1) * nx + i])(bx_err[(j+1) nx + i]) + (bx_err[(j-1) * nx + i])(bx_err[(j-1) nx + i]);
598 xudong 1.27 } 599 }
600 mbobra 1.33
601 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below. 602 ignore the edges for the error terms as those arrays have been initialized to zero. 603 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.*/ 604
605 xudong 1.27 i=0; 606 for (j = 0; j <= ny-1; j++) 607 {
608 mbobra 1.33 derx[j * nx + i] = ( (-3by[j nx + i]) + (4by[j nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
609 xudong 1.27 } 610 611 i=nx-1; 612 for (j = 0; j <= ny-1; j++) 613 {
614 mbobra 1.33 derx[j * nx + i] = ( (3by[j nx + i]) + (-4by[j nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
615 xudong 1.27 } 616 617 j=0; 618 for (i = 0; i <= nx-1; i++) 619 {
620 mbobra 1.33 dery[j * nx + i] = ( (-3bx[jnx + i]) + (4bx[(j+1) nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
621 xudong 1.27 } 622 623 j=ny-1; 624 for (i = 0; i <= nx-1; i++) 625 {
626 mbobra 1.33 dery[j * nx + i] = ( (3bx[j nx + i]) + (-4bx[(j-1) nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
627 xudong 1.27 }
628 mbobra 1.33
629 xudong 1.27
630 mbobra 1.33 for (i = 0; i <= nx-1; i++)
631 xudong 1.27 {
632 mbobra 1.33 for (j = 0; j <= ny-1; j++)
633 xudong 1.27 {
634 mbobra 1.33 // calculate jz at all points
635 xudong 1.27 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix
636 mbobra 1.33 jz_err[j * nx + i] = 0.5sqrt( err_term1[j nx + i] + err_term2[j * nx + i] ) ;
637 xudong 1.27 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]jz_err[j nx + i]);
638 mbobra 1.33 count_mask++;
639 xudong 1.27 } 640 } 641 return 0; 642 }
643 mbobra 1.33
644 mbobra 1.5 /===========================================/ 645
646 mbobra 1.11 /* Example function 9: Compute quantities on Jz array */
647 xudong 1.27 // Compute mean and total current on Jz array.
648 mbobra 1.6
649 mbobra 1.9 int computeJzsmooth(float bx, float by, int dims, float jz, float jz_smooth, float jz_err, float jz_rms_err, float jz_err_squared_smooth,
650 xudong 1.27 float mean_jz_ptr, float mean_jz_err_ptr, float us_i_ptr, float us_i_err_ptr, int mask, int bitmask,
651 mbobra 1.9 float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery)
652 mbobra 1.5 653 {
654 xudong 1.27 655 int nx = dims[0]; 656 int ny = dims[1]; 657 int i = 0; 658 int j = 0; 659 int count_mask = 0;
660 mbobra 1.15 double curl = 0.0;
661 xudong 1.27 double us_i = 0.0; 662 double err = 0.0; 663
664 mbobra 1.5 if (nx <= 0 \|\| ny <= 0) return 1;
665 xudong 1.27 666 /* At this point, use the smoothed Jz array with a Gaussian (FWHM of 4 pix and truncation width of 12 pixels) but keep the original array dimensions/ 667 for (i = 0; i <= nx-1; i++) 668 { 669 for (j = 0; j <= ny-1; j++) 670 { 671 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 672 if isnan(derx[j * nx + i]) continue; 673 if isnan(dery[j * nx + i]) continue; 674 if isnan(jz[j * nx + i]) continue; 675 curl += (jz[j * nx + i])(1/cdelt1)(rsun_obs/rsun_ref)(0.00010)(1/MUNAUGHT)(1000.); / curl is in units of mA / m^2 / 676 us_i += fabs(jz[j nx + i])(cdelt1/1)(rsun_ref/rsun_obs)(0.00010)(1/MUNAUGHT); /* us_i is in units of A / 677 err += (jz_err[j nx + i]jz_err[j nx + i]); 678 count_mask++; 679 } 680 } 681 682 /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above / 683 mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD / 684 mean_jz_err_ptr = (sqrt(err)/count_mask)((1/cdelt1)(rsun_obs/rsun_ref)(0.00010)(1/MUNAUGHT)(1000.)); // error in the quantity MEANJZD 685 686 xudong 1.27 us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ / 687 us_i_err_ptr = (sqrt(err))((cdelt1/1)(rsun_ref/rsun_obs)(0.00010)(1/MUNAUGHT)); // error in the quantity TOTUSJZ 688 689 //printf("MEANJZD=%f\n",mean_jz_ptr); 690 //printf("MEANJZD_err=%f\n",mean_jz_err_ptr); 691 692 //printf("TOTUSJZ=%g\n",us_i_ptr); 693 //printf("TOTUSJZ_err=%g\n",us_i_err_ptr); 694
695 xudong 1.1 return 0;
696 xudong 1.27
697 xudong 1.1 } 698
699 mbobra 1.5 /===========================================/ 700 701 /* Example function 10: Twist Parameter, alpha */
702 xudong 1.1
703 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
704 mbobra 1.19 // for alpha is weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):
705 xudong 1.27 706 // numerator = sum of all JzBz 707 // denominator = sum of BzBz 708 // alpha = numerator/denominator
709 mbobra 1.6
710 mbobra 1.5 // The units of alpha are in 1/Mm
711 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss. 712 //
713 xudong 1.27 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
714 xudong 1.1 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6) 715 // = 1/Mm 716
717 mbobra 1.9 int computeAlpha(float jz_err, float bz_err, float bz, int dims, float jz, float jz_smooth, float mean_alpha_ptr, float mean_alpha_err_ptr, int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs)
718 mbobra 1.5
719 xudong 1.1 {
720 xudong 1.27 int nx = dims[0]; 721 int ny = dims[1]; 722 int i = 0; 723 int j = 0;
724 mbobra 1.19 double alpha_total = 0.0;
725 xudong 1.27 double C = ((1/cdelt1)(rsun_obs/rsun_ref)(1000000.)); 726 double total = 0.0; 727 double A = 0.0; 728 double B = 0.0; 729
730 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
731 xudong 1.27 for (i = 1; i < nx-1; i++) 732 { 733 for (j = 1; j < ny-1; j++) 734 { 735 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 736 if isnan(jz[j * nx + i]) continue; 737 if isnan(bz[j * nx + i]) continue; 738 if (jz[j * nx + i] == 0.0) continue; 739 if (bz[j * nx + i] == 0.0) continue; 740 A += jz[jnx+i]bz[jnx+i]; 741 B += bz[jnx+i]bz[jnx+i]; 742 } 743 } 744 745 for (i = 1; i < nx-1; i++) 746 { 747 for (j = 1; j < ny-1; j++) 748 { 749 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 750 if isnan(jz[j * nx + i]) continue; 751 if isnan(bz[j * nx + i]) continue; 752 xudong 1.27 if (jz[j * nx + i] == 0.0) continue; 753 if (bz[j * nx + i] == 0.0) continue; 754 total += bz[jnx+i]bz[jnx+i]jz_err[jnx+i]jz_err[jnx+i] + (jz[jnx+i]-2bz[jnx+i]A/B)(jz[jnx+i]-2bz[jnx+i]A/B)bz_err[jnx+i]bz_err[jnx+i]; 755 } 756 } 757 758 /* Determine the absolute value of alpha. The units for alpha are 1/Mm / 759 alpha_total = ((A/B)C); 760 mean_alpha_ptr = alpha_total; 761 mean_alpha_err_ptr = (C/B)*(sqrt(total)); 762
763 xudong 1.1 return 0; 764 } 765 766 /===========================================/
767 mbobra 1.9 /* Example function 11: Helicity (mean current helicty, total unsigned current helicity, absolute value of net current helicity) */
768 xudong 1.1 769 // The current helicity is defined as Bz*Jz and the units are G^2 / m 770 // The units of Jz are in G/pix; the units of Bz are in G.
771 xudong 1.27 // Therefore, the units of BzJz = (Gauss)(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter) 772 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
773 mbobra 1.9 // = G^2 / m.
774 xudong 1.1
775 xudong 1.27 int computeHelicity(float jz_err, float jz_rms_err, float bz_err, float bz, int dims, float jz, float mean_ih_ptr, 776 float mean_ih_err_ptr, float total_us_ih_ptr, float total_abs_ih_ptr,
777 mbobra 1.9 float total_us_ih_err_ptr, float total_abs_ih_err_ptr, int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs)
778 xudong 1.1 779 {
780 xudong 1.27 781 int nx = dims[0]; 782 int ny = dims[1]; 783 int i = 0; 784 int j = 0; 785 int count_mask = 0;
786 mbobra 1.20 double sum = 0.0; 787 double sum2 = 0.0; 788 double err = 0.0;
789 xudong 1.1 790 if (nx <= 0 \|\| ny <= 0) return 1;
791 xudong 1.27 792 for (i = 0; i < nx; i++)
793 xudong 1.1 {
794 xudong 1.27 for (j = 0; j < ny; j++)
795 xudong 1.1 {
796 xudong 1.27 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 797 if isnan(jz[j * nx + i]) continue; 798 if isnan(bz[j * nx + i]) continue; 799 if isnan(jz_err[j * nx + i]) continue; 800 if isnan(bz_err[j * nx + i]) continue; 801 if (bz[j * nx + i] == 0.0) continue; 802 if (jz[j * nx + i] == 0.0) continue; 803 sum += (jz[j * nx + i]bz[j nx + i])(1/cdelt1)(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH 804 sum2 += fabs(jz[j * nx + i]bz[j nx + i])(1/cdelt1)(rsun_obs/rsun_ref); // contributes to TOTUSJH 805 err += (jz_err[j * nx + i]jz_err[j nx + i]bz[j nx + i]bz[j nx + i]) + (bz_err[j * nx + i]bz_err[j nx + i]jz[j nx + i]jz[j nx + i]); 806 count_mask++; 807 } 808 } 809 810 mean_ih_ptr = sum/count_mask ; / Units are G^2 / m ; keyword is MEANJZH */
811 mbobra 1.9 total_us_ih_ptr = sum2 ; / Units are G^2 / m ; keyword is TOTUSJH / 812 total_abs_ih_ptr = fabs(sum) ; /* Units are G^2 / m ; keyword is ABSNJZH */
813 xudong 1.27 814 mean_ih_err_ptr = (sqrt(err)/count_mask)(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity MEANJZH 815 total_us_ih_err_ptr = (sqrt(err))(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity TOTUSJH 816 total_abs_ih_err_ptr = (sqrt(err))(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity ABSNJZH 817 818 //printf("MEANJZH=%f\n",mean_ih_ptr); 819 //printf("MEANJZH_err=%f\n",mean_ih_err_ptr); 820 821 //printf("TOTUSJH=%f\n",total_us_ih_ptr); 822 //printf("TOTUSJH_err=%f\n",total_us_ih_err_ptr); 823 824 //printf("ABSNJZH=%f\n",total_abs_ih_ptr); 825 //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr); 826
827 xudong 1.1 return 0; 828 } 829 830 /===========================================/
831 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
832 xudong 1.1 833 // The Sum of the Absolute Value per polarity is defined as the following:
834 mbobra 1.36 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per square arcsecond.
835 xudong 1.1 // The units of jz are in G/pix. In this case, we would have the following: 836 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF), 837 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
838 mbobra 1.9 // 839 // The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
840 xudong 1.1
841 xudong 1.27 int computeSumAbsPerPolarity(float jz_err, float bz_err, float bz, float jz, int dims, float totaljzptr, float *totaljz_err_ptr,
842 mbobra 1.3 int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs)
843 xudong 1.1
844 xudong 1.27 { 845 int nx = dims[0]; 846 int ny = dims[1]; 847 int i=0; 848 int j=0; 849 int count_mask=0;
850 mbobra 1.32 double sum1=0.0;
851 xudong 1.27 double sum2=0.0; 852 double err=0.0;
853 mbobra 1.33 *totaljzptr=0.0;
854 xudong 1.27
855 mbobra 1.33 if (nx <= 0 \|\| ny <= 0) return 1;
856 xudong 1.27
857 mbobra 1.33 for (i = 0; i < nx; i++)
858 xudong 1.27 {
859 mbobra 1.33 for (j = 0; j < ny; j++)
860 xudong 1.27 { 861 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 862 if isnan(bz[j * nx + i]) continue; 863 if isnan(jz[j * nx + i]) continue; 864 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])(1/cdelt1)(0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs); 865 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])(1/cdelt1)(0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs); 866 err += (jz_err[j * nx + i]jz_err[j nx + i]); 867 count_mask++; 868 } 869 }
870 xudong 1.1
871 mbobra 1.32 totaljzptr = fabs(sum1) + fabs(sum2); / Units are Amperes per arcsecond */
872 xudong 1.27 totaljz_err_ptr = sqrt(err)(1/cdelt1)fabs((0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs)); 873 //printf("SAVNCPP=%g\n",totaljzptr); 874 //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr); 875
876 mbobra 1.33 return 0;
877 xudong 1.1 } 878 879 /===========================================/
880 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
881 xudong 1.1 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
882 xudong 1.27 // automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8PI can come out of the integral; thus, 883 // the integral is over B^2 dV and the 8PI is divided at the end.
884 xudong 1.1 // 885 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert 886 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
887 mbobra 1.9 // erg/cm^3(CDELT1^2)(RSUN_REF/RSUN_OBS ^2)(100.^2) 888 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 889 // = erg/cm^3*(1.30501e15)
890 xudong 1.1 // = erg/cm(1/pix^2) 891
892 xudong 1.27 int computeFreeEnergy(float bx_err, float by_err, float bx, float by, float bpx, float bpy, int dims, 893 float meanpotptr, float meanpot_err_ptr, float totpotptr, float totpot_err_ptr, int mask, int *bitmask,
894 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs) 895 896 {
897 xudong 1.27 int nx = dims[0]; 898 int ny = dims[1]; 899 int i = 0; 900 int j = 0; 901 int count_mask = 0;
902 mbobra 1.31 double sum = 0.0;
903 xudong 1.27 double sum1 = 0.0; 904 double err = 0.0; 905 *totpotptr = 0.0;
906 mbobra 1.31 *meanpotptr = 0.0;
907 xudong 1.27
908 mbobra 1.31 if (nx <= 0 \|\| ny <= 0) return 1;
909 xudong 1.27
910 mbobra 1.31 for (i = 0; i < nx; i++)
911 xudong 1.27 {
912 mbobra 1.31 for (j = 0; j < ny; j++)
913 xudong 1.27 { 914 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 915 if isnan(bx[j * nx + i]) continue; 916 if isnan(by[j * nx + i]) continue; 917 sum += ( ((bx[j * nx + i] - bpx[j * nx + i])(bx[j nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])(by[j nx + i] - bpy[j * nx + i])) )(cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0); 918 sum1 += ( ((bx[j * nx + i] - bpx[j * nx + i])(bx[j nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])(by[j nx + i] - bpy[j * nx + i])) ); 919 err += 4.0(bx[j nx + i] - bpx[j * nx + i])(bx[j nx + i] - bpx[j * nx + i])(bx_err[j nx + i]bx_err[j nx + i]) + 920 4.0(by[j nx + i] - bpy[j * nx + i])(by[j nx + i] - bpy[j * nx + i])(by_err[j nx + i]by_err[j nx + i]); 921 count_mask++; 922 } 923 } 924 925 /* Units of meanpotptr are ergs per centimeter */
926 mbobra 1.20 meanpotptr = (sum1) / (count_mask8.PI) ; / Units are ergs per cubic centimeter */
927 xudong 1.27 meanpot_err_ptr = (sqrt(err)) / (count_mask8.PI); // error in the quantity (sum)/(count_mask) 928 929 / Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter / 930 totpotptr = (sum)/(8.PI); 931 totpot_err_ptr = (sqrt(err))fabs(cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0(1/(8.PI))); 932 933 //printf("MEANPOT=%g\n",meanpotptr); 934 //printf("MEANPOT_err=%g\n",meanpot_err_ptr); 935 936 //printf("TOTPOT=%g\n",totpotptr); 937 //printf("TOTPOT_err=%g\n",totpot_err_ptr); 938
939 mbobra 1.31 return 0;
940 xudong 1.1 } 941 942 /===========================================/
943 mbobra 1.5 /* Example function 14: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
944 xudong 1.1
945 mbobra 1.20 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,
946 mbobra 1.9 float meanshear_angleptr, float meanshear_angle_err_ptr, float area_w_shear_gt_45ptr, int mask, int *bitmask)
947 mbobra 1.21 948
949 xudong 1.27 { 950 int nx = dims[0]; 951 int ny = dims[1]; 952 int i = 0; 953 int j = 0; 954 float count_mask = 0; 955 float count = 0; 956 double dotproduct = 0.0; 957 double magnitude_potential = 0.0; 958 double magnitude_vector = 0.0; 959 double shear_angle = 0.0; 960 double denominator = 0.0; 961 double term1 = 0.0; 962 double term2 = 0.0; 963 double term3 = 0.0; 964 double sumsum = 0.0; 965 double err = 0.0; 966 double part1 = 0.0; 967 double part2 = 0.0; 968 double part3 = 0.0; 969 *area_w_shear_gt_45ptr = 0.0;
970 mbobra 1.31 *meanshear_angleptr = 0.0;
971 xudong 1.1
972 mbobra 1.31 if (nx <= 0 \|\| ny <= 0) return 1;
973 xudong 1.27
974 mbobra 1.31 for (i = 0; i < nx; i++)
975 xudong 1.27 {
976 mbobra 1.31 for (j = 0; j < ny; j++)
977 xudong 1.27 { 978 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 979 if isnan(bpx[j * nx + i]) continue; 980 if isnan(bpy[j * nx + i]) continue; 981 if isnan(bpz[j * nx + i]) continue; 982 if isnan(bz[j * nx + i]) continue; 983 if isnan(bx[j * nx + i]) continue; 984 if isnan(by[j * nx + i]) continue; 985 if isnan(bx_err[j * nx + i]) continue; 986 if isnan(by_err[j * nx + i]) continue; 987 if isnan(bz_err[j * nx + i]) continue; 988 989 /* For mean 3D shear angle, percentage with shear greater than 45/ 990 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]); 991 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])); 992 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]) ); 993 //printf("dotproduct=%f\n",dotproduct); 994 //printf("magnitude_potential=%f\n",magnitude_potential); 995 //printf("magnitude_vector=%f\n",magnitude_vector); 996 997 shear_angle = acos(dotproduct/(magnitude_potentialmagnitude_vector))(180./PI); 998 xudong 1.27 sumsum += shear_angle; 999 //printf("shear_angle=%f\n",shear_angle); 1000 count ++; 1001 1002 if (shear_angle > 45) count_mask ++; 1003 1004 // For the error analysis 1005 1006 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]; 1007 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]; 1008 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]; 1009 1010 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]; 1011 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]; 1012 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]; 1013 1014 // denominator is squared 1015 denominator = part1part1part1part2(1.0-((part3part3)/(part1part2))); 1016 1017 err = (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) + 1018 (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) + 1019 xudong 1.27 (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) ; 1020 1021 } 1022 } 1023 /* For mean 3D shear angle, area with shear greater than 45/ 1024 meanshear_angleptr = (sumsum)/(count); /* Units are degrees */
1025 mbobra 1.39 1026 // For the error in the mean 3D shear angle: 1027 // If count_mask is 0, then we run into a divide by zero error. In this case, set meanshear_angle_err_ptr to NAN 1028 // If count_mask is greater than zero, then compute the error. 1029 if (count_mask == 0) 1030 meanshear_angle_err_ptr = NAN; 1031 else 1032 meanshear_angle_err_ptr = (sqrt(err)/count_mask)(180./PI);
1033 xudong 1.27 1034 /* The area here is a fractional area -- the % of the total area. This has no error associated with it. / 1035 area_w_shear_gt_45ptr = (count_mask/(count))(100.0); 1036 1037 //printf("MEANSHR=%f\n",meanshear_angleptr);
1038 mbobra 1.39 //printf("ERRMSHA=%f\n",*meanshear_angle_err_ptr);
1039 xudong 1.27 //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);
1040 mbobra 1.39 return 0;
1041 xudong 1.27 }
1042 xudong 1.1
1043 xudong 1.27 /===========================================/
1044 mbobra 1.29 /* Example function 15: R parameter as defined in Schrijver, 2007 / 1045 // 1046 // Note that there is a restriction on the function fsample() 1047 // If the following occurs: 1048 // nx_out > floor((ny_in-1)/scale + 1) 1049 // ny_out > floor((ny_in-1)/scale + 1), 1050 // where n_out are the dimensions of the output array and n*_in 1051 // are the dimensions of the input array, fsample() will usually result 1052 // in a segfault (though not always, depending on how the segfault was accessed.) 1053
1054 xudong 1.27 int computeR(float bz_err, float los, int dims, float Rparam, float cdelt1, 1055 float rim, float p1p0, float p1n0, float p1p, float p1n, float p1,
1056 mbobra 1.30 float pmap, int nx1, int ny1, 1057 int scale, float p1pad, int nxp, int nyp, float *pmapn) 1058
1059 mbobra 1.29 {
1060 xudong 1.27 int nx = dims[0]; 1061 int ny = dims[1]; 1062 int i = 0; 1063 int j = 0;
1064 mbobra 1.30 int index, index1;
1065 xudong 1.27 double sum = 0.0; 1066 double err = 0.0; 1067 *Rparam = 0.0; 1068 struct fresize_struct fresboxcar, fresgauss;
1069 mbobra 1.30 struct fint_struct fints;
1070 xudong 1.27 float sigma = 10.0/2.3548; 1071
1072 mbobra 1.30 // set up convolution kernels
1073 xudong 1.27 init_fresize_boxcar(&fresboxcar,1,1); 1074 init_fresize_gaussian(&fresgauss,sigma,20,1);
1075 mbobra 1.29
1076 mbobra 1.30 // =============== [STEP 1] =============== 1077 // bin the line-of-sight magnetogram down by a factor of scale
1078 xudong 1.27 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
1079 mbobra 1.29
1080 mbobra 1.30 // =============== [STEP 2] =============== 1081 // identify positive and negative pixels greater than +/- 150 gauss 1082 // and label those pixels with a 1.0 in arrays p1p0 and p1n0
1083 xudong 1.27 for (i = 0; i < nx1; i++) 1084 { 1085 for (j = 0; j < ny1; j++) 1086 { 1087 index = j * nx1 + i; 1088 if (rim[index] > 150) 1089 p1p0[index]=1.0; 1090 else 1091 p1p0[index]=0.0; 1092 if (rim[index] < -150) 1093 p1n0[index]=1.0; 1094 else 1095 p1n0[index]=0.0; 1096 } 1097 }
1098 mbobra 1.30 1099 // =============== [STEP 3] =============== 1100 // smooth each of the negative and positive pixel bitmaps
1101 xudong 1.27 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0); 1102 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1103 mbobra 1.30 1104 // =============== [STEP 4] =============== 1105 // find the pixels for which p1p and p1n are both equal to 1. 1106 // this defines the polarity inversion line
1107 xudong 1.27 for (i = 0; i < nx1; i++) 1108 { 1109 for (j = 0; j < ny1; j++) 1110 { 1111 index = j * nx1 + i;
1112 mbobra 1.30 if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
1113 xudong 1.27 p1[index]=1.0; 1114 else 1115 p1[index]=0.0; 1116 } 1117 }
1118 mbobra 1.30 1119 // pad p1 with zeroes so that the gaussian colvolution in step 5 1120 // does not cut off data within hwidth of the edge 1121 1122 // step i: zero p1pad 1123 for (i = 0; i < nxp; i++) 1124 { 1125 for (j = 0; j < nyp; j++) 1126 { 1127 index = j * nxp + i; 1128 p1pad[index]=0.0; 1129 } 1130 } 1131 1132 // step ii: place p1 at the center of p1pad 1133 for (i = 0; i < nx1; i++) 1134 { 1135 for (j = 0; j < ny1; j++) 1136 { 1137 index = j * nx1 + i; 1138 index1 = (j+20) * nxp + (i+20); 1139 mbobra 1.30 p1pad[index1]=p1[index]; 1140 } 1141 } 1142 1143 // =============== [STEP 5] =============== 1144 // convolve the polarity inversion line map with a gaussian 1145 // to identify the region near the plarity inversion line 1146 // the resultant array is called pmap 1147 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0); 1148 1149 1150 // select out the nx1 x ny1 non-padded array within pmap 1151 for (i = 0; i < nx1; i++) 1152 { 1153 for (j = 0; j < ny1; j++) 1154 { 1155 index = j * nx1 + i; 1156 index1 = (j+20) * nxp + (i+20); 1157 pmapn[index]=pmap[index1]; 1158 } 1159 } 1160 mbobra 1.30 1161 // =============== [STEP 6] =============== 1162 // the R parameter is calculated
1163 xudong 1.27 for (i = 0; i < nx1; i++) 1164 { 1165 for (j = 0; j < ny1; j++) 1166 { 1167 index = j * nx1 + i;
1168 mbobra 1.32 if isnan(pmapn[index]) continue; 1169 if isnan(rim[index]) continue;
1170 mbobra 1.30 sum += pmapn[index]*abs(rim[index]);
1171 xudong 1.27 } 1172 } 1173 1174 if (sum < 1.0) 1175 Rparam = 0.0; 1176 else 1177 Rparam = log10(sum);
1178 mbobra 1.30
1179 mbobra 1.32 //printf("R_VALUE=%f\n",*Rparam); 1180
1181 xudong 1.27 free_fresize(&fresboxcar); 1182 free_fresize(&fresgauss);
1183 mbobra 1.30
1184 xudong 1.27 return 0;
1185 mbobra 1.30
1186 xudong 1.1 } 1187
1188 mbobra 1.31 /===========================================/ 1189 /* Example function 16: Lorentz force as defined in Fisher, 2012 / 1190 // 1191 // This calculation is adapted from Xudong's code 1192 // at /proj/cgem/lorentz/apps/lorentz.c 1193 1194 int computeLorentz(float bx, float by, float bz, float fx, float fy, float fz, int dims, 1195 float totfx_ptr, float totfy_ptr, float totfz_ptr, float totbsq_ptr, 1196 float epsx_ptr, float epsy_ptr, float epsz_ptr, int mask, int bitmask, 1197 float cdelt1, double rsun_ref, double rsun_obs) 1198 1199 { 1200 1201 int nx = dims[0]; 1202 int ny = dims[1]; 1203 int nxny = nxny; 1204 int j = 0; 1205 int index; 1206 double totfx = 0, totfy = 0, totfz = 0; 1207 double bsq = 0, totbsq = 0; 1208 double epsx = 0, epsy = 0, epsz = 0; 1209 mbobra 1.31 double area = cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 1210 double k_h = -1.0 area / (4. * PI) / 1.0e20; 1211 double k_z = area / (8. * PI) / 1.0e20; 1212 1213 if (nx <= 0 \|\| ny <= 0) return 1; 1214 1215 for (int i = 0; i < nxny; i++) 1216 { 1217 if ( mask[i] < 70 \|\| bitmask[i] < 30 ) continue; 1218 if isnan(bx[i]) continue; 1219 if isnan(by[i]) continue; 1220 if isnan(bz[i]) continue;
1221 mbobra 1.32 fx[i] = bx[i] * bz[i] * k_h; 1222 fy[i] = by[i] * bz[i] * k_h;
1223 mbobra 1.31 fz[i] = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z; 1224 bsq = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i]; 1225 totfx += fx[i]; totfy += fy[i]; totfz += fz[i]; 1226 totbsq += bsq; 1227 } 1228 1229 totfx_ptr = totfx; 1230 totfy_ptr = totfy; 1231 totfz_ptr = totfz; 1232 totbsq_ptr = totbsq; 1233 epsx_ptr = (totfx / k_h) / totbsq; 1234 epsy_ptr = (totfy / k_h) / totbsq; 1235 *epsz_ptr = (totfz / k_z) / totbsq; 1236
1237 mbobra 1.32 //printf("TOTBSQ=%f\n",*totbsq_ptr); 1238
1239 mbobra 1.31 return 0; 1240 1241 } 1242
1243 mbobra 1.38 /===========================================/ 1244 1245 /* Example function 17: Compute total unsigned flux in units of G/cm^2 on the LOS field / 1246 1247 // To compute the unsigned flux, we simply calculate 1248 // flux = surface integral [(vector LOS) dot (normal vector)], 1249 // = surface integral [(magnitude LOS)(magnitude normal)(cos theta)]. 1250 // However, since the field is radial, we will assume cos theta = 1. 1251 // Therefore the pixels only need to be corrected for the projection. 1252 1253 // To convert G to Gcm^2, simply multiply by the number of square centimeters per pixel. 1254 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 1255 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2 1256 // =Gausscm^2 1257 1258 int computeAbsFlux_los(float los, int dims, float absFlux, 1259 float mean_vf_ptr, float count_mask_ptr, 1260 int bitmask, float cdelt1, double rsun_ref, double rsun_obs) 1261 1262 { 1263 1264 mbobra 1.38 int nx = dims[0]; 1265 int ny = dims[1]; 1266 int i = 0; 1267 int j = 0; 1268 int count_mask = 0; 1269 double sum = 0.0; 1270 absFlux = 0.0; 1271 mean_vf_ptr = 0.0; 1272 1273 1274 if (nx <= 0 \|\| ny <= 0) return 1; 1275 1276 for (i = 0; i < nx; i++) 1277 { 1278 for (j = 0; j < ny; j++) 1279 { 1280 if ( bitmask[j nx + i] < 30 ) continue; 1281 if isnan(los[j * nx + i]) continue; 1282 sum += (fabs(los[j * nx + i])); 1283 count_mask++; 1284 } 1285 mbobra 1.38 } 1286 1287 mean_vf_ptr = sumcdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 1288 count_mask_ptr = count_mask; 1289 1290 return 0; 1291 } 1292 1293 /===========================================/ 1294 /* Example function 18: Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) / 1295 1296 int computeLOSderivative(float los, int dims, float mean_derivative_los_ptr, int bitmask, float derx_los, float dery_los) 1297 { 1298 1299 int nx = dims[0]; 1300 int ny = dims[1]; 1301 int i = 0; 1302 int j = 0; 1303 int count_mask = 0; 1304 double sum = 0.0; 1305 mean_derivative_los_ptr = 0.0; 1306 mbobra 1.38 1307 if (nx <= 0 \|\| ny <= 0) return 1; 1308 1309 /* brute force method of calculating the derivative (no consideration for edges) / 1310 for (i = 1; i <= nx-2; i++) 1311 { 1312 for (j = 0; j <= ny-1; j++) 1313 { 1314 derx_los[j nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])0.5; 1315 } 1316 } 1317 1318 / brute force method of calculating the derivative (no consideration for edges) / 1319 for (i = 0; i <= nx-1; i++) 1320 { 1321 for (j = 1; j <= ny-2; j++) 1322 { 1323 dery_los[j nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])0.5; 1324 } 1325 } 1326 1327 mbobra 1.38 / consider the edges for the arrays that contribute to the variable "sum" in the computation below. 1328 ignore the edges for the error terms as those arrays have been initialized to zero. 1329 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./ 1330 i=0; 1331 for (j = 0; j <= ny-1; j++) 1332 { 1333 derx_los[j nx + i] = ( (-3los[j nx + i]) + (4los[j nx + (i+1)]) - (los[j * nx + (i+2)]) )0.5; 1334 } 1335 1336 i=nx-1; 1337 for (j = 0; j <= ny-1; j++) 1338 { 1339 derx_los[j nx + i] = ( (3los[j nx + i]) + (-4los[j nx + (i-1)]) - (-los[j * nx + (i-2)]) )0.5; 1340 } 1341 1342 j=0; 1343 for (i = 0; i <= nx-1; i++) 1344 { 1345 dery_los[j nx + i] = ( (-3los[jnx + i]) + (4los[(j+1) nx + i]) - (los[(j+2) * nx + i]) )0.5; 1346 } 1347 1348 mbobra 1.38 j=ny-1; 1349 for (i = 0; i <= nx-1; i++) 1350 { 1351 dery_los[j nx + i] = ( (3los[j nx + i]) + (-4los[(j-1) nx + i]) - (-los[(j-2) * nx + i]) )0.5; 1352 } 1353 1354 1355 for (i = 0; i <= nx-1; i++) 1356 { 1357 for (j = 0; j <= ny-1; j++) 1358 { 1359 if ( bitmask[j nx + i] < 30 ) continue; 1360 if ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) continue; 1361 if isnan(los[j * nx + i]) continue; 1362 if isnan(los[(j+1) * nx + i]) continue; 1363 if isnan(los[(j-1) * nx + i]) continue; 1364 if isnan(los[j * nx + i-1]) continue; 1365 if isnan(los[j * nx + i+1]) continue; 1366 if isnan(derx_los[j * nx + i]) continue; 1367 if isnan(dery_los[j * nx + i]) continue; 1368 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 / 1369 mbobra 1.38 count_mask++; 1370 } 1371 } 1372 1373 mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)(ny-2)) if shape of count_mask = shape of magnetogram 1374 //printf("mean_derivative_los_ptr=%f\n",mean_derivative_los_ptr); 1375 1376 return 0; 1377 } 1378
1379 xudong 1.1 /==================KEIJI'S CODE =========================/ 1380 1381 // #include <omp.h> 1382 #include <math.h> 1383 1384 void greenpot(float bx, float by, float *bz, int nnx, int nny) 1385 {
1386 xudong 1.27 /* local workings / 1387 int inx, iny, i, j, n; 1388 / local array / 1389 float pfpot, rdist; 1390 pfpot=(float )malloc(sizeof(float) nnxnny); 1391 rdist=(float )malloc(sizeof(float) nnxnny); 1392 float bztmp; 1393 bztmp=(float )malloc(sizeof(float) nnxnny); 1394 / make nan / 1395 // unsigned long long llnan = 0x7ff0000000000000; 1396 // float NAN = (float)(llnan); 1397 1398 // #pragma omp parallel for private (inx) 1399 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnxiny+inx] = 0.0;}} 1400 // #pragma omp parallel for private (inx) 1401 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnxiny+inx] = 0.0;}} 1402 // #pragma omp parallel for private (inx) 1403 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnxiny+inx] = 0.0;}} 1404 // #pragma omp parallel for private (inx) 1405 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnxiny+inx] = 0.0;}} 1406 // #pragma omp parallel for private (inx) 1407 xudong 1.27 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1408 { 1409 float val0 = bz[nnxiny + inx]; 1410 if (isnan(val0)){bztmp[nnxiny + inx] = 0.0;}else{bztmp[nnxiny + inx] = val0;} 1411 }} 1412 1413 // dz is the monopole depth 1414 float dz = 0.001; 1415 1416 // #pragma omp parallel for private (inx) 1417 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1418 { 1419 float rdd, rdd1, rdd2; 1420 float r; 1421 rdd1 = (float)(inx); 1422 rdd2 = (float)(iny); 1423 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz; 1424 rdist[nnxiny+inx] = 1.0/sqrt(rdd); 1425 }} 1426 1427 int iwindow; 1428 xudong 1.27 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;} 1429 float rwindow; 1430 rwindow = (float)(iwindow); 1431 rwindow = rwindow rwindow + 0.01; // must be of square 1432 1433 rwindow = 1.0e2; // limit the window size to be 10. 1434 1435 rwindow = sqrt(rwindow); 1436 iwindow = (int)(rwindow); 1437 1438 // #pragma omp parallel for private(inx) 1439 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
1440 xudong 1.1 {
1441 xudong 1.27 float val0 = bz[nnxiny + inx]; 1442 if (isnan(val0)) 1443 { 1444 pfpot[nnxiny + inx] = 0.0; // hmmm.. NAN; 1445 } 1446 else 1447 { 1448 float sum; 1449 sum = 0.0; 1450 int j2, i2; 1451 int j2s, j2e, i2s, i2e; 1452 j2s = iny - iwindow; 1453 j2e = iny + iwindow; 1454 if (j2s < 0){j2s = 0;} 1455 if (j2e > nny){j2e = nny;} 1456 i2s = inx - iwindow; 1457 i2e = inx + iwindow; 1458 if (i2s < 0){i2s = 0;}
1459 mbobra 1.37 if (i2e > nnx){i2e = nnx;} 1460 1461 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++) 1462 {
1463 xudong 1.27 float val1 = bztmp[nnx*j2 + i2];
1464 mbobra 1.37 float rr, r1, r2; 1465 // r1 = (float)(i2 - inx); 1466 // r2 = (float)(j2 - iny); 1467 // rr = r1r1 + r2r2; 1468 // if (rr < rwindow) 1469 // {
1470 xudong 1.27 int di, dj; 1471 di = abs(i2 - inx); 1472 dj = abs(j2 - iny); 1473 sum = sum + val1 * rdist[nnx * dj + di] * dz;
1474 mbobra 1.37 // }
1475 xudong 1.27 } } 1476 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition. 1477 } 1478 } } // end of OpenMP parallelism 1479 1480 // #pragma omp parallel for private(inx) 1481 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1482 xudong 1.1 {
1483 xudong 1.27 bx[nnxiny + inx] = -(pfpot[nnxiny + (inx+1)]-pfpot[nnxiny + (inx-1)]) 0.5; 1484 by[nnxiny + inx] = -(pfpot[nnx(iny+1) + inx]-pfpot[nnx(iny-1) + inx]) 0.5; 1485 } } // end of OpenMP parallelism
1486 mbobra 1.37
1487 xudong 1.27 free(rdist); 1488 free(pfpot); 1489 free(bztmp);
1490 xudong 1.1 } // end of void func. greenpot 1491 1492 1493 /===========END OF KEIJI'S CODE =========================/
1494 mbobra 1.14 1495 char *sw_functions_version() // Returns CVS version of sw_functions.c 1496 {
1497 xudong 1.27 return strdup("$Id");
1498 mbobra 1.14 } 1499
1500 xudong 1.1 /* ---------------- end of this file ----------------*/

Karen Tian