1 mbobra 1.1
2 /*===========================================
3
4 The following 3 functions calculate the following spaceweather indices:
5
6 USFLUX Total unsigned flux in Maxwells
7 MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
8 R_VALUE Karel Schrijver's R parameter
9
10 The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
11 pixels in which the bitmap segment is greater than 42. These ranges are selected because the CCD
12 coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data
13 prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI
14 contain nearest-neighbor bitmaps).
15
16 In the CCD coordinates, this means that we are selecting the pixels that that equal 33 or 34 in bitmap. Here are the definitions of the pixel values:
17
18
19 For bitmap:
20 1 : weak field outside smooth bounding curve
21 2 : strong field outside smooth bounding curve
22 mbobra 1.1 33 : weak field inside smooth bounding curve
23 34 : strong field inside smooth bounding curve
24
25 Written by Monica Bobra 15 August 2012
26 Potential Field code (appended) written by Keiji Hayashi
27 Error analysis modification 21 October 2013
28
29 ===========================================*/
30 #include <math.h>
31 #include <mkl.h>
32
33 #define PI (M_PI)
34 #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
35
36 /*===========================================*/
37
38 /* Example function 1: Compute total unsigned flux in units of G/cm^2 */
39
40 // To compute the unsigned flux, we simply calculate
41 // flux = surface integral [(vector Bz) dot (normal vector)],
42 // = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
43 mbobra 1.1 // However, since the field is radial, we will assume cos theta = 1.
44 // Therefore the pixels only need to be corrected for the projection.
45
46 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
47 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
48 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
49 // =Gauss*cm^2
50
51 int computeAbsFlux(float *bz, int *dims, float *absFlux,
52 float *mean_vf_ptr, float *count_mask_ptr,
53 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
54
55 {
56
57 int nx = dims[0];
58 int ny = dims[1];
59 int i = 0;
60 int j = 0;
61 int count_mask = 0;
62 double sum = 0.0;
63 *absFlux = 0.0;
64 mbobra 1.1 *mean_vf_ptr = 0.0;
65
66
67 if (nx <= 0 || ny <= 0) return 1;
68
69 for (i = 0; i < nx; i++)
70 {
71 for (j = 0; j < ny; j++)
72 {
73 if ( bitmask[j * nx + i] < 42 ) continue;
74 if isnan(bz[j * nx + i]) continue;
75 sum += (fabs(bz[j * nx + i]));
76 count_mask++;
77 }
78 }
79
80 *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
81 *count_mask_ptr = count_mask;
82 //printf("mean_vf_ptr=%f\n",*mean_vf_ptr);
83 //printf("count_mask_ptr=%f\n",*count_mask_ptr);
84 //printf("cdelt1=%f\n",cdelt1);
85 mbobra 1.1 //printf("rsun_ref=%f\n",rsun_ref);
86 //printf("rsun_obs=%f\n",rsun_obs);
87
88 return 0;
89 }
90
91 /*===========================================*/
92 /* Example function 2: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
93
94 int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *bitmask, float *derx_bz, float *dery_bz)
95 {
96
97 int nx = dims[0];
98 int ny = dims[1];
99 int i = 0;
100 int j = 0;
101 int count_mask = 0;
102 double sum = 0.0;
103 *mean_derivative_bz_ptr = 0.0;
104
105 if (nx <= 0 || ny <= 0) return 1;
106 mbobra 1.1
107 /* brute force method of calculating the derivative (no consideration for edges) */
108 for (i = 1; i <= nx-2; i++)
109 {
110 for (j = 0; j <= ny-1; j++)
111 {
112 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
113 }
114 }
115
116 /* brute force method of calculating the derivative (no consideration for edges) */
117 for (i = 0; i <= nx-1; i++)
118 {
119 for (j = 1; j <= ny-2; j++)
120 {
121 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
122 }
123 }
124
125 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
126 ignore the edges for the error terms as those arrays have been initialized to zero.
127 mbobra 1.1 this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
128 i=0;
129 for (j = 0; j <= ny-1; j++)
130 {
131 derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
132 }
133
134 i=nx-1;
135 for (j = 0; j <= ny-1; j++)
136 {
137 derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
138 }
139
140 j=0;
141 for (i = 0; i <= nx-1; i++)
142 {
143 dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
144 }
145
146 j=ny-1;
147 for (i = 0; i <= nx-1; i++)
148 mbobra 1.1 {
149 dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
150 }
151
152
153 for (i = 0; i <= nx-1; i++)
154 {
155 for (j = 0; j <= ny-1; j++)
156 {
157 if ( bitmask[j * nx + i] < 42 ) continue;
158 if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
159 if isnan(bz[j * nx + i]) continue;
160 if isnan(bz[(j+1) * nx + i]) continue;
161 if isnan(bz[(j-1) * nx + i]) continue;
162 if isnan(bz[j * nx + i-1]) continue;
163 if isnan(bz[j * nx + i+1]) continue;
164 if isnan(derx_bz[j * nx + i]) continue;
165 if isnan(dery_bz[j * nx + i]) continue;
166 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 */
167 count_mask++;
168 }
169 mbobra 1.1 }
170
171 *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
172 //printf("mean_derivative_bz_ptr=%f\n",*mean_derivative_bz_ptr);
173
174 return 0;
175 }
176
177 /*===========================================*/
178 /* Example function 3: R parameter as defined in Schrijver, 2007 */
179 //
180 // Note that there is a restriction on the function fsample()
181 // If the following occurs:
182 // nx_out > floor((ny_in-1)/scale + 1)
183 // ny_out > floor((ny_in-1)/scale + 1),
184 // where n*_out are the dimensions of the output array and n*_in
185 // are the dimensions of the input array, fsample() will usually result
186 // in a segfault (though not always, depending on how the segfault was accessed.)
187
188 int computeR(float *los, int *dims, float *Rparam, float cdelt1,
189 float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
190 mbobra 1.1 float *pmap, int nx1, int ny1,
191 int scale, float *p1pad, int nxp, int nyp, float *pmapn)
192
193 {
194 int nx = dims[0];
195 int ny = dims[1];
196 int i = 0;
197 int j = 0;
198 int index, index1;
199 double sum = 0.0;
200 *Rparam = 0.0;
201 struct fresize_struct fresboxcar, fresgauss;
202 struct fint_struct fints;
203 float sigma = 10.0/2.3548;
204
205 // set up convolution kernels
206 init_fresize_boxcar(&fresboxcar,1,1);
207 init_fresize_gaussian(&fresgauss,sigma,20,1);
208
209 // =============== [STEP 1] ===============
210 // bin the line-of-sight magnetogram down by a factor of scale
211 mbobra 1.1 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
212
213 // =============== [STEP 2] ===============
214 // identify positive and negative pixels greater than +/- 150 gauss
215 // and label those pixels with a 1.0 in arrays p1p0 and p1n0
216 for (i = 0; i < nx1; i++)
217 {
218 for (j = 0; j < ny1; j++)
219 {
220 index = j * nx1 + i;
221 if (rim[index] > 150)
222 p1p0[index]=1.0;
223 else
224 p1p0[index]=0.0;
225 if (rim[index] < -150)
226 p1n0[index]=1.0;
227 else
228 p1n0[index]=0.0;
229 }
230 }
231
232 mbobra 1.1 // =============== [STEP 3] ===============
233 // smooth each of the negative and positive pixel bitmaps
234 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
235 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
236
237 // =============== [STEP 4] ===============
238 // find the pixels for which p1p and p1n are both equal to 1.
239 // this defines the polarity inversion line
240 for (i = 0; i < nx1; i++)
241 {
242 for (j = 0; j < ny1; j++)
243 {
244 index = j * nx1 + i;
245 if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
246 p1[index]=1.0;
247 else
248 p1[index]=0.0;
249 }
250 }
251
252 // pad p1 with zeroes so that the gaussian colvolution in step 5
253 mbobra 1.1 // does not cut off data within hwidth of the edge
254
255 // step i: zero p1pad
256 for (i = 0; i < nxp; i++)
257 {
258 for (j = 0; j < nyp; j++)
259 {
260 index = j * nxp + i;
261 p1pad[index]=0.0;
262 }
263 }
264
265 // step ii: place p1 at the center of p1pad
266 for (i = 0; i < nx1; i++)
267 {
268 for (j = 0; j < ny1; j++)
269 {
270 index = j * nx1 + i;
271 index1 = (j+20) * nxp + (i+20);
272 p1pad[index1]=p1[index];
273 }
274 mbobra 1.1 }
275
276 // =============== [STEP 5] ===============
277 // convolve the polarity inversion line map with a gaussian
278 // to identify the region near the plarity inversion line
279 // the resultant array is called pmap
280 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0);
281
282
283 // select out the nx1 x ny1 non-padded array within pmap
284 for (i = 0; i < nx1; i++)
285 {
286 for (j = 0; j < ny1; j++)
287 {
288 index = j * nx1 + i;
289 index1 = (j+20) * nxp + (i+20);
290 pmapn[index]=pmap[index1];
291 }
292 }
293
294 // =============== [STEP 6] ===============
295 mbobra 1.1 // the R parameter is calculated
296 for (i = 0; i < nx1; i++)
297 {
298 for (j = 0; j < ny1; j++)
299 {
300 index = j * nx1 + i;
301 if isnan(pmapn[index]) continue;
302 if isnan(rim[index]) continue;
303 sum += pmapn[index]*abs(rim[index]);
304 }
305 }
306
307 if (sum < 1.0)
308 *Rparam = 0.0;
309 else
310 *Rparam = log10(sum);
311
312 //printf("R_VALUE=%f\n",*Rparam);
313
314 free_fresize(&fresboxcar);
315 free_fresize(&fresgauss);
316 mbobra 1.1
317 return 0;
318
319 }
320
321 /*===========================================*/
322
323 char *smarp_functions_version() // Returns CVS version of smarp_functions.c
324 {
325 return strdup("$Id");
326 }
|
Return to
smarp_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps