1 xudong 1.1 /*===========================================
2
|
3 mbobra 1.5 The following 14 functions calculate the following spaceweather indices:
|
4 xudong 1.1
5 USFLUX Total unsigned flux in Maxwells
6 MEANGAM Mean inclination angle, gamma, in degrees
7 MEANGBT Mean value of the total field gradient, in Gauss/Mm
8 MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
9 MEANGBH Mean value of the horizontal field gradient, in Gauss/Mm
10 MEANJZD Mean vertical current density, in mA/m2
11 TOTUSJZ Total unsigned vertical current, in Amperes
12 MEANALP Mean twist parameter, alpha, in 1/Mm
13 MEANJZH Mean current helicity in G2/m
14 TOTUSJH Total unsigned current helicity in G2/m
15 ABSNJZH Absolute value of the net current helicity in G2/m
16 SAVNCPP Sum of the Absolute Value of the Net Currents Per Polarity in Amperes
17 MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
18 TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
19 MEANSHR Mean shear angle (measured using Btotal) in degrees
20
|
21 mbobra 1.3 The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
22 pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
|
23 mbobra 1.20 coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data
24 prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI
25 contain nearest-neighbor bitmaps).
|
26 mbobra 1.3
27 In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
|
28 mbobra 1.20 and the pixels that equal 33 or 34 in bitmap. Here are the definitions of the pixel values:
|
29 mbobra 1.3
30 For conf_disambig:
31 50 : not all solutions agree (weak field method applied)
32 60 : not all solutions agree (weak field + annealed)
33 90 : all solutions agree (strong field + annealed)
34 0 : not disambiguated
35
36 For bitmap:
37 1 : weak field outside smooth bounding curve
38 2 : strong field outside smooth bounding curve
39 33 : weak field inside smooth bounding curve
40 34 : strong field inside smooth bounding curve
|
41 xudong 1.1
|
42 mbobra 1.20 Written by Monica Bobra 15 August 2012
|
43 xudong 1.1 Potential Field code (appended) written by Keiji Hayashi
|
44 mbobra 1.20 Error analysis modification 21 October 2013
|
45 xudong 1.1
46 ===========================================*/
47 #include <math.h>
|
48 mbobra 1.9 #include <mkl.h>
|
49 xudong 1.1
|
50 mbobra 1.9 #define PI (M_PI)
|
51 xudong 1.1 #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
52
53 /*===========================================*/
54
55 /* Example function 1: Compute total unsigned flux in units of G/cm^2 */
56
57 // To compute the unsigned flux, we simply calculate
58 // flux = surface integral [(vector Bz) dot (normal vector)],
59 // = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
60 // However, since the field is radial, we will assume cos theta = 1.
61 // Therefore the pixels only need to be corrected for the projection.
62
63 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
64 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
65 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
|
66 mbobra 1.20 // =Gauss*cm^2
|
67 xudong 1.1
|
68 mbobra 1.9 int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
69 float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
70 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
71 xudong 1.1
72 {
73
|
74 mbobra 1.14 int nx = dims[0];
75 int ny = dims[1];
|
76 mbobra 1.15 int i = 0;
77 int j = 0;
78 int count_mask = 0;
79 double sum = 0.0;
80 double err = 0.0;
|
81 mbobra 1.14 *absFlux = 0.0;
82 *mean_vf_ptr = 0.0;
83
84
|
85 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
86
|
87 mbobra 1.5 for (i = 0; i < nx; i++)
|
88 xudong 1.1 {
|
89 mbobra 1.5 for (j = 0; j < ny; j++)
|
90 xudong 1.1 {
|
91 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
92 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
93 xudong 1.1 sum += (fabs(bz[j * nx + i]));
|
94 mbobra 1.9 err += bz_err[j * nx + i]*bz_err[j * nx + i];
|
95 xudong 1.1 count_mask++;
96 }
97 }
98
|
99 mbobra 1.9 *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
100 *mean_vf_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); // error in the unsigned flux
101 *count_mask_ptr = count_mask;
|
102 xudong 1.1 return 0;
103 }
104
105 /*===========================================*/
|
106 mbobra 1.9 /* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
|
107 xudong 1.1 // Native units of Bh are Gauss
108
|
109 mbobra 1.9 int computeBh(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
|
110 mbobra 1.3 float *mean_hf_ptr, int *mask, int *bitmask)
|
111 xudong 1.1
112 {
113
|
114 mbobra 1.14 int nx = dims[0];
115 int ny = dims[1];
|
116 mbobra 1.15 int i = 0;
117 int j = 0;
118 int count_mask = 0;
119 double sum = 0.0;
|
120 mbobra 1.9 *mean_hf_ptr = 0.0;
|
121 xudong 1.1
122 if (nx <= 0 || ny <= 0) return 1;
123
|
124 mbobra 1.5 for (i = 0; i < nx; i++)
|
125 xudong 1.1 {
|
126 mbobra 1.5 for (j = 0; j < ny; j++)
|
127 mbobra 1.20 {
128 if isnan(bx[j * nx + i])
129 {
130 bh[j * nx + i] = NAN;
131 bh_err[j * nx + i] = NAN;
132 continue;
133 }
134 if isnan(by[j * nx + i])
135 {
136 bh[j * nx + i] = NAN;
137 bh_err[j * nx + i] = NAN;
138 continue;
139 }
140 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
141 sum += bh[j * nx + i];
142 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];
143 count_mask++;
|
144 xudong 1.1 }
145 }
146
147 *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
|
148 mbobra 1.9
|
149 xudong 1.1 return 0;
150 }
151
152 /*===========================================*/
153 /* Example function 3: Calculate Gamma in units of degrees */
154 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
|
155 mbobra 1.20 //
156 // Error analysis calculations are done in radians (since derivatives are only true in units of radians),
157 // and multiplied by (180./PI) at the end for consistency in units.
|
158 xudong 1.1
|
159 mbobra 1.9 int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
160 float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)
|
161 xudong 1.1 {
|
162 mbobra 1.14 int nx = dims[0];
163 int ny = dims[1];
|
164 mbobra 1.15 int i = 0;
165 int j = 0;
166 int count_mask = 0;
167 double sum = 0.0;
168 double err = 0.0;
169 *mean_gamma_ptr = 0.0;
|
170 xudong 1.1
171 if (nx <= 0 || ny <= 0) return 1;
172
173 for (i = 0; i < nx; i++)
174 {
175 for (j = 0; j < ny; j++)
176 {
177 if (bh[j * nx + i] > 100)
178 {
|
179 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
180 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
181 mbobra 1.9 if isnan(bz_err[j * nx + i]) continue;
|
182 mbobra 1.20 if isnan(bh_err[j * nx + i]) continue;
183 if isnan(bh[j * nx + i]) continue;
|
184 mbobra 1.9 if (bz[j * nx + i] == 0) continue;
|
185 mbobra 1.19 sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))*(180./PI);
|
186 mbobra 1.20 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])))) *
187 ( ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) +
188 ((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])) );
|
189 mbobra 1.9 count_mask++;
|
190 xudong 1.1 }
191 }
192 }
193
194 *mean_gamma_ptr = sum/count_mask;
|
195 mbobra 1.20 *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);
|
196 mbobra 1.16 //printf("MEANGAM=%f\n",*mean_gamma_ptr);
197 //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
|
198 xudong 1.1 return 0;
199 }
200
201 /*===========================================*/
202 /* Example function 4: Calculate B_Total*/
203 // Native units of B_Total are in gauss
204
|
205 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)
|
206 xudong 1.1 {
207
|
208 mbobra 1.14 int nx = dims[0];
209 int ny = dims[1];
|
210 mbobra 1.15 int i = 0;
211 int j = 0;
212 int count_mask = 0;
|
213 xudong 1.1
214 if (nx <= 0 || ny <= 0) return 1;
215
216 for (i = 0; i < nx; i++)
217 {
218 for (j = 0; j < ny; j++)
219 {
|
220 mbobra 1.20 if isnan(bx[j * nx + i])
221 {
222 bt[j * nx + i] = NAN;
223 bt_err[j * nx + i] = NAN;
224 continue;
225 }
226 if isnan(by[j * nx + i])
227 {
228 bt[j * nx + i] = NAN;
229 bt_err[j * nx + i] = NAN;
230 continue;
231 }
232 if isnan(bz[j * nx + i])
233 {
234 bt[j * nx + i] = NAN;
235 bt_err[j * nx + i] = NAN;
236 continue;
237 }
238 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]);
239 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];
|
240 xudong 1.1 }
241 }
242 return 0;
243 }
244
245 /*===========================================*/
246 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
247
|
248 mbobra 1.9 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)
|
249 xudong 1.1 {
250
|
251 mbobra 1.14 int nx = dims[0];
252 int ny = dims[1];
|
253 mbobra 1.15 int i = 0;
254 int j = 0;
255 int count_mask = 0;
256 double sum = 0.0;
257 double err = 0.0;
|
258 mbobra 1.14 *mean_derivative_btotal_ptr = 0.0;
|
259 xudong 1.1
260 if (nx <= 0 || ny <= 0) return 1;
261
262 /* brute force method of calculating the derivative (no consideration for edges) */
|
263 mbobra 1.10 for (i = 1; i <= nx-2; i++)
|
264 xudong 1.1 {
265 for (j = 0; j <= ny-1; j++)
266 {
|
267 mbobra 1.10 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
|
268 xudong 1.1 }
269 }
270
271 /* brute force method of calculating the derivative (no consideration for edges) */
272 for (i = 0; i <= nx-1; i++)
273 {
|
274 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
275 xudong 1.1 {
|
276 mbobra 1.10 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
|
277 xudong 1.1 }
278 }
279
280
|
281 mbobra 1.10 /* consider the edges */
|
282 xudong 1.1 i=0;
283 for (j = 0; j <= ny-1; j++)
284 {
|
285 mbobra 1.10 derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
|
286 xudong 1.1 }
287
288 i=nx-1;
289 for (j = 0; j <= ny-1; j++)
290 {
|
291 mbobra 1.10 derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
|
292 mbobra 1.9 }
293
|
294 xudong 1.1 j=0;
295 for (i = 0; i <= nx-1; i++)
296 {
|
297 mbobra 1.10 dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5;
|
298 xudong 1.1 }
299
300 j=ny-1;
301 for (i = 0; i <= nx-1; i++)
302 {
|
303 mbobra 1.10 dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
|
304 xudong 1.1 }
305
306
|
307 mbobra 1.20 for (i = 1; i <= nx-2; i++)
|
308 xudong 1.1 {
|
309 mbobra 1.20 for (j = 1; j <= ny-2; j++)
|
310 xudong 1.1 {
|
311 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
312 mbobra 1.21 if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue;
313 if isnan(bt[j * nx + i]) continue;
314 if isnan(bt[(j+1) * nx + i]) continue;
315 if isnan(bt[(j-1) * nx + i]) continue;
316 if isnan(bt[j * nx + i-1]) continue;
317 if isnan(bt[j * nx + i+1]) continue;
318 if isnan(bt_err[j * nx + i]) continue;
|
319 mbobra 1.5 if isnan(derx_bt[j * nx + i]) continue;
320 if isnan(dery_bt[j * nx + i]) continue;
|
321 xudong 1.1 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 */
|
322 mbobra 1.20 err += (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] ))+
323 (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] )) ;
|
324 xudong 1.1 count_mask++;
325 }
326 }
327
|
328 mbobra 1.20 *mean_derivative_btotal_ptr = (sum)/(count_mask);
329 *mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask);
|
330 mbobra 1.16 //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
331 //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
|
332 mbobra 1.20
|
333 xudong 1.1 return 0;
334 }
335
336
337 /*===========================================*/
338 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
339
|
340 mbobra 1.9 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)
|
341 xudong 1.1 {
342
|
343 mbobra 1.14 int nx = dims[0];
344 int ny = dims[1];
|
345 mbobra 1.15 int i = 0;
346 int j = 0;
347 int count_mask = 0;
348 double sum= 0.0;
349 double err =0.0;
|
350 mbobra 1.14 *mean_derivative_bh_ptr = 0.0;
|
351 xudong 1.1
352 if (nx <= 0 || ny <= 0) return 1;
353
354 /* brute force method of calculating the derivative (no consideration for edges) */
|
355 mbobra 1.10 for (i = 1; i <= nx-2; i++)
|
356 xudong 1.1 {
357 for (j = 0; j <= ny-1; j++)
358 {
|
359 mbobra 1.10 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
|
360 xudong 1.1 }
361 }
362
363 /* brute force method of calculating the derivative (no consideration for edges) */
364 for (i = 0; i <= nx-1; i++)
365 {
|
366 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
367 xudong 1.1 {
|
368 mbobra 1.10 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
|
369 xudong 1.1 }
370 }
371
372
|
373 mbobra 1.10 /* consider the edges */
|
374 xudong 1.1 i=0;
375 for (j = 0; j <= ny-1; j++)
376 {
|
377 mbobra 1.10 derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
|
378 xudong 1.1 }
379
380 i=nx-1;
381 for (j = 0; j <= ny-1; j++)
382 {
|
383 mbobra 1.10 derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
|
384 mbobra 1.9 }
385
|
386 xudong 1.1 j=0;
387 for (i = 0; i <= nx-1; i++)
388 {
|
389 mbobra 1.10 dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5;
|
390 xudong 1.1 }
391
392 j=ny-1;
393 for (i = 0; i <= nx-1; i++)
394 {
|
395 mbobra 1.10 dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
|
396 xudong 1.1 }
397
398
399 for (i = 0; i <= nx-1; i++)
400 {
401 for (j = 0; j <= ny-1; j++)
402 {
|
403 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
404 mbobra 1.21 if ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue;
405 if isnan(bh[j * nx + i]) continue;
406 if isnan(bh[(j+1) * nx + i]) continue;
407 if isnan(bh[(j-1) * nx + i]) continue;
408 if isnan(bh[j * nx + i-1]) continue;
409 if isnan(bh[j * nx + i+1]) continue;
410 if isnan(bh_err[j * nx + i]) continue;
|
411 mbobra 1.5 if isnan(derx_bh[j * nx + i]) continue;
412 if isnan(dery_bh[j * nx + i]) continue;
|
413 xudong 1.1 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 */
|
414 mbobra 1.20 err += (((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])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] ))+
415 (((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)])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] )) ;
|
416 xudong 1.1 count_mask++;
417 }
418 }
419
|
420 mbobra 1.9 *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
421 *mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
|
422 mbobra 1.16 //printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
423 //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
|
424 mbobra 1.9
|
425 xudong 1.1 return 0;
426 }
427
428 /*===========================================*/
429 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
430
|
431 mbobra 1.9 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)
|
432 xudong 1.1 {
433
|
434 mbobra 1.14 int nx = dims[0];
435 int ny = dims[1];
|
436 mbobra 1.15 int i = 0;
437 int j = 0;
438 int count_mask = 0;
439 double sum = 0.0;
440 double err = 0.0;
|
441 mbobra 1.14 *mean_derivative_bz_ptr = 0.0;
|
442 xudong 1.1
443 if (nx <= 0 || ny <= 0) return 1;
444
445 /* brute force method of calculating the derivative (no consideration for edges) */
|
446 mbobra 1.10 for (i = 1; i <= nx-2; i++)
|
447 xudong 1.1 {
448 for (j = 0; j <= ny-1; j++)
449 {
|
450 mbobra 1.10 if isnan(bz[j * nx + i]) continue;
451 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
|
452 xudong 1.1 }
453 }
454
455 /* brute force method of calculating the derivative (no consideration for edges) */
456 for (i = 0; i <= nx-1; i++)
457 {
|
458 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
459 xudong 1.1 {
|
460 mbobra 1.10 if isnan(bz[j * nx + i]) continue;
461 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
|
462 xudong 1.1 }
463 }
464
465
|
466 mbobra 1.10 /* consider the edges */
|
467 xudong 1.1 i=0;
468 for (j = 0; j <= ny-1; j++)
469 {
|
470 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
471 mbobra 1.10 derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
|
472 xudong 1.1 }
473
474 i=nx-1;
475 for (j = 0; j <= ny-1; j++)
476 {
|
477 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
478 mbobra 1.10 derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
|
479 xudong 1.1 }
480
481 j=0;
482 for (i = 0; i <= nx-1; i++)
483 {
|
484 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
485 mbobra 1.10 dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
|
486 xudong 1.1 }
487
488 j=ny-1;
489 for (i = 0; i <= nx-1; i++)
490 {
|
491 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
492 mbobra 1.10 dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
|
493 xudong 1.1 }
494
495
496 for (i = 0; i <= nx-1; i++)
497 {
498 for (j = 0; j <= ny-1; j++)
499 {
|
500 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
501 mbobra 1.21 if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
502 if isnan(bz[j * nx + i]) continue;
503 if isnan(bz[(j+1) * nx + i]) continue;
504 if isnan(bz[(j-1) * nx + i]) continue;
505 if isnan(bz[j * nx + i-1]) continue;
506 if isnan(bz[j * nx + i+1]) continue;
507 if isnan(bz_err[j * nx + i]) continue;
|
508 mbobra 1.4 if isnan(derx_bz[j * nx + i]) continue;
509 if isnan(dery_bz[j * nx + i]) continue;
|
510 xudong 1.1 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 */
|
511 mbobra 1.21 err += (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i])) /
512 (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] )) +
513 (((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)])) /
514 (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] )) ;
|
515 xudong 1.1 count_mask++;
516 }
517 }
518
519 *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
|
520 mbobra 1.9 *mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
|
521 mbobra 1.16 //printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
522 //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
|
523 mbobra 1.9
|
524 xudong 1.1 return 0;
525 }
526
527 /*===========================================*/
528 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */
529
530 // In discretized space like data pixels,
531 // the current (or curl of B) is calculated as
532 // the integration of the field Bx and By along
533 // the circumference of the data pixel divided by the area of the pixel.
534 // One form of differencing (a word for the differential operator
|
535 mbobra 1.10 // in the discretized space) of the curl is expressed as
|
536 xudong 1.1 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2
537 // +dy * (By(i+1,j)+By(i,j)) / 2
538 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
539 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
540 //
541 //
542 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
543 // one must perform the following unit conversions:
|
544 mbobra 1.8 // (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
545 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or
546 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
|
547 xudong 1.1 // where a Tesla is represented as a Newton/Ampere*meter.
|
548 mbobra 1.4 //
|
549 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
550 // In that case, we would have the following:
551 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
552 // jz * (35.0)
553 //
554 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
|
555 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)
556 // = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
|
557 xudong 1.1
558
|
559 mbobra 1.9 // Comment out random number generator, which can only run on solar3
560 //int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
561 // int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
562 // float *noiseby, float *noisebz)
563
564 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
565 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
566
|
567 xudong 1.1
|
568 mbobra 1.10 {
|
569 mbobra 1.14 int nx = dims[0];
570 int ny = dims[1];
|
571 mbobra 1.15 int i = 0;
572 int j = 0;
573 int count_mask = 0;
|
574 mbobra 1.10
575 if (nx <= 0 || ny <= 0) return 1;
|
576 xudong 1.1
|
577 mbobra 1.9 /* Calculate the derivative*/
|
578 xudong 1.1 /* brute force method of calculating the derivative (no consideration for edges) */
|
579 mbobra 1.10
580 for (i = 1; i <= nx-2; i++)
|
581 xudong 1.1 {
582 for (j = 0; j <= ny-1; j++)
583 {
|
584 mbobra 1.12 if isnan(by[j * nx + i]) continue;
|
585 mbobra 1.10 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
|
586 xudong 1.1 }
587 }
588
589 for (i = 0; i <= nx-1; i++)
590 {
|
591 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
592 xudong 1.1 {
|
593 mbobra 1.12 if isnan(bx[j * nx + i]) continue;
|
594 mbobra 1.10 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
|
595 xudong 1.1 }
596 }
597
|
598 mbobra 1.10 // consider the edges
|
599 xudong 1.1 i=0;
600 for (j = 0; j <= ny-1; j++)
601 {
|
602 mbobra 1.4 if isnan(by[j * nx + i]) continue;
|
603 mbobra 1.10 derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
|
604 xudong 1.1 }
605
606 i=nx-1;
607 for (j = 0; j <= ny-1; j++)
608 {
|
609 mbobra 1.4 if isnan(by[j * nx + i]) continue;
|
610 mbobra 1.10 derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
611 }
|
612 mbobra 1.9
|
613 xudong 1.1 j=0;
614 for (i = 0; i <= nx-1; i++)
615 {
|
616 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
|
617 mbobra 1.10 dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
|
618 xudong 1.1 }
619
620 j=ny-1;
|
621 mbobra 1.11 for (i = 0; i <= nx-1; i++)
|
622 xudong 1.1 {
|
623 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
|
624 mbobra 1.10 dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
|
625 mbobra 1.9 }
626
|
627 mbobra 1.17 for (i = 1; i <= nx-2; i++)
|
628 xudong 1.1 {
|
629 mbobra 1.17 for (j = 1; j <= ny-2; j++)
|
630 xudong 1.1 {
|
631 mbobra 1.17 // calculate jz at all points
632
|
633 mbobra 1.15 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix
634 jz_err[j * nx + i] = 0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +
|
635 mbobra 1.19 (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;
|
636 mbobra 1.15 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
|
637 mbobra 1.10 count_mask++;
|
638 mbobra 1.17
|
639 mbobra 1.5 }
|
640 mbobra 1.10 }
641 return 0;
642 }
|
643 mbobra 1.5
644 /*===========================================*/
645
|
646 mbobra 1.11 /* Example function 9: Compute quantities on Jz array */
647 // 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 float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
651 float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
|
652 mbobra 1.5
653 {
654
|
655 mbobra 1.14 int nx = dims[0];
656 int ny = dims[1];
|
657 mbobra 1.15 int i = 0;
658 int j = 0;
659 int count_mask = 0;
660 double curl = 0.0;
661 double us_i = 0.0;
662 double err = 0.0;
|
663 mbobra 1.5
664 if (nx <= 0 || ny <= 0) return 1;
665
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 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
672 mbobra 1.4 if isnan(derx[j * nx + i]) continue;
673 if isnan(dery[j * nx + i]) continue;
|
674 mbobra 1.9 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 xudong 1.1 count_mask++;
|
679 mbobra 1.9 }
|
680 xudong 1.1 }
681
|
682 mbobra 1.15 /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
|
683 xudong 1.1 *mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD */
|
684 mbobra 1.20 *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 mbobra 1.9
|
686 mbobra 1.4 *us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ */
|
687 mbobra 1.20 *us_i_err_ptr = (sqrt(err))*((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ
|
688 mbobra 1.9
|
689 mbobra 1.16 //printf("MEANJZD=%f\n",*mean_jz_ptr);
690 //printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
|
691 mbobra 1.9
|
692 mbobra 1.16 //printf("TOTUSJZ=%g\n",*us_i_ptr);
693 //printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
|
694 mbobra 1.9
|
695 xudong 1.1 return 0;
696
697 }
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 mbobra 1.5
|
706 mbobra 1.19 // numerator = sum of all Jz*Bz
707 // denominator = sum of Bz*Bz
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 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
714 // = (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 mbobra 1.19 int nx = dims[0];
721 int ny = dims[1];
722 int i = 0;
723 int j = 0;
724 double alpha_total = 0.0;
725 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 xudong 1.1
730 if (nx <= 0 || ny <= 0) return 1;
|
731 mbobra 1.19 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[j*nx+i]*bz[j*nx+i];
741 B += bz[j*nx+i]*bz[j*nx+i];
742 }
743 }
|
744 xudong 1.1
745 for (i = 1; i < nx-1; i++)
746 {
747 for (j = 1; j < ny-1; j++)
748 {
|
749 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
750 mbobra 1.19 if isnan(jz[j * nx + i]) continue;
751 if isnan(bz[j * nx + i]) continue;
752 if (jz[j * nx + i] == 0.0) continue;
753 if (bz[j * nx + i] == 0.0) continue;
754 total += bz[j*nx+i]*bz[j*nx+i]*jz_err[j*nx+i]*jz_err[j*nx+i] + (jz[j*nx+i]-2*bz[j*nx+i]*A/B)*(jz[j*nx+i]-2*bz[j*nx+i]*A/B)*bz_err[j*nx+i]*bz_err[j*nx+i];
755 }
|
756 xudong 1.1 }
|
757 mbobra 1.14
|
758 mbobra 1.19 /* 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 mbobra 1.9
|
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 mbobra 1.9 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter)
|
772 xudong 1.1 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
|
773 mbobra 1.9 // = G^2 / m.
|
774 xudong 1.1
|
775 mbobra 1.9 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 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
|
781 mbobra 1.20 int nx = dims[0];
782 int ny = dims[1];
783 int i = 0;
784 int j = 0;
|
785 mbobra 1.15 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
|
792 mbobra 1.5 for (i = 0; i < nx; i++)
|
793 xudong 1.1 {
|
794 mbobra 1.5 for (j = 0; j < ny; j++)
|
795 xudong 1.1 {
|
796 mbobra 1.9 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
797 mbobra 1.20 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 mbobra 1.9 sum += (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
|
804 mbobra 1.14 sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
|
805 mbobra 1.20 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 mbobra 1.9 count_mask++;
|
807 xudong 1.1 }
808 }
809
|
810 mbobra 1.9 *mean_ih_ptr = sum/count_mask ; /* Units are G^2 / m ; keyword is MEANJZH */
811 *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
|
814 mbobra 1.20 *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 mbobra 1.9
|
818 mbobra 1.21 //printf("MEANJZH=%f\n",*mean_ih_ptr);
819 //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
|
820 mbobra 1.9
|
821 mbobra 1.21 //printf("TOTUSJH=%f\n",*total_us_ih_ptr);
822 //printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
|
823 mbobra 1.9
|
824 mbobra 1.21 //printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
825 //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
|
826 xudong 1.1
827 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 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
835 // 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 mbobra 1.9 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 {
|
845 mbobra 1.14 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.15 double sum1=0.0;
851 double sum2=0.0;
852 double err=0.0;
|
853 mbobra 1.14 *totaljzptr=0.0;
|
854 xudong 1.1
855 if (nx <= 0 || ny <= 0) return 1;
856
857 for (i = 0; i < nx; i++)
858 {
859 for (j = 0; j < ny; j++)
860 {
|
861 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
862 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
863 mbobra 1.18 if isnan(jz[j * nx + i]) continue;
|
864 mbobra 1.9 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 xudong 1.1 }
869 }
870
|
871 mbobra 1.9 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */
872 *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
|
873 mbobra 1.16 //printf("SAVNCPP=%g\n",*totaljzptr);
874 //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
|
875 mbobra 1.9
|
876 xudong 1.1 return 0;
877 }
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 mbobra 1.11 // automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8*PI can come out of the integral; thus,
883 // the integral is over B^2 dV and the 8*PI 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 mbobra 1.9 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 mbobra 1.14 int nx = dims[0];
898 int ny = dims[1];
|
899 mbobra 1.15 int i = 0;
900 int j = 0;
901 int count_mask = 0;
902 double sum = 0.0;
903 double sum1 = 0.0;
904 double err = 0.0;
905 *totpotptr = 0.0;
906 *meanpotptr = 0.0;
|
907 mbobra 1.14
908 if (nx <= 0 || ny <= 0) return 1;
|
909 xudong 1.1
910 for (i = 0; i < nx; i++)
911 {
912 for (j = 0; j < ny; j++)
913 {
|
914 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
915 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
916 if isnan(by[j * nx + i]) continue;
|
917 mbobra 1.20 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])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.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 xudong 1.1 count_mask++;
922 }
923 }
924
|
925 mbobra 1.20 /* Units of meanpotptr are ergs per centimeter */
926 *meanpotptr = (sum1) / (count_mask*8.*PI) ; /* Units are ergs per cubic centimeter */
927 *meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
|
928 mbobra 1.9
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 mbobra 1.11 *totpotptr = (sum)/(8.*PI);
931 *totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
|
932 mbobra 1.9
|
933 mbobra 1.16 //printf("MEANPOT=%g\n",*meanpotptr);
934 //printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
|
935 mbobra 1.9
|
936 mbobra 1.16 //printf("TOTPOT=%g\n",*totpotptr);
937 //printf("TOTPOT_err=%g\n",*totpot_err_ptr);
|
938 mbobra 1.9
|
939 xudong 1.1 return 0;
940 }
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.1 {
|
950 mbobra 1.14 int nx = dims[0];
951 int ny = dims[1];
|
952 mbobra 1.15 int i = 0;
953 int j = 0;
|
954 mbobra 1.21 float count_mask = 0;
955 float count = 0;
|
956 mbobra 1.15 double dotproduct = 0.0;
957 double magnitude_potential = 0.0;
958 double magnitude_vector = 0.0;
959 double shear_angle = 0.0;
|
960 mbobra 1.20 double denominator = 0.0;
961 double term1 = 0.0;
962 double term2 = 0.0;
963 double term3 = 0.0;
|
964 mbobra 1.21 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 mbobra 1.15 *area_w_shear_gt_45ptr = 0.0;
970 *meanshear_angleptr = 0.0;
|
971 xudong 1.1
972 if (nx <= 0 || ny <= 0) return 1;
973
974 for (i = 0; i < nx; i++)
975 {
976 for (j = 0; j < ny; j++)
977 {
|
978 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
979 xudong 1.1 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 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
984 if isnan(by[j * nx + i]) continue;
|
985 mbobra 1.21 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 mbobra 1.20 /* For mean 3D shear angle, percentage with shear greater than 45*/
|
990 xudong 1.1 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 mbobra 1.9 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 mbobra 1.21 //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_potential*magnitude_vector))*(180./PI);
998 sumsum += shear_angle;
999 //printf("shear_angle=%f\n",shear_angle);
|
1000 xudong 1.1 count ++;
|
1001 mbobra 1.21
|
1002 xudong 1.1 if (shear_angle > 45) count_mask ++;
|
1003 mbobra 1.20
|
1004 mbobra 1.21 // 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 mbobra 1.20 // denominator is squared
|
1015 mbobra 1.21 denominator = part1*part1*part1*part2*(1.0-((part3*part3)/(part1*part2)));
1016
1017 err = (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
1018 (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
1019 (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) ;
1020
1021 }
1022 }
|
1023 xudong 1.1 /* For mean 3D shear angle, area with shear greater than 45*/
|
1024 mbobra 1.21 *meanshear_angleptr = (sumsum)/(count); /* Units are degrees */
1025 *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);
|
1026 mbobra 1.20
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 mbobra 1.9
|
1030 mbobra 1.22 //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 xudong 1.1
1034 return 0;
1035 }
1036
|
1037 mbobra 1.23 /*===========================================*/
1038 int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1,
1039 float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
1040 float *pmap, int nx1, int ny1)
1041 {
1042
1043 int nx = dims[0];
1044 int ny = dims[1];
1045 int i = 0;
1046 int j = 0;
1047 int index;
1048 double sum = 0.0;
1049 double err = 0.0;
1050 *Rparam = 0.0;
1051 struct fresize_struct fresboxcar, fresgauss;
1052 int scale = round(2.0/cdelt1);
1053 float sigma = 10.0/2.3548;
1054
1055 init_fresize_boxcar(&fresboxcar,1,1);
1056
1057 // setup convolution kernel
1058 mbobra 1.23 //init_fresize_gaussian(&fresgauss,sigma,4,1);
1059 init_fresize_gaussian(&fresgauss,sigma,20,1);
1060
1061 if ((nx || ny) < 40.) return -1;
1062
1063 //float *test = (float *)malloc(nx1*ny1*sizeof(float));
1064
1065 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
1066 for (i = 0; i < nx1; i++)
1067 {
1068 for (j = 0; j < ny1; j++)
1069 {
1070 index = j * nx1 + i;
1071 if (rim[index] > 150)
1072 p1p0[index]=1.0;
1073 else
1074 p1p0[index]=0.0;
1075 if (rim[index] < -150)
1076 p1n0[index]=1.0;
1077 else
1078 p1n0[index]=0.0;
1079 mbobra 1.23 }
1080 }
1081
1082 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1083 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1084
1085 for (i = 0; i < nx1; i++)
1086 {
1087 for (j = 0; j < ny1; j++)
1088 {
1089 index = j * nx1 + i;
1090 if (p1p[index] > 0 && p1n[index] > 0)
1091 p1[index]=1.0;
1092 else
1093 p1[index]=0.0;
1094 }
1095 }
1096
1097 fresize(&fresgauss, p1, pmap, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1098
1099 for (i = 0; i < nx1; i++)
1100 mbobra 1.23 {
1101 for (j = 0; j < ny1; j++)
1102 {
1103 index = j * nx1 + i;
1104 sum += pmap[index]*abs(rim[index]);
1105 }
1106 }
1107
1108 if (sum < 1.0)
1109 *Rparam = 0.0;
1110 else
1111 *Rparam = log10(sum);
1112
1113 printf("Rparam=%f",*Rparam);
1114
1115 free_fresize(&fresboxcar);
1116 free_fresize(&fresgauss);
1117
1118 return 0;
1119 }
1120
|
1121 xudong 1.1
1122 /*==================KEIJI'S CODE =========================*/
1123
1124 // #include <omp.h>
1125 #include <math.h>
1126
1127 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
1128 {
1129 /* local workings */
1130 int inx, iny, i, j, n;
1131 /* local array */
1132 float *pfpot, *rdist;
1133 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
1134 rdist=(float *)malloc(sizeof(float) *nnx*nny);
1135 float *bztmp;
1136 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
1137 /* make nan */
1138 // unsigned long long llnan = 0x7ff0000000000000;
1139 // float NAN = (float)(llnan);
1140
1141 // #pragma omp parallel for private (inx)
1142 xudong 1.1 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
1143 // #pragma omp parallel for private (inx)
1144 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
1145 // #pragma omp parallel for private (inx)
1146 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
1147 // #pragma omp parallel for private (inx)
1148 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
1149 // #pragma omp parallel for private (inx)
1150 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1151 {
1152 float val0 = bz[nnx*iny + inx];
1153 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
1154 }}
1155
1156 // dz is the monopole depth
1157 float dz = 0.001;
1158
1159 // #pragma omp parallel for private (inx)
1160 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1161 {
1162 float rdd, rdd1, rdd2;
1163 xudong 1.1 float r;
1164 rdd1 = (float)(inx);
1165 rdd2 = (float)(iny);
1166 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
1167 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
1168 }}
1169
1170 int iwindow;
1171 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
1172 float rwindow;
1173 rwindow = (float)(iwindow);
1174 rwindow = rwindow * rwindow + 0.01; // must be of square
1175
1176 rwindow = 1.0e2; // limit the window size to be 10.
1177
1178 rwindow = sqrt(rwindow);
1179 iwindow = (int)(rwindow);
1180
1181 // #pragma omp parallel for private(inx)
1182 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
1183 {
1184 xudong 1.1 float val0 = bz[nnx*iny + inx];
1185 if (isnan(val0))
1186 {
1187 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
1188 }
1189 else
1190 {
1191 float sum;
1192 sum = 0.0;
1193 int j2, i2;
1194 int j2s, j2e, i2s, i2e;
1195 j2s = iny - iwindow;
1196 j2e = iny + iwindow;
1197 if (j2s < 0){j2s = 0;}
1198 if (j2e > nny){j2e = nny;}
1199 i2s = inx - iwindow;
1200 i2e = inx + iwindow;
1201 if (i2s < 0){i2s = 0;}
1202 if (i2e > nnx){i2e = nnx;}
1203
1204 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
1205 xudong 1.1 {
1206 float val1 = bztmp[nnx*j2 + i2];
1207 float rr, r1, r2;
1208 // r1 = (float)(i2 - inx);
1209 // r2 = (float)(j2 - iny);
1210 // rr = r1*r1 + r2*r2;
1211 // if (rr < rwindow)
1212 // {
1213 int di, dj;
1214 di = abs(i2 - inx);
1215 dj = abs(j2 - iny);
1216 sum = sum + val1 * rdist[nnx * dj + di] * dz;
1217 // }
1218 } }
1219 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
1220 }
1221 } } // end of OpenMP parallelism
1222
1223 // #pragma omp parallel for private(inx)
1224 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1225 {
1226 xudong 1.1 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
1227 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
1228 } } // end of OpenMP parallelism
1229
1230 free(rdist);
1231 free(pfpot);
1232 free(bztmp);
1233 } // end of void func. greenpot
1234
1235
1236 /*===========END OF KEIJI'S CODE =========================*/
|
1237 mbobra 1.14
1238 char *sw_functions_version() // Returns CVS version of sw_functions.c
1239 {
|
1240 mbobra 1.23 return strdup("$Id: sw_functions.c,v 1.22 2013/11/11 23:21:21 mbobra Exp $");
|
1241 mbobra 1.14 }
1242
|
1243 xudong 1.1 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps