1 mbobra 1.1 /* sharp_functions.c */
2 #define PI (3.141592653589793)
3
4 /*===========================================*/
5
6 /* Example function 1: Compute total unsigned flux in units of G/cm^2, Mean Vertical Field */
7
8 // To convert G to G/cm^2, simply multiply by the number of square centimeters per pixel.
9 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
10 // Gauss(CDELT1)(RSUN_REF/RSUN_OBS)(100.)
11 // = Gauss(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
12 // = Gauss(1.30501e15)
13
14 int computeAbsFlux(float *bz, int *dims, float *absFlux, float *mean_vf_ptr, int *mask)
15 {
16
17 int nx = dims[0], ny = dims[1];
18 int i, j, count_mask=0;
19 double sum=0.0;
20
21 if (nx <= 0 || ny <= 0) return 1;
22 mbobra 1.1
23 *absFlux = 0.0;
24 *mean_vf_ptr =0.0;
25
26 for (j = 0; j < ny; j++)
27 {
28 for (i = 0; i < nx; i++)
29 {
30 if (mask[j * nx + i] <= 1) continue;
31 sum += (fabs(bz[j * nx + i]));
32 count_mask++;
33 }
34 }
35
36 printf("nx=%d,ny=%d,count_mask=%d,sum=%f\n",nx,ny,count_mask,sum);
37 *mean_vf_ptr = sum*(1.30501e15);
38 return 0;
39 }
40
41 /*===========================================*/
42 /* Example function 2: Calculate Bh in units of Gauss */
43 mbobra 1.1 // Native units of Bh are Gauss
44
45 int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_hf_ptr, int *mask)
46 {
47
48 int nx = dims[0], ny = dims[1];
49 int i, j, count_mask=0;
50 float sum=0.0;
51 *mean_hf_ptr =0.0;
52
53 if (nx <= 0 || ny <= 0) return 1;
54
55 for (j = 0; j < ny; j++)
56 {
57 for (i = 0; i < nx; i++)
58 {
59 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
60 sum += bh[j * nx + i];
61 count_mask++;
62 }
63 }
64 mbobra 1.1
65 *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
66 printf("*mean_hf_ptr=%f\n",*mean_hf_ptr);
67 return 0;
68 }
69
70 /*===========================================*/
71
72 /* Example function 3: Calculate Gamma in units of degrees
73 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
74
75 this is wrong, also the divide by nx*ny is wrong */
76
77 int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_gamma_ptr, int *mask)
78 {
79 int nx = dims[0], ny = dims[1];
80 int i, j, count_mask=0;
81
82 if (nx <= 0 || ny <= 0) return 1;
83
84 *mean_gamma_ptr=0.0;
85 mbobra 1.1 float sum=0.0;
86 int count=0;
87
88 for (i = 0; i < nx; i++)
89 {
90 for (j = 0; j < ny; j++)
91 {
92 if (bh[j * nx + i] > 100)
93 {
94 if (mask[j * nx + i] <= 1) continue;
95 sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
96 count_mask++;
97 }
98 }
99 }
100
101 *mean_gamma_ptr = sum/count_mask;
102 printf("*mean_gamma_ptr=%f\n",*mean_gamma_ptr);
103 return 0;
104 }
105
106 mbobra 1.1 /*===========================================*/
107
108 /* Example function 4: Calculate B_Total*/
109 // Native units of B_Total are in gauss
110
111 int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask)
112 {
113
114 int nx = dims[0], ny = dims[1];
115 int i, j, count_mask=0;
116
117 if (nx <= 0 || ny <= 0) return 1;
118
119 for (i = 0; i < nx; i++)
120 {
121 for (j = 0; j < ny; j++)
122 {
123 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]);
124 }
125 }
126 return 0;
127 mbobra 1.1 }
128
129 /*===========================================*/
130
131 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
132
133 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask)
134 {
135
136 int nx = dims[0], ny = dims[1];
137 int i, j, count_mask=0;
138
139 if (nx <= 0 || ny <= 0) return 1;
140
141 *mean_derivative_btotal_ptr = 0.0;
142 float derx, dery, sum = 0.0;
143
144 for (i = 1; i < nx-1; i++)
145 {
146 for (j = 1; j < ny-1; j++)
147 {
148 mbobra 1.1 if (mask[j * nx + i] <= 1) continue;
149 /* Missing a (*0.5) */
150 derx = bt[j * nx + i+1] - bt[j * nx + i-1];
151 dery = bt[(j+1) * nx + i] - bt[(j-1) * nx + i];
152 sum += sqrt(derx*derx + dery*dery);
153 count_mask++;
154 }
155 }
156
157 *mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
158 printf("*mean_derivative_btotal_ptr=%f\n",*mean_derivative_btotal_ptr);
159 return 0;
160 }
161
162 /*===========================================*/
163 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
164
165 int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask)
166 {
167
168 int nx = dims[0], ny = dims[1];
169 mbobra 1.1 int i, j, count_mask=0;
170
171 if (nx <= 0 || ny <= 0) return 1;
172
173 *mean_derivative_bh_ptr = 0.0;
174 float derx, dery, sum = 0.0;
175
176 for (i = 1; i < nx-1; i++)
177 {
178 for (j = 1; j < ny-1; j++)
179 {
180 if (mask[j * nx + i] <= 1) continue;
181 /* Missing a (*0.5) */
182 derx = bh[j * nx + i+1] - bh[j * nx + i-1];
183 //printf("derx=%f\n",derx);
184 dery = bh[(j+1) * nx + i] - bh[(j-1) * nx + i];
185 // printf("dery=%f\n",dery);
186 sum += sqrt(derx*derx + dery*dery);
187 //printf("sum=%f\n",sum);
188 count_mask++;
189 //printf("count_mask=%d\n",count_mask);
190 mbobra 1.1 }
191 }
192
193 *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
194 printf("*mean_derivative_bh_ptr=%f,nx=%d,ny=%d,sum=%f\n",*mean_derivative_bh_ptr,nx,ny,sum);
195 return 0;
196 }
197
198 /*===========================================*/
199 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
200
201 int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask)
202 {
203
204 int nx = dims[0], ny = dims[1];
205 int i, j, count_mask=0;
206
207 if (nx <= 0 || ny <= 0) return 1;
208
209 *mean_derivative_bz_ptr = 0.0;
210 float derx, dery, sum = 0.0;
211 mbobra 1.1
212 for (i = 1; i < nx-1; i++)
213 {
214 for (j = 1; j < ny-1; j++)
215 {
216 if (mask[j * nx + i] <= 1) continue;
217 /* Missing a (*0.5) */
218 derx = bz[j * nx + i+1] - bz[j * nx + i-1];
219 dery = bz[(j+1) * nx + i] - bz[(j-1) * nx + i];
220 sum += sqrt(derx*derx + dery*dery);
221 count_mask++;
222 }
223 }
224
225 *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
226
227 return 0;
228 }
229
230 /*===========================================*/
231 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */
232 mbobra 1.1
233 // In discretized space like data pixels,
234 // the current (or curl of B) is calculated as
235 // the integration of the field Bx and By along
236 // the circumference of the data pixel divided by the area of the pixel.
237 // One form of differencing (a word for the differential operator
238 // in the discretized space) of the curl is expressed as
239 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2
240 // +dy * (By(i+1,j)+By(i,j)) / 2
241 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
242 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
243 //
244 //
245 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
246 // one must perform the following unit conversions:
247 // (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
248 // (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),
249 // where a Tesla is represented as a Newton/Ampere*meter.
250 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
251 // In that case, we would have the following:
252 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
253 mbobra 1.1 // jz * (35.0)
254 //
255 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
256 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(722500)(722500), or
257 // (Gauss/pix)(1.81*10^10)
258 // (Gauss/pix)(18100000000.)
259
260 int computeJz(float *bx, float *by, int *dims, float *jz, float *mean_jz_ptr, float *us_i_ptr, int *mask)
261 {
262
263 int nx = dims[0], ny = dims[1];
264 int i, j, count_mask=0;
265
266 if (nx <= 0 || ny <= 0) return 1;
267
268 *mean_jz_ptr = 0.0;
269 float derx, dery, curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
270
271 for (i = 1; i < nx-1; i++)
272 {
273 for (j = 1; j < ny-1; j++)
274 mbobra 1.1 {
275 if (mask[j * nx + i] <= 1) continue;
276 derx = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
277 dery = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
278 curl += (derx-dery)*(35.0); /* curl is in units of mA/m^2 */
279 jz[j * nx + i] = (derx-dery); /* jz is in units of Gauss per pixel */
280 us_i += fabs(derx-dery)*(18100000000.) ; /* us_i is in units of A */
281 count_mask++;
282 }
283 }
284
285 mean_curl = (curl/count_mask);
286 printf("mean_curl=%f\n",mean_curl);
287 *mean_jz_ptr = curl/(count_mask);
288 printf("count_mask=%d\n",count_mask);
289 *us_i_ptr = (us_i);
290 return 0;
291 }
292
293
294
295 mbobra 1.1 /*===========================================*/
296 /* Example function 9: Twist Parameter, alpha */
297
298 // The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm
299 // The units of Jz are in G/pix; the units of Bz are in G.
300 // Therefore, the units of Jz/Bz = (pix/arcsec)(arcsec/meter)(meter/Mm), or
301 // Jz/Bz = (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6).
302 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
303 // In that case, we would have the following:
304 // (Gauss/pix)(1/0.5)(1/722500)(10^6), or
305 // (Gauss/pix)*2.7
306
307 int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask)
308 {
309 int nx = dims[0], ny = dims[1];
310 int i, j, count_mask=0;
311
312 if (nx <= 0 || ny <= 0) return 1;
313
314 *mean_alpha_ptr = 0.0;
315 float aa, bb, cc, bznew, alpha2, sum=0.0;
316 mbobra 1.1
317 for (i = 1; i < nx-1; i++)
318 {
319 for (j = 1; j < ny-1; j++)
320 {
321 if (mask[j * nx + i] <= 1) continue;
322 if isnan(jz[j * nx + i]) continue;
323 if isnan(bz[j * nx + i]) continue;
324 if (bz[j * nx + i] == 0.0) continue;
325 sum += (jz[j * nx + i] / bz[j * nx + i])*2.7 ; /* the units for (jz) Gauss/pix; the units for bz are Gauss */
326 //printf("sum=%f\n",sum);
327 //printf("jz[j * nx + i]=%f\n",jz[j * nx + i]);
328 //printf("bz[j * nx + i]=%f\n",bz[j * nx + i]);
329 //printf("(jz[j * nx + i] / bz[j * nx + i])*2.7=%f\n",(jz[j * nx + i] / bz[j * nx + i])*2.7);
330 count_mask++;
331 //printf("count_mask=%d\n",count_mask);
332 }
333 }
334
335 printf("count_mask=%d\n",count_mask);
336 printf("sum=%f\n",sum);
337 mbobra 1.1 *mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */
338 return 0;
339 }
340
341 /*===========================================*/
342
343 /* Example function 10: Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */
344
345 // The current helicity is defined as Bz*Jz and the units are G^2 / m
346 // The units of Jz are in G/pix; the units of Bz are in G.
347 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF) = G^2 / m.
348 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
349 // In that case, we would have the following:
350 // (Gauss/pix)(1/0.5)(1/722500), or
351 // (Gauss/pix)*0.000002768166
352
353 int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr, int *mask)
354 {
355
356
357 int nx = dims[0], ny = dims[1];
358 mbobra 1.1 int i, j, count_mask=0;
359
360 if (nx <= 0 || ny <= 0) return 1;
361
362 *mean_ih_ptr = 0.0;
363 float sum=0.0, sum2=0.0;
364
365 for (j = 0; j < ny; j++)
366 {
367 for (i = 0; i < nx; i++)
368 {
369 if (mask[j * nx + i] <= 1) continue;
370 if isnan(jz[j * nx + i]) continue;
371 if isnan(bz[j * nx + i]) continue;
372 if (bz[j * nx + i] == 0.0) continue;
373 if (jz[j * nx + i] == 0.0) continue;
374 sum += (jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166);
375 sum2 += fabs(jz[j * nx + i]*(bz[j * nx + i]))*(0.000002768166);
376 count_mask++;
377 printf("(jz[j * nx + i])=%f\n",(jz[j * nx + i]));
378 printf("(bz[j * nx + i])=%f\n",(bz[j * nx + i]));
379 mbobra 1.1 printf("(jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166)=%f\n",(jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166));
380 printf("sum=%f\n",sum);
381 printf("sum2=%f\n",sum2);
382 printf("count_mask=%d\n",count_mask);
383 }
384 }
385
386 printf("sum/count_mask=%f\n",sum/count_mask);
387 *mean_ih_ptr = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */
388 *total_us_ih_ptr = sum2; /* Units are G^2 / m */
389 *total_abs_ih_ptr= fabs(sum); /* Units are G^2 / m */
390
391 return 0;
392 }
393
394 /*===========================================*/
395
396 /* Example function 11: Sum of Absolute Value per polarity */
397
398 // The Sum of the Absolute Value per polarity is defined as the following:
399 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
400 mbobra 1.1 // The units of jz are in G/pix. In this case, we would have the following:
401 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(722500)(722500), or
402 // (Gauss/pix)(1.81*10^10)
403 // (Gauss/pix)(18100000000.)
404
405 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr, int *mask)
406 {
407
408
409 int nx = dims[0], ny = dims[1];
410 int i, j, count_mask=0;
411
412 if (nx <= 0 || ny <= 0) return 1;
413
414 *totaljzptr=0.0;
415 float sum1=0.0, sum2=0.0;
416
417 for (i = 0; i < nx; i++)
418 {
419 for (j = 0; j < ny; j++)
420 {
421 mbobra 1.1 if (mask[j * nx + i] <= 1) continue;
422 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i]*18100000000.);
423 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i]*18100000000.);
424 }
425 }
426
427 *totaljzptr = fabs(sum1)+fabs(sum2); /* Units are A */
428 return 0;
429 }
430
431 /*===========================================*/
432
433 /* Example function 12: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
434 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
435 // automatically yields erg per cubic centimeter for an input B in Gauss.
436 //
437 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
438 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
439 // erg/cm^3(CDELT1)(RSUN_REF/RSUN_OBS)(100.)
440 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
441 // = erg/cm^3(1.30501e15)
442 mbobra 1.1 // = erg/cm(1/pix^2)
443
444 int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, float *meanpotptr, float *totpotptr, int *mask)
445 {
446 int nx = dims[0], ny = dims[1];
447 int i, j, count_mask=0;
448
449 if (nx <= 0 || ny <= 0) return 1;
450
451 *totpotptr=0.0;
452 *meanpotptr=0.0;
453 float sum=0.0;
454
455 for (i = 0; i < nx; i++)
456 {
457 for (j = 0; j < ny; j++)
458 {
459 if (mask[j * nx + i] <= 1) continue;
460 sum += (( ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) - ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i])) )/8.*PI);
461 count_mask++;
462 }
463 mbobra 1.1 }
464
465 *meanpotptr = (sum)/(count_mask); /* Units are ergs per cubic centimeter */
466 *totpotptr = sum*(1.30501e15)*(count_mask); /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2, units of count_mask are pix^2; therefore, units of totpotptr are ergs per centimeter */
467 return 0;
468 }
469
470 /*===========================================*/
471 /* Example function 13: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
472
473 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, float *meanshear_angleptr, float *area_w_shear_gt_45ptr, float *meanshear_anglehptr, float *area_w_shear_gt_45hptr, int *mask)
474 {
475 int nx = dims[0], ny = dims[1];
476 int i, j, count_mask=0;
477
478 if (nx <= 0 || ny <= 0) return 1;
479
480 *area_w_shear_gt_45ptr=0.0;
481 *meanshear_angleptr=0.0;
482 float dotproduct, magnitude_potential, magnitude_vector, shear_angle=0.0, sum = 0.0, count=0.0;
483 float dotproducth, magnitude_potentialh, magnitude_vectorh, shear_angleh=0.0, sum1 = 0.0, counth = 0.0;
484 mbobra 1.1
485 for (i = 0; i < nx; i++)
486 {
487 for (j = 0; j < ny; j++)
488 {
489 if (mask[j * nx + i] <= 1) continue;
490 /* For mean 3D shear angle, area with shear greater than 45*/
491 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]);
492 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]));
493 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]) );
494 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
495 sum += shear_angle ;
496 if (shear_angle > 45) count++;
497
498 /* For mean horizontal shear angle, area with horizontal shear angle greater than 45 */
499 dotproducth = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]);
500 magnitude_potentialh = sqrt((bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]));
501 magnitude_vectorh = sqrt( (bx[j * nx + i]*bx[j * nx + i]) + (by[j * nx + i]*by[j * nx + i]) );
502 shear_angleh = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
503 sum1 += shear_angleh ;
504 if (shear_angleh > 45) counth++;
505 mbobra 1.1 count_mask++;
506 }
507 }
508
509 /* For mean 3D shear angle, area with shear greater than 45*/
510 *meanshear_angleptr = (sum)/(count_mask); /* Units are degrees */
511 *area_w_shear_gt_45ptr = (count/(count_mask))*(100.); /* The area here is a fractional area -- the % of the total area */
512
513 // /* For mean horizontal shear angle, area with horizontal shear angle greater than 45 */
514 // *meanshear_anglehptr = (sum1)/(count_mask); /* Units are degrees */
515 // *area_w_shear_gt_45hptr = (counth/(count_mask))*(100.); /* The area here is a fractional area -- the % of the total area */
516 return 0;
517 }
518
519 /*==================KEIJI'S CODE =========================*/
520
521 // #include <omp.h>
522 #include <math.h>
523
524 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
525 {
526 mbobra 1.1 /* local workings */
527 int inx, iny, i, j, n;
528 /* local array */
529 float *pfpot, *rdist;
530 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
531 rdist=(float *)malloc(sizeof(float) *nnx*nny);
532 /* make nan */
533 unsigned long long llnan = 0x7ff0000000000000;
534 // float NAN = (float)(llnan);
535
536 // #pragma omp parallel for private (inx)
537 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
538 // #pragma omp parallel for private (inx)
539 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
540 // #pragma omp parallel for private (inx)
541 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
542 // #pragma omp parallel for private (inx)
543 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
544
545 float dz;
546 dz = params_get_float(&cmdparams, "dzvalue");
547 mbobra 1.1
548 // float dz = 0.0001;
549
550 // #pragma omp parallel for private (inx)
551 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
552 {
553 float rdd, rdd1, rdd2;
554 float r;
555 rdd1 = (float)(inx);
556 rdd2 = (float)(iny);
557 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
558 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
559 }}
560
561 int iwindow;
562 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
563 float rwindow;
564 rwindow = (float)(iwindow);
565 rwindow = rwindow * rwindow + 0.01;
566 // #pragma omp parallel for private(inx)
567 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
568 mbobra 1.1 {
569 float val0 = bz[nnx*iny + inx];
570 if (isnan(val0))
571 {
572 pfpot[nnx*iny + inx] = NAN;
573 }
574 else
575 {
576 float sum;
577 sum = 0.0;
578 int j2, i2;
579 for (j2=0;j2<nny;j2++){for (i2=0;i2<nnx;i2++)
580 {
581 float val1 = bz[nnx*j2 + i2];
582 float rr, r1, r2;
583 r1 = (float)(i2 - inx);
584 r2 = (float)(j2 - iny);
585 rr = r1*r1 + r2*r2;
586 if ((!isnan(val1)) && (rr < rwindow))
587 {
588 int di, dj;
589 mbobra 1.1 di = abs(i2 - inx);
590 dj = abs(j2 - iny);
591 sum = sum + val1 * rdist[nnx * dj + di] * dz;
592 }
593 } }
594 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
595 }
596 } } // end of OpenMP parallelism
597
598 // #pragma omp parallel for private(inx)
599 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
600 {
601 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) / 2.0;
602 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) / 2.0;
603 } }
604 free(pfpot);
605 } // end of void func. greenpot
606
607
608 /*===========END OF KEIJI'S CODE =========================*/
609 /* ---------------- end of this file ----------------*/
|
Return to
sharp_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps