ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing proj/sharp/apps/smarp_functions.c (file contents):
Revision 1.2 by mbobra, Mon May 4 21:11:07 2020 UTC vs.
Revision 1.3 by mbobra, Mon Jun 29 22:19:51 2020 UTC

#	Line 1 | Line 1
1
2		/*===========================================
3
4	<	The following 3 functions calculate the following spaceweather indices:
4	>	The following 3 functions calculate the following spaceweather indices on line-of-sight
5	>	magnetic field data:
6
7	<	USFLUX Total unsigned flux in Maxwells
8	<	MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
9	<	R_VALUE Karel Schrijver's R parameter
7	>	USFLUX  Total unsigned flux in Maxwells
8	>	MEANGBZ Mean value of the line-of-sight (approximately vertical in remapped and reprojected data) field gradient, in Gauss/Mm
9	>	R_VALUE R parameter (Schrijver, 2007)
10
11		The indices are calculated on the pixels in which the bitmap segment is greater than 36.
12
13		The SMARP bitmap has 13 unique values because they describe three different characteristics:
14		the location of the pixel (whether the pixel is off disk or a member of the patch), the
15		character of the magnetic field (quiet or active), and the character of the continuum
16	<	intensity data (quiet, part of faculae, part of a sunspot).
16	>	intensity data (quiet, faculae, sunspot).
17
18		Here are all the possible values:
19
#	Line 58 | Line 59
59		/* Example function 1: Compute total unsigned flux in units of G/cm^2 */
60
61		//  To compute the unsigned flux, we simply calculate
62	<	//  flux = surface integral [(vector Bz) dot (normal vector)],
63	<	//       = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
62	>	//  flux = surface integral [(vector LOS) dot (normal vector)],
63	>	//       = surface integral [(magnitude LOS)*(magnitude normal)*(cos theta)].
64		//  However, since the field is radial, we will assume cos theta = 1.
65		//  Therefore the pixels only need to be corrected for the projection.
66
#	Line 68 | Line 69
69		//  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
70		//  =Gauss*cm^2
71
72	<	int computeAbsFlux(float *bz, int *dims, float *absFlux,
73	<	float *mean_vf_ptr, float *count_mask_ptr,
74	<	int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
72	>	int computeAbsFlux_los(float *los, int *dims, float *absFlux,
73	>	float *mean_vf_ptr, float *count_mask_ptr,
74	>	int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
75
76		{
77
#	Line 91 | Line 92 | int computeAbsFlux(float *bz, int *dims,
92		for (j = 0; j < ny; j++)
93		{
94		if ( bitmask[j * nx + i] < 36 ) continue;
95	<	if isnan(bz[j * nx + i]) continue;
96	<	sum += (fabs(bz[j * nx + i]));
95	>	if isnan(los[j * nx + i]) continue;
96	>	sum += (fabs(los[j * nx + i]));
97		count_mask++;
98		}
99		}
#	Line 109 | Line 110 | int computeAbsFlux(float *bz, int *dims,
110		}
111
112		/*===========================================*/
113	<	/* Example function 2:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
113	>	/* Example function 2:  Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */
114
115	<	int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *bitmask, float *derx_bz, float *dery_bz)
115	>	int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los)
116		{
117
118		int nx = dims[0];
#	Line 120 | Line 121 | int computeBzderivative(float *bz, int *
121		int j = 0;
122		int count_mask = 0;
123		double sum = 0.0;
124	<	*mean_derivative_bz_ptr = 0.0;
124	>	*mean_derivative_los_ptr = 0.0;
125
126		if (nx <= 0 || ny <= 0) return 1;
127
#	Line 129 | Line 130 | int computeBzderivative(float *bz, int *
130		{
131		for (j = 0; j <= ny-1; j++)
132		{
133	<	derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
133	>	derx_los[j * nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])*0.5;
134		}
135		}
136
#	Line 138 | Line 139 | int computeBzderivative(float *bz, int *
139		{
140		for (j = 1; j <= ny-2; j++)
141		{
142	<	dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
142	>	dery_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
143		}
144		}
145
#	Line 148 | Line 149 | int computeBzderivative(float *bz, int *
149		i=0;
150		for (j = 0; j <= ny-1; j++)
151		{
152	<	derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
152	>	derx_los[j * nx + i] = ( (-3*los[j * nx + i]) + (4*los[j * nx + (i+1)]) - (los[j * nx + (i+2)]) )*0.5;
153		}
154
155		i=nx-1;
156		for (j = 0; j <= ny-1; j++)
157		{
158	<	derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
158	>	derx_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[j * nx + (i-1)]) - (-los[j * nx + (i-2)]) )*0.5;
159		}
160
161		j=0;
162		for (i = 0; i <= nx-1; i++)
163		{
164	<	dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
164	>	dery_los[j * nx + i] = ( (-3*los[j*nx + i]) + (4*los[(j+1) * nx + i]) - (los[(j+2) * nx + i]) )*0.5;
165		}
166
167		j=ny-1;
168		for (i = 0; i <= nx-1; i++)
169		{
170	<	dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
170	>	dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
171		}
172
173
#	Line 175 | Line 176 | int computeBzderivative(float *bz, int *
176		for (j = 0; j <= ny-1; j++)
177		{
178		if ( bitmask[j * nx + i] < 36 ) continue;
179	<	if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
180	<	if isnan(bz[j * nx + i])      continue;
181	<	if isnan(bz[(j+1) * nx + i])  continue;
182	<	if isnan(bz[(j-1) * nx + i])  continue;
183	<	if isnan(bz[j * nx + i-1])    continue;
184	<	if isnan(bz[j * nx + i+1])    continue;
185	<	if isnan(derx_bz[j * nx + i]) continue;
186	<	if isnan(dery_bz[j * nx + i]) continue;
187	<	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 */
179	>	if ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) continue;
180	>	if isnan(los[j * nx + i])      continue;
181	>	if isnan(los[(j+1) * nx + i])  continue;
182	>	if isnan(los[(j-1) * nx + i])  continue;
183	>	if isnan(los[j * nx + i-1])    continue;
184	>	if isnan(los[j * nx + i+1])    continue;
185	>	if isnan(derx_los[j * nx + i]) continue;
186	>	if isnan(dery_los[j * nx + i]) continue;
187	>	sum += sqrt( derx_los[j * nx + i]*derx_los[j * nx + i]  + dery_los[j * nx + i]*dery_los[j * nx + i] ); /* Units of Gauss */
188		count_mask++;
189		}
190		}
191
192	<	*mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
193	<	//printf("mean_derivative_bz_ptr=%f\n",*mean_derivative_bz_ptr);
192	>	*mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
193	>	//printf("mean_derivative_los_ptr=%f\n",*mean_derivative_los_ptr);
194
195		return 0;
196		}
#	Line 205 | Line 206 | int computeBzderivative(float *bz, int *
206		// are the dimensions of the input array, fsample() will usually result
207		// in a segfault (though not always, depending on how the segfault was accessed.)
208
209	<	int computeR(float *los, int *dims, float *Rparam, float cdelt1,
210	<	float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
211	<	float *pmap, int nx1, int ny1,
212	<	int scale, float *p1pad, int nxp, int nyp, float *pmapn)
209	>	int computeR_los(float *los, int *dims, float *Rparam, float cdelt1,
210	>	float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
211	>	float *pmap, int nx1, int ny1,
212	>	int scale, float *p1pad, int nxp, int nyp, float *pmapn)
213
214		{
215		int nx = dims[0];

Diff Legend

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