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 mbobra 1.2  The indices are calculated on the pixels in which the bitmap segment is greater than 36. 

 11 mbobra 1.1  

 12 mbobra 1.2  The SMARP bitmap has 13 unique values because they describe three different characteristics: 
 13             the location of the pixel (whether the pixel is off disk or a member of the patch), the 
 14             character of the magnetic field (quiet or active), and the character of the continuum 
 15             intensity data (quiet, part of faculae, part of a sunspot). 
 16            
 17             Here are all the possible values:

 18 mbobra 1.1  

 19 mbobra 1.2  Value 	 Keyword         Definition
 20             0       OFFDISK         The pixel is located somewhere off the solar disk.
 21             1       QUIET           The pixel is associated with weak line-of-sight magnetic field.
 22             2       ACTIVE          The pixel is associated with strong line-of-sight magnetic field.
 23             4       SPTQUIET        The pixel is associated with no features in the continuum intensity data.
 24             8       SPTFACUL        The pixel is associated with faculae in the continuum intensity data.
 25             16      SPTSPOT         The pixel is associated with sunspots in the continuum intensity data.
 26             32      ON_PATCH        The pixel is inside the patch.
 27            
 28             Here are all the possible permutations:
 29            
 30             Value   Definition
 31             0       The pixel is located somewhere off the solar disk.
 32             5       The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data.
 33             9       The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data.
 34             17      The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and sunspots in the continuum intensity data.
 35             6       The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and no features in the continuum intensity data.
 36             10      The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data.
 37             18      The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data.
 38             37      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data.
 39             41      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data.
 40 mbobra 1.2  49      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and sunspots in the continuum intensity data.
 41             38      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and no features in the continuum intensity data.
 42             42      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data.
 43             50      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data.
 44            
 45             Thus pixels with a value greater than 36 are located inside the patch.
 46            
 47             Written by Monica Bobra

 48 mbobra 1.1  
 49             ===========================================*/
 50            #include <math.h>
 51            #include <mkl.h>
 52            
 53            #define PI       (M_PI)
 54            #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
 55            
 56            /*===========================================*/
 57            
 58            /* Example function 1: Compute total unsigned flux in units of G/cm^2 */
 59            
 60            //  To compute the unsigned flux, we simply calculate
 61            //  flux = surface integral [(vector Bz) dot (normal vector)],
 62            //       = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
 63            //  However, since the field is radial, we will assume cos theta = 1.
 64            //  Therefore the pixels only need to be corrected for the projection.
 65            
 66            //  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
 67            //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
 68            //  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
 69 mbobra 1.1 //  =Gauss*cm^2
 70            
 71            int computeAbsFlux(float *bz, int *dims, float *absFlux,
 72                               float *mean_vf_ptr, float *count_mask_ptr,
 73                               int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
 74            
 75            {
 76                
 77                int nx = dims[0];
 78                int ny = dims[1];
 79                int i = 0;
 80                int j = 0;
 81                int count_mask = 0;
 82                double sum = 0.0;
 83                *absFlux = 0.0;
 84                *mean_vf_ptr = 0.0;
 85                
 86                 
 87                if (nx <= 0 || ny <= 0) return 1;
 88                
 89            	for (i = 0; i < nx; i++)
 90 mbobra 1.1 	{
 91            	   for (j = 0; j < ny; j++)
 92            	   {

 93 mbobra 1.2 	    if ( bitmask[j * nx + i] < 36 ) continue;

 94 mbobra 1.1             if isnan(bz[j * nx + i]) continue;
 95                        sum += (fabs(bz[j * nx + i]));
 96                        count_mask++;
 97            	   }
 98            	}
 99                
100                *mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
101                *count_mask_ptr  = count_mask;
102                //printf("mean_vf_ptr=%f\n",*mean_vf_ptr);
103                //printf("count_mask_ptr=%f\n",*count_mask_ptr);
104                //printf("cdelt1=%f\n",cdelt1);
105                //printf("rsun_ref=%f\n",rsun_ref);
106                //printf("rsun_obs=%f\n",rsun_obs);
107            
108                return 0;
109            }
110            
111            /*===========================================*/
112            /* Example function 2:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
113            
114            int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *bitmask, float *derx_bz, float *dery_bz)
115 mbobra 1.1 {
116                
117                int nx = dims[0];
118                int ny = dims[1];
119                int i = 0;
120                int j = 0;
121                int count_mask = 0;
122                double sum = 0.0;
123                *mean_derivative_bz_ptr = 0.0;
124                
125                if (nx <= 0 || ny <= 0) return 1;
126                
127                /* brute force method of calculating the derivative (no consideration for edges) */
128                for (i = 1; i <= nx-2; i++)
129                {
130            	for (j = 0; j <= ny-1; j++)
131                    {
132                       derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
133                    }
134                }
135                
136 mbobra 1.1     /* brute force method of calculating the derivative (no consideration for edges) */
137                for (i = 0; i <= nx-1; i++)
138                {
139                    for (j = 1; j <= ny-2; j++)
140                    {
141                       dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
142                    }
143                }
144                
145                /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
146                ignore the edges for the error terms as those arrays have been initialized to zero. 
147                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.*/
148                i=0;
149                for (j = 0; j <= ny-1; j++)
150                {
151                    derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
152                }
153                
154                i=nx-1;
155                for (j = 0; j <= ny-1; j++)
156                {
157 mbobra 1.1         derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
158                }
159                
160                j=0;
161                for (i = 0; i <= nx-1; i++)
162                {
163                    dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
164                }
165                
166                j=ny-1;
167                for (i = 0; i <= nx-1; i++)
168                {
169                    dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
170                }
171                
172                
173                for (i = 0; i <= nx-1; i++)
174                {
175                    for (j = 0; j <= ny-1; j++)
176                    {

177 mbobra 1.2             if ( bitmask[j * nx + i] < 36 ) continue;

178 mbobra 1.1             if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
179                        if isnan(bz[j * nx + i])      continue;
180                        if isnan(bz[(j+1) * nx + i])  continue;
181                        if isnan(bz[(j-1) * nx + i])  continue;
182                        if isnan(bz[j * nx + i-1])    continue;
183                        if isnan(bz[j * nx + i+1])    continue;
184                        if isnan(derx_bz[j * nx + i]) continue;
185                        if isnan(dery_bz[j * nx + i]) continue;
186                        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 */
187                        count_mask++;
188                    }
189                }
190                
191                *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
192                //printf("mean_derivative_bz_ptr=%f\n",*mean_derivative_bz_ptr);
193                
194            	return 0;
195            }
196            
197            /*===========================================*/
198            /* Example function 3: R parameter as defined in Schrijver, 2007 */
199 mbobra 1.1 //
200            // Note that there is a restriction on the function fsample()
201            // If the following occurs:
202            //      nx_out > floor((ny_in-1)/scale + 1) 
203            //      ny_out > floor((ny_in-1)/scale + 1),
204            // where n*_out are the dimensions of the output array and n*_in 
205            // are the dimensions of the input array, fsample() will usually result 
206            // in a segfault (though not always, depending on how the segfault was accessed.) 
207            
208            int computeR(float *los, int *dims, float *Rparam, float cdelt1,
209                         float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
210                         float *pmap, int nx1, int ny1, 
211                         int scale, float *p1pad, int nxp, int nyp, float *pmapn)
212            
213            { 
214                int nx = dims[0];
215                int ny = dims[1];
216                int i = 0;
217                int j = 0;
218                int index, index1;
219                double sum = 0.0;
220 mbobra 1.1     *Rparam = 0.0;
221                struct fresize_struct fresboxcar, fresgauss;
222                struct fint_struct fints;
223                float sigma = 10.0/2.3548;
224                
225                // set up convolution kernels
226                init_fresize_boxcar(&fresboxcar,1,1);
227                init_fresize_gaussian(&fresgauss,sigma,20,1);
228            
229                // =============== [STEP 1] =============== 
230                // bin the line-of-sight magnetogram down by a factor of scale 
231                fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
232            
233                // =============== [STEP 2] =============== 
234                // identify positive and negative pixels greater than +/- 150 gauss
235                // and label those pixels with a 1.0 in arrays p1p0 and p1n0
236                for (i = 0; i < nx1; i++)
237                {
238                    for (j = 0; j < ny1; j++)
239                    {
240                        index = j * nx1 + i;
241 mbobra 1.1             if (rim[index] > 150)
242                            p1p0[index]=1.0;
243                        else
244                            p1p0[index]=0.0;
245                        if (rim[index] < -150)
246                            p1n0[index]=1.0;
247                        else
248                            p1n0[index]=0.0;
249                    }
250                }
251            
252                // =============== [STEP 3] =============== 
253                // smooth each of the negative and positive pixel bitmaps      
254                fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
255                fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
256            
257                // =============== [STEP 4] =============== 
258                // find the pixels for which p1p and p1n are both equal to 1. 
259                // this defines the polarity inversion line
260                for (i = 0; i < nx1; i++)
261                {
262 mbobra 1.1         for (j = 0; j < ny1; j++)
263                    {
264                        index = j * nx1 + i;
265                        if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
266                            p1[index]=1.0;
267                        else
268                            p1[index]=0.0;
269                    }
270                }
271            
272                // pad p1 with zeroes so that the gaussian colvolution in step 5
273                // does not cut off data within hwidth of the edge
274               
275                // step i: zero p1pad
276                for (i = 0; i < nxp; i++)
277                {
278                    for (j = 0; j < nyp; j++)
279                    {
280                        index = j * nxp + i;
281                        p1pad[index]=0.0;
282                    }
283 mbobra 1.1     }
284            
285                // step ii: place p1 at the center of p1pad
286                for (i = 0; i < nx1; i++)
287                {
288                   for (j = 0; j < ny1; j++)
289                   {
290                        index  = j * nx1 + i; 
291                        index1 = (j+20) * nxp + (i+20);
292                        p1pad[index1]=p1[index];
293                    }
294                }
295            
296                // =============== [STEP 5] =============== 
297                // convolve the polarity inversion line map with a gaussian
298                // to identify the region near the plarity inversion line
299                // the resultant array is called pmap
300                fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0);
301            
302            
303               // select out the nx1 x ny1 non-padded array  within pmap
304 mbobra 1.1     for (i = 0; i < nx1; i++)
305                {
306                   for (j = 0; j < ny1; j++)
307                   {
308                        index  = j * nx1 + i; 
309                        index1 = (j+20) * nxp + (i+20);
310                        pmapn[index]=pmap[index1];
311                    }
312                }
313            
314                // =============== [STEP 6] =============== 
315                // the R parameter is calculated
316                for (i = 0; i < nx1; i++)
317                {
318                    for (j = 0; j < ny1; j++)
319                    {
320                        index = j * nx1 + i;
321                        if isnan(pmapn[index]) continue;
322                        if isnan(rim[index]) continue;
323                        sum += pmapn[index]*abs(rim[index]);
324                    }
325 mbobra 1.1     }
326                
327                if (sum < 1.0)
328                    *Rparam = 0.0;
329                else
330                    *Rparam = log10(sum);
331            
332                //printf("R_VALUE=%f\n",*Rparam);
333            
334                free_fresize(&fresboxcar);
335                free_fresize(&fresgauss);
336            
337                return 0;
338            
339            }
340            
341            /*===========================================*/
342            
343            char *smarp_functions_version() // Returns CVS version of smarp_functions.c
344            {
345                return strdup("$Id");
346 mbobra 1.1 }