ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/JSOC/proj/sharp/apps/sw_functions.c

Comparing proj/sharp/apps/sw_functions.c (file contents):
Revision 1.1 by xudong, Thu Aug 23 08:38:33 2012 UTC vs.
Revision 1.13 by mbobra, Fri May 31 22:47:15 2013 UTC

#	Line 1 | Line 1
1		/*===========================================
2
3	<	The following 13 functions calculate the following spaceweather indices:
3	>	The following 14 functions calculate the following spaceweather indices:
4
5		USFLUX Total unsigned flux in Maxwells
6		MEANGAM Mean inclination angle, gamma, in degrees
#	Line 18 | Line 18
18		TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
19		MEANSHR Mean shear angle (measured using Btotal) in degrees
20
21	<	The indices are calculated on the pixels in which the disambig bitmap equals 5 or 7:
22	<	5: pixels for which the radial acute disambiguation solution was chosen
23	<	7: pixels for which the radial acute and NRWA disambiguation agree
21	>	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	>	coordinate bitmaps are interpolated.
24	>
25	>	In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
26	>	and the pixels that equal 33 or 44 in bitmap. Here are the definitions of the pixel values:
27	>
28	>	For conf_disambig:
29	>	50 : not all solutions agree (weak field method applied)
30	>	60 : not all solutions agree (weak field + annealed)
31	>	90 : all solutions agree (strong field + annealed)
32	>	0 : not disambiguated
33	>
34	>	For bitmap:
35	>	1  : weak field outside smooth bounding curve
36	>	2  : strong field outside smooth bounding curve
37	>	33 : weak field inside smooth bounding curve
38	>	34 : strong field inside smooth bounding curve
39
40		Written by Monica Bobra 15 August 2012
41		Potential Field code (appended) written by Keiji Hayashi
42
43		===========================================*/
44		#include <math.h>
45	+	#include <mkl.h>
46
47	<	#define PI              (M_PI)
47	>	#define PI       (M_PI)
48		#define MUNAUGHT (0.0000012566370614) /* magnetic constant */
49
50		/*===========================================*/
#	Line 51 | Line 67
67		//  5: pixels for which the radial acute disambiguation solution was chosen
68		//  7: pixels for which the radial acute and NRWA disambiguation agree
69
70	<	int computeAbsFlux(float *bz, int *dims, float *absFlux,
71	<	float *mean_vf_ptr, int *mask,
72	<	float cdelt1, double rsun_ref, double rsun_obs)
70	>	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	>	int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
73
74		{
75
76		int nx = dims[0], ny = dims[1];
77		int i, j, count_mask=0;
78	<	double sum=0.0;
78	>	double sum,err=0.0;
79
80		if (nx <= 0 || ny <= 0) return 1;
81
82		*absFlux = 0.0;
83		*mean_vf_ptr =0.0;
84
85	<	for (j = 0; j < ny; j++)
85	>	for (i = 0; i < nx; i++)
86		{
87	<	for (i = 0; i < nx; i++)
87	>	for (j = 0; j < ny; j++)
88		{
89	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
89	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
90	>	if isnan(bz[j * nx + i]) continue;
91		sum += (fabs(bz[j * nx + i]));
92	<	//sum += (fabs(bz[j * nx + i]))*inverseMu[j * nx + i]; // use this with mu function
92	>	err += bz_err[j * nx + i]*bz_err[j * nx + i];
93		count_mask++;
94		}
95		}
96
97	<	printf("nx=%d,ny=%d,count_mask=%d,sum=%f\n",nx,ny,count_mask,sum);
98	<	printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);
99	<	*mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
97	>	*mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
98	>	*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
99	>	*count_mask_ptr  = count_mask;
100	>	printf("CMASK=%g\n",*count_mask_ptr);
101	>	printf("USFLUX=%g\n",*mean_vf_ptr);
102	>	printf("sum=%f\n",sum);
103	>	printf("USFLUX_err=%g\n",*mean_vf_err_ptr);
104		return 0;
105		}
106
107		/*===========================================*/
108	<	/* Example function 2: Calculate Bh in units of Gauss */
108	>	/* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
109		// Native units of Bh are Gauss
110
111	<	int computeBh(float *bx, float *by, float *bz, float *bh, int *dims,
112	<	float *mean_hf_ptr, int *mask)
111	>	int computeBh(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
112	>	float *mean_hf_ptr, int *mask, int *bitmask)
113
114		{
115
116		int nx = dims[0], ny = dims[1];
117		int i, j, count_mask=0;
118		float sum=0.0;
119	<	*mean_hf_ptr =0.0;
119	>	*mean_hf_ptr = 0.0;
120
121		if (nx <= 0 || ny <= 0) return 1;
122
123	<	for (j = 0; j < ny; j++)
123	>	for (i = 0; i < nx; i++)
124		{
125	<	for (i = 0; i < nx; i++)
125	>	for (j = 0; j < ny; j++)
126		{
127	+	if isnan(bx[j * nx + i]) continue;
128	+	if isnan(by[j * nx + i]) continue;
129		bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
130		sum += bh[j * nx + i];
131	+	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];
132		count_mask++;
133		}
134		}
135
136		*mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
137	<	printf("*mean_hf_ptr=%f\n",*mean_hf_ptr);
137	>
138		return 0;
139		}
140
141		/*===========================================*/
142		/* Example function 3: Calculate Gamma in units of degrees */
143		// Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
144	+	// Redo calculation in radians for error analysis (since derivatives are only true in units of radians).
145
146	<
147	<	int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims,
123	<	float *mean_gamma_ptr, int *mask)
146	>	int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
147	>	float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)
148		{
149		int nx = dims[0], ny = dims[1];
150		int i, j, count_mask=0;
#	Line 128 | Line 152 | int computeGamma(float *bx, float *by, f
152		if (nx <= 0 || ny <= 0) return 1;
153
154		*mean_gamma_ptr=0.0;
155	<	float sum=0.0;
156	<	int count=0;
155	>	float sum,err,err_value=0.0;
156	>
157
158		for (i = 0; i < nx; i++)
159		{
#	Line 137 | Line 161 | int computeGamma(float *bx, float *by, f
161		{
162		if (bh[j * nx + i] > 100)
163		{
164	<	//if (mask[j * nx + i] != 90 ) continue;
165	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
166	<	sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
167	<	count_mask++;
164	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
165	>	if isnan(bz[j * nx + i]) continue;
166	>	if isnan(bz_err[j * nx + i]) continue;
167	>	if isnan(bh_err[j * nx + i]) continue;
168	>	if (bz[j * nx + i] == 0) continue;
169	>	sum += (atan(fabs(bz[j * nx + i]/bh[j * nx + i] )))*(180./PI);
170	>	err += (( sqrt ( ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) + ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bh[j * nx + i]*bh[j * nx + i])))  * fabs(bz[j * nx + i]/bh[j * nx + i]) ) / (1 + (bz[j * nx + i]/bh[j * nx + i])*(bz[j * nx + i]/bh[j * nx + i]))) *(180./PI);
171	>	count_mask++;
172		}
173		}
174		}
175
176		*mean_gamma_ptr = sum/count_mask;
177	<	printf("*mean_gamma_ptr=%f\n",*mean_gamma_ptr);
177	>	*mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.); // error in the quantity (sum)/(count_mask)
178	>	printf("MEANGAM=%f\n",*mean_gamma_ptr);
179	>	printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
180		return 0;
181		}
182
#	Line 154 | Line 184 | int computeGamma(float *bx, float *by, f
184		/* Example function 4: Calculate B_Total*/
185		// Native units of B_Total are in gauss
186
187	<	int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask)
187	>	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)
188		{
189
190		int nx = dims[0], ny = dims[1];
#	Line 166 | Line 196 | int computeB_total(float *bx, float *by,
196		{
197		for (j = 0; j < ny; j++)
198		{
199	+	if isnan(bx[j * nx + i]) continue;
200	+	if isnan(by[j * nx + i]) continue;
201	+	if isnan(bz[j * nx + i]) continue;
202		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]);
203	+	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];
204		}
205		}
206		return 0;
#	Line 175 | Line 209 | int computeB_total(float *bx, float *by,
209		/*===========================================*/
210		/* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
211
212	<	int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, float *derx_bt, float *dery_bt)
212	>	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)
213		{
214
215		int nx = dims[0], ny = dims[1];
#	Line 184 | Line 218 | int computeBtotalderivative(float *bt, i
218		if (nx <= 0 || ny <= 0) return 1;
219
220		*mean_derivative_btotal_ptr = 0.0;
221	<	float sum = 0.0;
221	>	float sum, err = 0.0;
222
223
224		/* brute force method of calculating the derivative (no consideration for edges) */
#	Line 232 | Line 266 | int computeBtotalderivative(float *bt, i
266		}
267
268
235	–	/* Just some print statements
236	–	for (i = 0; i < nx; i++)
237	–	{
238	–	for (j = 0; j < ny; j++)
239	–	{
240	–	printf("j=%d\n",j);
241	–	printf("i=%d\n",i);
242	–	printf("dery_bt[j*nx+i]=%f\n",dery_bt[j*nx+i]);
243	–	printf("derx_bt[j*nx+i]=%f\n",derx_bt[j*nx+i]);
244	–	printf("bt[j*nx+i]=%f\n",bt[j*nx+i]);
245	–	}
246	–	}
247	–	*/
248	–
269		for (i = 0; i <= nx-1; i++)
270		{
271		for (j = 0; j <= ny-1; j++)
272		{
273	<	// if ( (derx_bt[j * nx + i]-dery_bt[j * nx + i]) == 0) continue;
274	<	//if (mask[j * nx + i] != 90 ) continue;
275	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
273	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
274	>	if isnan(derx_bt[j * nx + i]) continue;
275	>	if isnan(dery_bt[j * nx + i]) continue;
276		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 */
277	+	err += (2.0)*bt_err[j * nx + i]*bt_err[j * nx + i];
278		count_mask++;
279		}
280		}
281
282	<	*mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
283	<	printf("*mean_derivative_btotal_ptr=%f\n",*mean_derivative_btotal_ptr);
282	>	*mean_derivative_btotal_ptr     = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
283	>	*mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
284	>	printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
285	>	printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
286		return 0;
287		}
288
#	Line 267 | Line 290 | int computeBtotalderivative(float *bt, i
290		/*===========================================*/
291		/* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
292
293	<	int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask, float *derx_bh, float *dery_bh)
293	>	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)
294		{
295
296		int nx = dims[0], ny = dims[1];
#	Line 276 | Line 299 | int computeBhderivative(float *bh, int *
299		if (nx <= 0 || ny <= 0) return 1;
300
301		*mean_derivative_bh_ptr = 0.0;
302	<	float sum = 0.0;
302	>	float sum,err = 0.0;
303
304		/* brute force method of calculating the derivative (no consideration for edges) */
305		for (i = 1; i <= nx-2; i++)
#	Line 323 | Line 346 | int computeBhderivative(float *bh, int *
346		}
347
348
326	–	/*Just some print statements
327	–	for (i = 0; i < nx; i++)
328	–	{
329	–	for (j = 0; j < ny; j++)
330	–	{
331	–	printf("j=%d\n",j);
332	–	printf("i=%d\n",i);
333	–	printf("dery_bh[j*nx+i]=%f\n",dery_bh[j*nx+i]);
334	–	printf("derx_bh[j*nx+i]=%f\n",derx_bh[j*nx+i]);
335	–	printf("bh[j*nx+i]=%f\n",bh[j*nx+i]);
336	–	}
337	–	}
338	–	*/
339	–
349		for (i = 0; i <= nx-1; i++)
350		{
351		for (j = 0; j <= ny-1; j++)
352		{
353	<	// if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;
354	<	//if (mask[j * nx + i] != 90 ) continue;
355	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
353	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
354	>	if isnan(derx_bh[j * nx + i]) continue;
355	>	if isnan(dery_bh[j * nx + i]) continue;
356		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 */
357	+	err += (2.0)*bh_err[j * nx + i]*bh_err[j * nx + i];
358		count_mask++;
359		}
360		}
361
362	<	*mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
362	>	*mean_derivative_bh_ptr     = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
363	>	*mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
364	>	printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
365	>	printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
366	>
367		return 0;
368		}
369
370		/*===========================================*/
371		/* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
372
373	<	int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask, float *derx_bz, float *dery_bz)
373	>	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)
374		{
375
376		int nx = dims[0], ny = dims[1];
#	Line 365 | Line 379 | int computeBzderivative(float *bz, int *
379		if (nx <= 0 || ny <= 0) return 1;
380
381		*mean_derivative_bz_ptr = 0.0;
382	<	float sum = 0.0;
382	>	float sum,err = 0.0;
383
384		/* brute force method of calculating the derivative (no consideration for edges) */
385		for (i = 1; i <= nx-2; i++)
386		{
387		for (j = 0; j <= ny-1; j++)
388		{
389	+	if isnan(bz[j * nx + i]) continue;
390		derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
391		}
392		}
#	Line 381 | Line 396 | int computeBzderivative(float *bz, int *
396		{
397		for (j = 1; j <= ny-2; j++)
398		{
399	+	if isnan(bz[j * nx + i]) continue;
400		dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
401		}
402		}
#	Line 390 | Line 406 | int computeBzderivative(float *bz, int *
406		i=0;
407		for (j = 0; j <= ny-1; j++)
408		{
409	+	if isnan(bz[j * nx + i]) continue;
410		derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
411		}
412
413		i=nx-1;
414		for (j = 0; j <= ny-1; j++)
415		{
416	+	if isnan(bz[j * nx + i]) continue;
417		derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
418		}
419
420		j=0;
421		for (i = 0; i <= nx-1; i++)
422		{
423	+	if isnan(bz[j * nx + i]) continue;
424		dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
425		}
426
427		j=ny-1;
428		for (i = 0; i <= nx-1; i++)
429		{
430	+	if isnan(bz[j * nx + i]) continue;
431		dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
432		}
433
434
415	–	/*Just some print statements
416	–	for (i = 0; i < nx; i++)
417	–	{
418	–	for (j = 0; j < ny; j++)
419	–	{
420	–	printf("j=%d\n",j);
421	–	printf("i=%d\n",i);
422	–	printf("dery_bz[j*nx+i]=%f\n",dery_bz[j*nx+i]);
423	–	printf("derx_bz[j*nx+i]=%f\n",derx_bz[j*nx+i]);
424	–	printf("bz[j*nx+i]=%f\n",bz[j*nx+i]);
425	–	}
426	–	}
427	–	*/
428	–
435		for (i = 0; i <= nx-1; i++)
436		{
437		for (j = 0; j <= ny-1; j++)
438		{
439		// if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;
440	<	//if (mask[j * nx + i] != 90 ) continue;
441	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
440	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
441	>	if isnan(bz[j * nx + i]) continue;
442	>	//if isnan(bz_err[j * nx + i]) continue;
443	>	if isnan(derx_bz[j * nx + i]) continue;
444	>	if isnan(dery_bz[j * nx + i]) continue;
445		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 */
446	+	err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];
447		count_mask++;
448		}
449		}
450
451		*mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
452	+	*mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
453	+	printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
454	+	printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
455	+
456		return 0;
457		}
458
459		/*===========================================*/
446	–
460		/* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */
461
462		//  In discretized space like data pixels,
#	Line 460 | Line 473 | int computeBzderivative(float *bz, int *
473		//
474		//  To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
475		//  one must perform the following unit conversions:
476	<	//  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
477	<	//  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),
476	>	//  (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
477	>	//  (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or
478	>	//  (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
479		//  where a Tesla is represented as a Newton/Ampere*meter.
480	+	//
481		//  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
482		//  In that case, we would have the following:
483		//  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
484		//  jz * (35.0)
485		//
486		//  The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
487	<	//  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(1000.)
488	<	//  =(Gauss/pix)(1/CDELT1)(0.0010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
474	<	//  =(Gauss/pix)(1/0.5)(10^-4)(4*PI*10^7)(722500)(1000.)
475	<	//  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
476	<
477	<	int computeJz(float *bx, float *by, int *dims, float *jz,
478	<	float *mean_jz_ptr, float *us_i_ptr, int *mask,
479	<	float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
487	>	//  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)
488	>	//  = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
489
490
491	+	// Comment out random number generator, which can only run on solar3
492	+	//int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
493	+	//            int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
494	+	//              float *noiseby, float *noisebz)
495
496	<	{
496	>	int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
497	>	int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
498
485	–	int nx = dims[0], ny = dims[1];
486	–	int i, j, count_mask=0;
499
500	<	if (nx <= 0 || ny <= 0) return 1;
500	>	{
501	>	int nx = dims[0], ny = dims[1];
502	>	int i, j, count_mask=0;
503	>
504	>	if (nx <= 0 || ny <= 0) return 1;
505	>	float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
506
490	–	*mean_jz_ptr = 0.0;
491	–	float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
492	–
507
508	+	/* Calculate the derivative*/
509		/* brute force method of calculating the derivative (no consideration for edges) */
510	+
511	+
512		for (i = 1; i <= nx-2; i++)
513		{
514		for (j = 0; j <= ny-1; j++)
515		{
516	+	if isnan(by[j * nx + i]) continue;
517		derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
518		}
519		}
520
503	–	/* brute force method of calculating the derivative (no consideration for edges) */
521		for (i = 0; i <= nx-1; i++)
522		{
523		for (j = 1; j <= ny-2; j++)
524		{
525	+	if isnan(bx[j * nx + i]) continue;
526		dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
527		}
528		}
529
530	<
513	<	/* consider the edges */
530	>	// consider the edges
531		i=0;
532		for (j = 0; j <= ny-1; j++)
533		{
534	+	if isnan(by[j * nx + i]) continue;
535		derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
536		}
537
538		i=nx-1;
539		for (j = 0; j <= ny-1; j++)
540		{
541	+	if isnan(by[j * nx + i]) continue;
542		derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
543	<	}
543	>	}
544
545		j=0;
546		for (i = 0; i <= nx-1; i++)
547		{
548	+	if isnan(bx[j * nx + i]) continue;
549		dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
550		}
551
552		j=ny-1;
553	<	for (i = 0; i <= nx-1; i++)
553	>	for (i = 0; i <= nx-1; i++)
554		{
555	+	if isnan(bx[j * nx + i]) continue;
556		dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
557		}
558
538	–	/* Just some print statements
539	–	for (i = 0; i < nx; i++)
540	–	{
541	–	for (j = 0; j < ny; j++)
542	–	{
543	–	printf("j=%d\n",j);
544	–	printf("i=%d\n",i);
545	–	printf("dery[j*nx+i]=%f\n",dery[j*nx+i]);
546	–	printf("derx[j*nx+i]=%f\n",derx[j*nx+i]);
547	–	printf("bx[j*nx+i]=%f\n",bx[j*nx+i]);
548	–	printf("by[j*nx+i]=%f\n",by[j*nx+i]);
549	–	}
550	–	}
551	–	*/
552	–
559
560		for (i = 0; i <= nx-1; i++)
561		{
562		for (j = 0; j <= ny-1; j++)
563		{
564	<	// if ( (derx[j * nx + i]-dery[j * nx + i]) == 0) continue;
565	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
566	<	//if (mask[j * nx + i] != 90 ) continue;
567	<	curl +=     (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
568	<	us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */
569	<	jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */
564	>	// calculate jz at all points
565	>	jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix
566	>	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]) +
567	>	(by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;
568	>	jz_err_squared[j * nx + i]=(jz_err[j * nx + i]*jz_err[j * nx + i]);
569	>	count_mask++;
570	>	}
571	>	}
572	>
573	>	return 0;
574	>	}
575
576	+	/*===========================================*/
577	+
578	+
579	+	/* Example function 9:  Compute quantities on Jz array */
580	+	// Compute mean and total current on Jz array.
581	+
582	+	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,
583	+	float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
584	+	float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
585	+
586	+	{
587	+
588	+	int nx = dims[0], ny = dims[1];
589	+	int i, j, count_mask=0;
590	+
591	+	if (nx <= 0 || ny <= 0) return 1;
592	+
593	+	float curl,us_i,test_perimeter,mean_curl,err=0.0;
594	+
595	+
596	+	/* 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*/
597	+	for (i = 0; i <= nx-1; i++)
598	+	{
599	+	for (j = 0; j <= ny-1; j++)
600	+	{
601	+	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
602	+	if isnan(derx[j * nx + i]) continue;
603	+	if isnan(dery[j * nx + i]) continue;
604	+	if isnan(jz[j * nx + i]) continue;
605	+	curl +=     (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
606	+	us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */
607	+	err  += (jz_err[j * nx + i]*jz_err[j * nx + i]);
608		count_mask++;
609	<	}
609	>	}
610		}
611
612	<	mean_curl        = (curl/count_mask);
570	<	printf("mean_curl=%f\n",mean_curl);
571	<	printf("cdelt1, what is it?=%f\n",cdelt1);
612	>	/* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
613		*mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */
614	<	printf("count_mask=%d\n",count_mask);
615	<	*us_i_ptr        = (us_i);                   /* us_i gets populated as MEANJZD */
614	>	*mean_jz_err_ptr = (sqrt(err))*fabs(((rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.))/(count_mask)); // error in the quantity MEANJZD
615	>
616	>	*us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */
617	>	*us_i_err_ptr    = (sqrt(err))*fabs((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ
618	>
619	>	printf("MEANJZD=%f\n",*mean_jz_ptr);
620	>	printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
621	>
622	>	printf("TOTUSJZ=%g\n",*us_i_ptr);
623	>	printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
624	>
625		return 0;
626
627		}
628
579	–
629		/*===========================================*/
581	–	/* Example function 9:  Twist Parameter, alpha */
630
631	<	// The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm
631	>	/* Example function 10:  Twist Parameter, alpha */
632	>
633	>	// The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
634	>	// for alpha is calculated in the following way (different from Leka and Barnes' approach):
635	>
636	>	// (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
637	>	// (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
638	>	// avg alpha = avg Jz / avg Bz
639	>
640	>	// The sign is assigned as follows:
641	>	// If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.
642	>	// If this value is > 0, then alpha is > 0.
643	>	// If this value is < 0, then alpha is <0.
644	>	//
645	>	// If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.
646	>	// If this value is > 0, then alpha is < 0.
647	>	// If this value is < 0, then alpha is > 0.
648	>
649	>	// The units of alpha are in 1/Mm
650		// The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
651		//
652		// Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
653		//                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
654		//                               = 1/Mm
655
656	<	int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask,
657	<	float cdelt1, double rsun_ref, double rsun_obs)
656	>	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)
657	>
658		{
659		int nx = dims[0], ny = dims[1];
660	<	int i, j, count_mask=0;
660	>	int i, j, count_mask, a,b,c,d=0;
661
662		if (nx <= 0 || ny <= 0) return 1;
663
664	<	*mean_alpha_ptr = 0.0;
599	<	float aa, bb, cc, bznew, alpha2, sum=0.0;
664	>	float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6, sum_err=0.0;
665
666		for (i = 1; i < nx-1; i++)
667		{
668		for (j = 1; j < ny-1; j++)
669		{
670	<	//if (mask[j * nx + i] != 90 ) continue;
671	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
670	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
671	>	//if isnan(jz_smooth[j * nx + i]) continue;
672		if isnan(jz[j * nx + i]) continue;
673		if isnan(bz[j * nx + i]) continue;
674	<	if (bz[j * nx + i] == 0.0) continue;
675	<	sum += (jz[j * nx + i] / bz[j * nx + i])*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */
674	>	//if (jz_smooth[j * nx + i] == 0) continue;
675	>	if (jz[j * nx + i]     == 0.0) continue;
676	>	if (bz_err[j * nx + i] == 0.0) continue;
677	>	if (bz[j * nx + i]     == 0.0) continue;
678	>	if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i] ); a++;
679	>	if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;
680	>	//if (bz[j * nx + i] >  0) sum3 += ( jz_smooth[j * nx + i]);
681	>	//if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
682	>	if (bz[j * nx + i] >  0) sum3 += ( jz[j * nx + i] ); c++;
683	>	if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;
684	>	sum5    += bz[j * nx + i];
685	>	/* sum_err is a fractional uncertainty */
686	>	sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs( ( (jz[j * nx + i]) / (bz[j * nx + i]) ) *(1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
687		count_mask++;
688		}
689		}
690	+
691	+	sum     = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)); /* the units for (jz/bz) are 1/Mm */
692	+
693	+	/* Determine the sign of alpha */
694	+	if ((sum5 > 0) && (sum3 >  0)) sum=sum;
695	+	if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
696	+	if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
697	+	if ((sum5 < 0) && (sum4 >  0)) sum=-sum;
698	+
699	+	*mean_alpha_ptr = sum; /* Units are 1/Mm */
700	+	*mean_alpha_err_ptr    = (sqrt(sum_err*sum_err)) / ((a+b+c+d)*100.); // error in the quantity (sum)/(count_mask); factor of 100 comes from converting percent
701	+
702	+	printf("a=%d\n",a);
703	+	printf("b=%d\n",b);
704	+	printf("d=%d\n",d);
705	+	printf("c=%d\n",c);
706	+
707	+	printf("MEANALP=%f\n",*mean_alpha_ptr);
708	+	printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);
709
615	–	printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);
616	–	printf("count_mask=%d\n",count_mask);
617	–	printf("sum=%f\n",sum);
618	–	*mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */
710		return 0;
711		}
712
713		/*===========================================*/
714	<	/* Example function 10:  Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */
714	>	/* Example function 11:  Helicity (mean current helicty, total unsigned current helicity, absolute value of net current helicity) */
715
716		//  The current helicity is defined as Bz*Jz and the units are G^2 / m
717		//  The units of Jz are in G/pix; the units of Bz are in G.
718	<	//  Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/m)
718	>	//  Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter)
719		//                                                      = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
720	<	//                                                      = G^2 / m.
630	<
720	>	//                                                      =  G^2 / m.
721
722	<	int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr,
723	<	float *total_abs_ih_ptr, int *mask, float cdelt1, double rsun_ref, double rsun_obs)
722	>	int computeHelicity(float *jz_err, float *jz_rms_err, float *bz_err, float *bz, int *dims, float *jz, float *mean_ih_ptr,
723	>	float *mean_ih_err_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr,
724	>	float *total_us_ih_err_ptr, float *total_abs_ih_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
725
726		{
727
#	Line 639 | Line 730 | int computeHelicity(float *bz, int *dims
730
731		if (nx <= 0 || ny <= 0) return 1;
732
733	<	*mean_ih_ptr = 0.0;
643	<	float sum=0.0, sum2=0.0;
733	>	float sum,sum2,sum_err=0.0;
734
735	<	for (j = 0; j < ny; j++)
735	>	for (i = 0; i < nx; i++)
736		{
737	<	for (i = 0; i < nx; i++)
737	>	for (j = 0; j < ny; j++)
738		{
739	<	//if (mask[j * nx + i] != 90 ) continue;
740	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
741	<	if isnan(jz[j * nx + i]) continue;
742	<	if isnan(bz[j * nx + i]) continue;
743	<	if (bz[j * nx + i] == 0.0) continue;
744	<	if (jz[j * nx + i] == 0.0) continue;
745	<	sum  +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
746	<	sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
747	<	count_mask++;
739	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
740	>	if isnan(jz[j * nx + i]) continue;
741	>	if isnan(bz[j * nx + i]) continue;
742	>	if (bz[j * nx + i] == 0.0) continue;
743	>	if (jz[j * nx + i] == 0.0) continue;
744	>	sum     +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
745	>	sum2    += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
746	>	sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs(jz[j * nx + i]*bz[j * nx + i]*(1/cdelt1)*(rsun_obs/rsun_ref));
747	>	count_mask++;
748		}
749		}
750
751	<	printf("sum/count_mask=%f\n",sum/count_mask);
752	<	printf("(1/cdelt1)*(rsun_obs/rsun_ref)=%f\n",(1/cdelt1)*(rsun_obs/rsun_ref));
753	<	*mean_ih_ptr     = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */
754	<	*total_us_ih_ptr = sum2;           /* Units are G^2 / m */
755	<	*total_abs_ih_ptr= fabs(sum);      /* Units are G^2 / m */
751	>	*mean_ih_ptr          = sum/count_mask ; /* Units are G^2 / m ; keyword is MEANJZH */
752	>	*total_us_ih_ptr      = sum2           ; /* Units are G^2 / m ; keyword is TOTUSJH */
753	>	*total_abs_ih_ptr     = fabs(sum)      ; /* Units are G^2 / m ; keyword is ABSNJZH */
754	>
755	>	*mean_ih_err_ptr      = (sqrt(sum_err*sum_err)) / (count_mask*100.)    ;  // error in the quantity MEANJZH
756	>	*total_us_ih_err_ptr  = (sqrt(sum_err*sum_err)) / (100.)               ;  // error in the quantity TOTUSJH
757	>	*total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.)               ;  // error in the quantity ABSNJZH
758	>
759	>	printf("MEANJZH=%f\n",*mean_ih_ptr);
760	>	printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
761	>
762	>	printf("TOTUSJH=%f\n",*total_us_ih_ptr);
763	>	printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
764	>
765	>	printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
766	>	printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
767
768		return 0;
769		}
770
771		/*===========================================*/
772	<	/* Example function 11:  Sum of Absolute Value per polarity  */
772	>	/* Example function 12:  Sum of Absolute Value per polarity  */
773
774		//  The Sum of the Absolute Value per polarity is defined as the following:
775		//  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
776		//  The units of jz are in G/pix. In this case, we would have the following:
777		//  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
778		//     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
779	+	//
780	+	//  The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
781
782	<	int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,
783	<	int *mask, float cdelt1, double rsun_ref, double rsun_obs)
782	>	int computeSumAbsPerPolarity(float *jz_err, float *bz_err, float *bz, float *jz, int *dims, float *totaljzptr, float *totaljz_err_ptr,
783	>	int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
784
785		{
786		int nx = dims[0], ny = dims[1];
#	Line 686 | Line 789 | int computeSumAbsPerPolarity(float *bz,
789		if (nx <= 0 || ny <= 0) return 1;
790
791		*totaljzptr=0.0;
792	<	float sum1=0.0, sum2=0.0;
792	>	float sum1,sum2,err=0.0;
793
794		for (i = 0; i < nx; i++)
795		{
796		for (j = 0; j < ny; j++)
797		{
798	<	//if (mask[j * nx + i] != 90 ) continue;
799	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
798	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
799	>	if isnan(bz[j * nx + i]) continue;
800		if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
801		if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
802	+	err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
803	+	count_mask++;
804		}
805		}
806
807	<	*totaljzptr = fabs(sum1) + fabs(sum2);  /* Units are A */
807	>	*totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */
808	>	*totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
809	>	printf("SAVNCPP=%g\n",*totaljzptr);
810	>	printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
811	>
812		return 0;
813		}
814
815		/*===========================================*/
816	<	/* Example function 12:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
816	>	/* Example function 13:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
817		// The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
818	<	// automatically yields erg per cubic centimeter for an input B in Gauss.
818	>	// automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8*PI can come out of the integral; thus,
819	>	// the integral is over B^2 dV and the 8*PI is divided at the end.
820		//
821		// Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
822		// ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
823	<	// erg/cm^3(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.)^2
824	<	// = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
825	<	// = erg/cm^3(1.30501e15)
823	>	//   erg/cm^3*(CDELT1^2)*(RSUN_REF/RSUN_OBS ^2)*(100.^2)
824	>	// = erg/cm^3*(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
825	>	// = erg/cm^3*(1.30501e15)
826		// = erg/cm(1/pix^2)
827
828	<	int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims,
829	<	float *meanpotptr, float *totpotptr, int *mask,
828	>	int computeFreeEnergy(float *bx_err, float *by_err, float *bx, float *by, float *bpx, float *bpy, int *dims,
829	>	float *meanpotptr, float *meanpot_err_ptr, float *totpotptr, float *totpot_err_ptr, int *mask, int *bitmask,
830		float cdelt1, double rsun_ref, double rsun_obs)
831
832		{
#	Line 727 | Line 837 | int computeFreeEnergy(float *bx, float *
837
838		*totpotptr=0.0;
839		*meanpotptr=0.0;
840	<	float sum=0.0;
840	>	float sum,sum1,err=0.0;
841
842		for (i = 0; i < nx; i++)
843		{
844		for (j = 0; j < ny; j++)
845		{
846	<	//if (mask[j * nx + i] != 90 ) continue;
847	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
848	<	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);
846	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
847	>	if isnan(bx[j * nx + i]) continue;
848	>	if isnan(by[j * nx + i]) continue;
849	>	sum  += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
850	>	sum1 += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) );
851	>	err  += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);
852		count_mask++;
853		}
854		}
855
856	<	*meanpotptr = (sum)/(count_mask);              /* Units are ergs per cubic centimeter */
857	<	*totpotptr  = sum*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0)*(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 */
856	>	*meanpotptr      = (sum1/(8.*PI)) / (count_mask);     /* Units are ergs per cubic centimeter */
857	>	*meanpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
858	>
859	>	/* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
860	>	*totpotptr       = (sum)/(8.*PI);
861	>	*totpot_err_ptr  = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
862	>
863	>	printf("MEANPOT=%g\n",*meanpotptr);
864	>	printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
865	>
866	>	printf("TOTPOT=%g\n",*totpotptr);
867	>	printf("TOTPOT_err=%g\n",*totpot_err_ptr);
868	>
869		return 0;
870		}
871
872		/*===========================================*/
873	<	/* 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 */
873	>	/* 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 */
874
875	<	int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
876	<	float *meanshear_angleptr, float *area_w_shear_gt_45ptr,
753	<	float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,
754	<	int *mask)
875	>	int computeShearAngle(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
876	>	float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
877		{
878		int nx = dims[0], ny = dims[1];
879		int i, j;
880
881		if (nx <= 0 || ny <= 0) return 1;
882
883	<	*area_w_shear_gt_45ptr=0.0;
884	<	*meanshear_angleptr=0.0;
885	<	float dotproduct, magnitude_potential, magnitude_vector, shear_angle=0.0, sum = 0.0, count=0.0, count_mask=0.0;
764	<	float dotproducth, magnitude_potentialh, magnitude_vectorh, shear_angleh=0.0, sum1 = 0.0, counth = 0.0;
883	>	//*area_w_shear_gt_45ptr=0.0;
884	>	//*meanshear_angleptr=0.0;
885	>	float dotproduct, magnitude_potential, magnitude_vector, shear_angle,err=0.0, sum = 0.0, count=0.0, count_mask=0.0;
886
887		for (i = 0; i < nx; i++)
888		{
889		for (j = 0; j < ny; j++)
890		{
891	<	//if (mask[j * nx + i] != 90 ) continue;
771	<	if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;
891	>	if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
892		if isnan(bpx[j * nx + i]) continue;
893		if isnan(bpy[j * nx + i]) continue;
894		if isnan(bpz[j * nx + i]) continue;
895		if isnan(bz[j * nx + i]) continue;
896	+	if isnan(bx[j * nx + i]) continue;
897	+	if isnan(by[j * nx + i]) continue;
898		/* For mean 3D shear angle, area with shear greater than 45*/
899		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]);
900	<	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]));
901	<	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]) );
900	>	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]));
901	>	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]) );
902		shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
903		count ++;
904		sum += shear_angle ;
905	+	err += -(1./(1.- sqrt(bx_err[j * nx + i]*bx_err[j * nx + i]+by_err[j * nx + i]*by_err[j * nx + i]+bh_err[j * nx + i]*bh_err[j * nx + i])));
906		if (shear_angle > 45) count_mask ++;
907		}
908		}
909
910		/* For mean 3D shear angle, area with shear greater than 45*/
911	<	*meanshear_angleptr = (sum)/(count);              /* Units are degrees */
912	<	printf("count_mask=%f\n",count_mask);
913	<	printf("count=%f\n",count);
914	<	*area_w_shear_gt_45ptr = (count_mask/(count))*(100.);  /* The area here is a fractional area -- the % of the total area */
911	>	*meanshear_angleptr = (sum)/(count);                 /* Units are degrees */
912	>	*meanshear_angle_err_ptr = (sqrt(err*err))/(count);  // error in the quantity (sum)/(count_mask)
913	>	*area_w_shear_gt_45ptr   = (count_mask/(count))*(100.);/* The area here is a fractional area -- the % of the total area */
914	>
915	>	printf("MEANSHR=%f\n",*meanshear_angleptr);
916	>	printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
917
918		return 0;
919		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines