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 meanshear_angle_err_ptr = (sqrt(err)/count_mask)(180./PI); 1026 1027 / The area here is a fractional area -- the % of the total area. This has no error associated with it. / 1028 area_w_shear_gt_45ptr = (count_mask/(count))(100.0); 1029 1030 //printf("MEANSHR=%f\n",meanshear_angleptr); 1031 //printf("MEANSHR_err=%f\n",meanshear_angle_err_ptr); 1032 //printf("SHRGT45=%f\n",area_w_shear_gt_45ptr); 1033 1034 return 0; 1035 }
1036 xudong 1.1
1037 xudong 1.27 /===========================================/
1038 mbobra 1.29 /* Example function 15: R parameter as defined in Schrijver, 2007 / 1039 // 1040 // Note that there is a restriction on the function fsample() 1041 // If the following occurs: 1042 // nx_out > floor((ny_in-1)/scale + 1) 1043 // ny_out > floor((ny_in-1)/scale + 1), 1044 // where n_out are the dimensions of the output array and n*_in 1045 // are the dimensions of the input array, fsample() will usually result 1046 // in a segfault (though not always, depending on how the segfault was accessed.) 1047
1048 xudong 1.27 int computeR(float bz_err, float los, int dims, float Rparam, float cdelt1, 1049 float rim, float p1p0, float p1n0, float p1p, float p1n, float p1,
1050 mbobra 1.30 float pmap, int nx1, int ny1, 1051 int scale, float p1pad, int nxp, int nyp, float *pmapn) 1052
1053 mbobra 1.29 {
1054 xudong 1.27 int nx = dims[0]; 1055 int ny = dims[1]; 1056 int i = 0; 1057 int j = 0;
1058 mbobra 1.30 int index, index1;
1059 xudong 1.27 double sum = 0.0; 1060 double err = 0.0; 1061 *Rparam = 0.0; 1062 struct fresize_struct fresboxcar, fresgauss;
1063 mbobra 1.30 struct fint_struct fints;
1064 xudong 1.27 float sigma = 10.0/2.3548; 1065
1066 mbobra 1.30 // set up convolution kernels
1067 xudong 1.27 init_fresize_boxcar(&fresboxcar,1,1); 1068 init_fresize_gaussian(&fresgauss,sigma,20,1);
1069 mbobra 1.29
1070 mbobra 1.30 // =============== [STEP 1] =============== 1071 // bin the line-of-sight magnetogram down by a factor of scale
1072 xudong 1.27 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
1073 mbobra 1.29
1074 mbobra 1.30 // =============== [STEP 2] =============== 1075 // identify positive and negative pixels greater than +/- 150 gauss 1076 // and label those pixels with a 1.0 in arrays p1p0 and p1n0
1077 xudong 1.27 for (i = 0; i < nx1; i++) 1078 { 1079 for (j = 0; j < ny1; j++) 1080 { 1081 index = j * nx1 + i; 1082 if (rim[index] > 150) 1083 p1p0[index]=1.0; 1084 else 1085 p1p0[index]=0.0; 1086 if (rim[index] < -150) 1087 p1n0[index]=1.0; 1088 else 1089 p1n0[index]=0.0; 1090 } 1091 }
1092 mbobra 1.30 1093 // =============== [STEP 3] =============== 1094 // smooth each of the negative and positive pixel bitmaps
1095 xudong 1.27 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0); 1096 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1097 mbobra 1.30 1098 // =============== [STEP 4] =============== 1099 // find the pixels for which p1p and p1n are both equal to 1. 1100 // this defines the polarity inversion line
1101 xudong 1.27 for (i = 0; i < nx1; i++) 1102 { 1103 for (j = 0; j < ny1; j++) 1104 { 1105 index = j * nx1 + i;
1106 mbobra 1.30 if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
1107 xudong 1.27 p1[index]=1.0; 1108 else 1109 p1[index]=0.0; 1110 } 1111 }
1112 mbobra 1.30 1113 // pad p1 with zeroes so that the gaussian colvolution in step 5 1114 // does not cut off data within hwidth of the edge 1115 1116 // step i: zero p1pad 1117 for (i = 0; i < nxp; i++) 1118 { 1119 for (j = 0; j < nyp; j++) 1120 { 1121 index = j * nxp + i; 1122 p1pad[index]=0.0; 1123 } 1124 } 1125 1126 // step ii: place p1 at the center of p1pad 1127 for (i = 0; i < nx1; i++) 1128 { 1129 for (j = 0; j < ny1; j++) 1130 { 1131 index = j * nx1 + i; 1132 index1 = (j+20) * nxp + (i+20); 1133 mbobra 1.30 p1pad[index1]=p1[index]; 1134 } 1135 } 1136 1137 // =============== [STEP 5] =============== 1138 // convolve the polarity inversion line map with a gaussian 1139 // to identify the region near the plarity inversion line 1140 // the resultant array is called pmap 1141 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0); 1142 1143 1144 // select out the nx1 x ny1 non-padded array within pmap 1145 for (i = 0; i < nx1; i++) 1146 { 1147 for (j = 0; j < ny1; j++) 1148 { 1149 index = j * nx1 + i; 1150 index1 = (j+20) * nxp + (i+20); 1151 pmapn[index]=pmap[index1]; 1152 } 1153 } 1154 mbobra 1.30 1155 // =============== [STEP 6] =============== 1156 // the R parameter is calculated
1157 xudong 1.27 for (i = 0; i < nx1; i++) 1158 { 1159 for (j = 0; j < ny1; j++) 1160 { 1161 index = j * nx1 + i;
1162 mbobra 1.32 if isnan(pmapn[index]) continue; 1163 if isnan(rim[index]) continue;
1164 mbobra 1.30 sum += pmapn[index]*abs(rim[index]);
1165 xudong 1.27 } 1166 } 1167 1168 if (sum < 1.0) 1169 Rparam = 0.0; 1170 else 1171 Rparam = log10(sum);
1172 mbobra 1.30
1173 mbobra 1.32 //printf("R_VALUE=%f\n",*Rparam); 1174
1175 xudong 1.27 free_fresize(&fresboxcar); 1176 free_fresize(&fresgauss);
1177 mbobra 1.30
1178 xudong 1.27 return 0;
1179 mbobra 1.30
1180 xudong 1.1 } 1181
1182 mbobra 1.31 /===========================================/ 1183 /* Example function 16: Lorentz force as defined in Fisher, 2012 / 1184 // 1185 // This calculation is adapted from Xudong's code 1186 // at /proj/cgem/lorentz/apps/lorentz.c 1187 1188 int computeLorentz(float bx, float by, float bz, float fx, float fy, float fz, int dims, 1189 float totfx_ptr, float totfy_ptr, float totfz_ptr, float totbsq_ptr, 1190 float epsx_ptr, float epsy_ptr, float epsz_ptr, int mask, int bitmask, 1191 float cdelt1, double rsun_ref, double rsun_obs) 1192 1193 { 1194 1195 int nx = dims[0]; 1196 int ny = dims[1]; 1197 int nxny = nxny; 1198 int j = 0; 1199 int index; 1200 double totfx = 0, totfy = 0, totfz = 0; 1201 double bsq = 0, totbsq = 0; 1202 double epsx = 0, epsy = 0, epsz = 0; 1203 mbobra 1.31 double area = cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 1204 double k_h = -1.0 area / (4. * PI) / 1.0e20; 1205 double k_z = area / (8. * PI) / 1.0e20; 1206 1207 if (nx <= 0 \|\| ny <= 0) return 1; 1208 1209 for (int i = 0; i < nxny; i++) 1210 { 1211 if ( mask[i] < 70 \|\| bitmask[i] < 30 ) continue; 1212 if isnan(bx[i]) continue; 1213 if isnan(by[i]) continue; 1214 if isnan(bz[i]) continue;
1215 mbobra 1.32 fx[i] = bx[i] * bz[i] * k_h; 1216 fy[i] = by[i] * bz[i] * k_h;
1217 mbobra 1.31 fz[i] = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z; 1218 bsq = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i]; 1219 totfx += fx[i]; totfy += fy[i]; totfz += fz[i]; 1220 totbsq += bsq; 1221 } 1222 1223 totfx_ptr = totfx; 1224 totfy_ptr = totfy; 1225 totfz_ptr = totfz; 1226 totbsq_ptr = totbsq; 1227 epsx_ptr = (totfx / k_h) / totbsq; 1228 epsy_ptr = (totfy / k_h) / totbsq; 1229 *epsz_ptr = (totfz / k_z) / totbsq; 1230
1231 mbobra 1.32 //printf("TOTBSQ=%f\n",*totbsq_ptr); 1232
1233 mbobra 1.31 return 0; 1234 1235 } 1236
1237 mbobra 1.38 /===========================================/ 1238 1239 /* Example function 17: Compute total unsigned flux in units of G/cm^2 on the LOS field / 1240 1241 // To compute the unsigned flux, we simply calculate 1242 // flux = surface integral [(vector LOS) dot (normal vector)], 1243 // = surface integral [(magnitude LOS)(magnitude normal)(cos theta)]. 1244 // However, since the field is radial, we will assume cos theta = 1. 1245 // Therefore the pixels only need to be corrected for the projection. 1246 1247 // To convert G to Gcm^2, simply multiply by the number of square centimeters per pixel. 1248 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 1249 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2 1250 // =Gausscm^2 1251 1252 int computeAbsFlux_los(float los, int dims, float absFlux, 1253 float mean_vf_ptr, float count_mask_ptr, 1254 int bitmask, float cdelt1, double rsun_ref, double rsun_obs) 1255 1256 { 1257 1258 mbobra 1.38 int nx = dims[0]; 1259 int ny = dims[1]; 1260 int i = 0; 1261 int j = 0; 1262 int count_mask = 0; 1263 double sum = 0.0; 1264 absFlux = 0.0; 1265 mean_vf_ptr = 0.0; 1266 1267 1268 if (nx <= 0 \|\| ny <= 0) return 1; 1269 1270 for (i = 0; i < nx; i++) 1271 { 1272 for (j = 0; j < ny; j++) 1273 { 1274 if ( bitmask[j nx + i] < 30 ) continue; 1275 if isnan(los[j * nx + i]) continue; 1276 sum += (fabs(los[j * nx + i])); 1277 count_mask++; 1278 } 1279 mbobra 1.38 } 1280 1281 mean_vf_ptr = sumcdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 1282 count_mask_ptr = count_mask; 1283 1284 return 0; 1285 } 1286 1287 /===========================================/ 1288 /* Example function 18: Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) / 1289 1290 int computeLOSderivative(float los, int dims, float mean_derivative_los_ptr, int bitmask, float derx_los, float dery_los) 1291 { 1292 1293 int nx = dims[0]; 1294 int ny = dims[1]; 1295 int i = 0; 1296 int j = 0; 1297 int count_mask = 0; 1298 double sum = 0.0; 1299 mean_derivative_los_ptr = 0.0; 1300 mbobra 1.38 1301 if (nx <= 0 \|\| ny <= 0) return 1; 1302 1303 /* brute force method of calculating the derivative (no consideration for edges) / 1304 for (i = 1; i <= nx-2; i++) 1305 { 1306 for (j = 0; j <= ny-1; j++) 1307 { 1308 derx_los[j nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])0.5; 1309 } 1310 } 1311 1312 / brute force method of calculating the derivative (no consideration for edges) / 1313 for (i = 0; i <= nx-1; i++) 1314 { 1315 for (j = 1; j <= ny-2; j++) 1316 { 1317 dery_los[j nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])0.5; 1318 } 1319 } 1320 1321 mbobra 1.38 / consider the edges for the arrays that contribute to the variable "sum" in the computation below. 1322 ignore the edges for the error terms as those arrays have been initialized to zero. 1323 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./ 1324 i=0; 1325 for (j = 0; j <= ny-1; j++) 1326 { 1327 derx_los[j nx + i] = ( (-3los[j nx + i]) + (4los[j nx + (i+1)]) - (los[j * nx + (i+2)]) )0.5; 1328 } 1329 1330 i=nx-1; 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 j=0; 1337 for (i = 0; i <= nx-1; i++) 1338 { 1339 dery_los[j nx + i] = ( (-3los[jnx + i]) + (4los[(j+1) nx + i]) - (los[(j+2) * nx + i]) )0.5; 1340 } 1341 1342 mbobra 1.38 j=ny-1; 1343 for (i = 0; i <= nx-1; i++) 1344 { 1345 dery_los[j nx + i] = ( (3los[j nx + i]) + (-4los[(j-1) nx + i]) - (-los[(j-2) * nx + i]) )0.5; 1346 } 1347 1348 1349 for (i = 0; i <= nx-1; i++) 1350 { 1351 for (j = 0; j <= ny-1; j++) 1352 { 1353 if ( bitmask[j nx + i] < 30 ) continue; 1354 if ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) continue; 1355 if isnan(los[j * nx + i]) continue; 1356 if isnan(los[(j+1) * nx + i]) continue; 1357 if isnan(los[(j-1) * nx + i]) continue; 1358 if isnan(los[j * nx + i-1]) continue; 1359 if isnan(los[j * nx + i+1]) continue; 1360 if isnan(derx_los[j * nx + i]) continue; 1361 if isnan(dery_los[j * nx + i]) continue; 1362 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 / 1363 mbobra 1.38 count_mask++; 1364 } 1365 } 1366 1367 mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)(ny-2)) if shape of count_mask = shape of magnetogram 1368 //printf("mean_derivative_los_ptr=%f\n",mean_derivative_los_ptr); 1369 1370 return 0; 1371 } 1372
1373 xudong 1.1 /==================KEIJI'S CODE =========================/ 1374 1375 // #include <omp.h> 1376 #include <math.h> 1377 1378 void greenpot(float bx, float by, float *bz, int nnx, int nny) 1379 {
1380 xudong 1.27 /* local workings / 1381 int inx, iny, i, j, n; 1382 / local array / 1383 float pfpot, rdist; 1384 pfpot=(float )malloc(sizeof(float) nnxnny); 1385 rdist=(float )malloc(sizeof(float) nnxnny); 1386 float bztmp; 1387 bztmp=(float )malloc(sizeof(float) nnxnny); 1388 / make nan / 1389 // unsigned long long llnan = 0x7ff0000000000000; 1390 // float NAN = (float)(llnan); 1391 1392 // #pragma omp parallel for private (inx) 1393 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnxiny+inx] = 0.0;}} 1394 // #pragma omp parallel for private (inx) 1395 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnxiny+inx] = 0.0;}} 1396 // #pragma omp parallel for private (inx) 1397 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnxiny+inx] = 0.0;}} 1398 // #pragma omp parallel for private (inx) 1399 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnxiny+inx] = 0.0;}} 1400 // #pragma omp parallel for private (inx) 1401 xudong 1.27 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1402 { 1403 float val0 = bz[nnxiny + inx]; 1404 if (isnan(val0)){bztmp[nnxiny + inx] = 0.0;}else{bztmp[nnxiny + inx] = val0;} 1405 }} 1406 1407 // dz is the monopole depth 1408 float dz = 0.001; 1409 1410 // #pragma omp parallel for private (inx) 1411 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1412 { 1413 float rdd, rdd1, rdd2; 1414 float r; 1415 rdd1 = (float)(inx); 1416 rdd2 = (float)(iny); 1417 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz; 1418 rdist[nnxiny+inx] = 1.0/sqrt(rdd); 1419 }} 1420 1421 int iwindow; 1422 xudong 1.27 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;} 1423 float rwindow; 1424 rwindow = (float)(iwindow); 1425 rwindow = rwindow rwindow + 0.01; // must be of square 1426 1427 rwindow = 1.0e2; // limit the window size to be 10. 1428 1429 rwindow = sqrt(rwindow); 1430 iwindow = (int)(rwindow); 1431 1432 // #pragma omp parallel for private(inx) 1433 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
1434 xudong 1.1 {
1435 xudong 1.27 float val0 = bz[nnxiny + inx]; 1436 if (isnan(val0)) 1437 { 1438 pfpot[nnxiny + inx] = 0.0; // hmmm.. NAN; 1439 } 1440 else 1441 { 1442 float sum; 1443 sum = 0.0; 1444 int j2, i2; 1445 int j2s, j2e, i2s, i2e; 1446 j2s = iny - iwindow; 1447 j2e = iny + iwindow; 1448 if (j2s < 0){j2s = 0;} 1449 if (j2e > nny){j2e = nny;} 1450 i2s = inx - iwindow; 1451 i2e = inx + iwindow; 1452 if (i2s < 0){i2s = 0;}
1453 mbobra 1.37 if (i2e > nnx){i2e = nnx;} 1454 1455 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++) 1456 {
1457 xudong 1.27 float val1 = bztmp[nnx*j2 + i2];
1458 mbobra 1.37 float rr, r1, r2; 1459 // r1 = (float)(i2 - inx); 1460 // r2 = (float)(j2 - iny); 1461 // rr = r1r1 + r2r2; 1462 // if (rr < rwindow) 1463 // {
1464 xudong 1.27 int di, dj; 1465 di = abs(i2 - inx); 1466 dj = abs(j2 - iny); 1467 sum = sum + val1 * rdist[nnx * dj + di] * dz;
1468 mbobra 1.37 // }
1469 xudong 1.27 } } 1470 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition. 1471 } 1472 } } // end of OpenMP parallelism 1473 1474 // #pragma omp parallel for private(inx) 1475 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1476 xudong 1.1 {
1477 xudong 1.27 bx[nnxiny + inx] = -(pfpot[nnxiny + (inx+1)]-pfpot[nnxiny + (inx-1)]) 0.5; 1478 by[nnxiny + inx] = -(pfpot[nnx(iny+1) + inx]-pfpot[nnx(iny-1) + inx]) 0.5; 1479 } } // end of OpenMP parallelism
1480 mbobra 1.37
1481 xudong 1.27 free(rdist); 1482 free(pfpot); 1483 free(bztmp);
1484 xudong 1.1 } // end of void func. greenpot 1485 1486 1487 /===========END OF KEIJI'S CODE =========================/
1488 mbobra 1.14 1489 char *sw_functions_version() // Returns CVS version of sw_functions.c 1490 {
1491 xudong 1.27 return strdup("$Id");
1492 mbobra 1.14 } 1493
1494 xudong 1.1 /* ---------------- end of this file ----------------*/

Karen Tian