sharp/apps/smarp_functions.c


/*===========================================
 
 The following 3 functions calculate the following spaceweather indices:
 
 USFLUX Total unsigned flux in Maxwells
 MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
 R_VALUE Karel Schrijver's R parameter
 
 The indices are calculated on the pixels in which the bitmap segment is greater than 36. 
 
 The SMARP bitmap has 13 unique values because they describe three different characteristics: 
 the location of the pixel (whether the pixel is off disk or a member of the patch), the 
 character of the magnetic field (quiet or active), and the character of the continuum 
 intensity data (quiet, part of faculae, part of a sunspot). 

 Here are all the possible values:
 
 Value   Keyword         Definition
 0       OFFDISK         The pixel is located somewhere off the solar disk.
 1       QUIET           The pixel is associated with weak line-of-sight magnetic field.
 2       ACTIVE          The pixel is associated with strong line-of-sight magnetic field.
 4       SPTQUIET        The pixel is associated with no features in the continuum intensity data.
 8       SPTFACUL        The pixel is associated with faculae in the continuum intensity data.
 16      SPTSPOT         The pixel is associated with sunspots in the continuum intensity data.
 32      ON_PATCH        The pixel is inside the patch.

 Here are all the possible permutations:

 Value   Definition
 0       The pixel is located somewhere off the solar disk.
 5       The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data.
 9       The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data.
 17      The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and sunspots in the continuum intensity data.
 6       The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and no features in the continuum intensity data.
 10      The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data.
 18      The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data.
 37      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data.
 41      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data.
 49      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and sunspots in the continuum intensity data.
 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      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data.
 50      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data.

 Thus pixels with a value greater than 36 are located inside the patch.

 Written by Monica Bobra
 
 ===========================================*/
#include <math.h>
#include <mkl.h>

#define PI       (M_PI)
#define MUNAUGHT (0.0000012566370614) /* magnetic constant */

/*===========================================*/

/* Example function 1: Compute total unsigned flux in units of G/cm^2 */

//  To compute the unsigned flux, we simply calculate
//  flux = surface integral [(vector Bz) dot (normal vector)],
//       = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
//  However, since the field is radial, we will assume cos theta = 1.
//  Therefore the pixels only need to be corrected for the projection.

//  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
//  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
//  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
//  =Gauss*cm^2

int computeAbsFlux(float *bz, int *dims, float *absFlux,
                   float *mean_vf_ptr, float *count_mask_ptr,
                   int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

{
    
    int nx = dims[0];
    int ny = dims[1];
    int i = 0;
    int j = 0;
    int count_mask = 0;
    double sum = 0.0;
    *absFlux = 0.0;
    *mean_vf_ptr = 0.0;
    
     
    if (nx <= 0 || ny <= 0) return 1;
    
        for (i = 0; i < nx; i++)
        {
           for (j = 0; j < ny; j++)
           {
            if ( bitmask[j * nx + i] < 36 ) continue;
            if isnan(bz[j * nx + i]) continue;
            sum += (fabs(bz[j * nx + i]));
            count_mask++;
           }
        }
    
    *mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
    *count_mask_ptr  = count_mask;
    //printf("mean_vf_ptr=%f\n",*mean_vf_ptr);
    //printf("count_mask_ptr=%f\n",*count_mask_ptr);
    //printf("cdelt1=%f\n",cdelt1);
    //printf("rsun_ref=%f\n",rsun_ref);
    //printf("rsun_obs=%f\n",rsun_obs);

    return 0;
}

/*===========================================*/
/* Example function 2:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */

int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *bitmask, float *derx_bz, float *dery_bz)
{
    
    int nx = dims[0];
    int ny = dims[1];
    int i = 0;
    int j = 0;
    int count_mask = 0;
    double sum = 0.0;
    *mean_derivative_bz_ptr = 0.0;
    
    if (nx <= 0 || ny <= 0) return 1;
    
    /* brute force method of calculating the derivative (no consideration for edges) */
    for (i = 1; i <= nx-2; i++)
    {
        for (j = 0; j <= ny-1; j++)
        {
           derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
        }
    }
    
    /* brute force method of calculating the derivative (no consideration for edges) */
    for (i = 0; i <= nx-1; i++)
    {
        for (j = 1; j <= ny-2; j++)
        {
           dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
        }
    }
    
    /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
    ignore the edges for the error terms as those arrays have been initialized to zero. 
    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.*/
    i=0;
    for (j = 0; j <= ny-1; j++)
    {
        derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
    }
    
    i=nx-1;
    for (j = 0; j <= ny-1; j++)
    {
        derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
    }
    
    j=0;
    for (i = 0; i <= nx-1; i++)
    {
        dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
    }
    
    j=ny-1;
    for (i = 0; i <= nx-1; i++)
    {
        dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
    }
    
    
    for (i = 0; i <= nx-1; i++)
    {
        for (j = 0; j <= ny-1; j++)
        {
            if ( bitmask[j * nx + i] < 36 ) continue;
            if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
            if isnan(bz[j * nx + i])      continue;
            if isnan(bz[(j+1) * nx + i])  continue;
            if isnan(bz[(j-1) * nx + i])  continue;
            if isnan(bz[j * nx + i-1])    continue;
            if isnan(bz[j * nx + i+1])    continue;
            if isnan(derx_bz[j * nx + i]) continue;
            if isnan(dery_bz[j * nx + i]) continue;
            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 */
            count_mask++;
        }
    }
    
    *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
    //printf("mean_derivative_bz_ptr=%f\n",*mean_derivative_bz_ptr);
    
        return 0;
}

/*===========================================*/
/* Example function 3: R parameter as defined in Schrijver, 2007 */
//
// Note that there is a restriction on the function fsample()
// If the following occurs:
//      nx_out > floor((ny_in-1)/scale + 1) 
//      ny_out > floor((ny_in-1)/scale + 1),
// where n*_out are the dimensions of the output array and n*_in 
// are the dimensions of the input array, fsample() will usually result 
// in a segfault (though not always, depending on how the segfault was accessed.) 

int computeR(float *los, int *dims, float *Rparam, float cdelt1,
             float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
             float *pmap, int nx1, int ny1, 
             int scale, float *p1pad, int nxp, int nyp, float *pmapn)

{ 
    int nx = dims[0];
    int ny = dims[1];
    int i = 0;
    int j = 0;
    int index, index1;
    double sum = 0.0;
    *Rparam = 0.0;
    struct fresize_struct fresboxcar, fresgauss;
    struct fint_struct fints;
    float sigma = 10.0/2.3548;
    
    // set up convolution kernels
    init_fresize_boxcar(&fresboxcar,1,1);
    init_fresize_gaussian(&fresgauss,sigma,20,1);

    // =============== [STEP 1] =============== 
    // bin the line-of-sight magnetogram down by a factor of scale 
    fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);

    // =============== [STEP 2] =============== 
    // identify positive and negative pixels greater than +/- 150 gauss
    // and label those pixels with a 1.0 in arrays p1p0 and p1n0
    for (i = 0; i < nx1; i++)
    {
        for (j = 0; j < ny1; j++)
        {
            index = j * nx1 + i;
            if (rim[index] > 150)
                p1p0[index]=1.0;
            else
                p1p0[index]=0.0;
            if (rim[index] < -150)
                p1n0[index]=1.0;
            else
                p1n0[index]=0.0;
        }
    }

    // =============== [STEP 3] =============== 
    // smooth each of the negative and positive pixel bitmaps      
    fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
    fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);

    // =============== [STEP 4] =============== 
    // find the pixels for which p1p and p1n are both equal to 1. 
    // this defines the polarity inversion line
    for (i = 0; i < nx1; i++)
    {
        for (j = 0; j < ny1; j++)
        {
            index = j * nx1 + i;
            if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
                p1[index]=1.0;
            else
                p1[index]=0.0;
        }
    }

    // pad p1 with zeroes so that the gaussian colvolution in step 5
    // does not cut off data within hwidth of the edge
   
    // step i: zero p1pad
    for (i = 0; i < nxp; i++)
    {
        for (j = 0; j < nyp; j++)
        {
            index = j * nxp + i;
            p1pad[index]=0.0;
        }
    }

    // step ii: place p1 at the center of p1pad
    for (i = 0; i < nx1; i++)
    {
       for (j = 0; j < ny1; j++)
       {
            index  = j * nx1 + i; 
            index1 = (j+20) * nxp + (i+20);
            p1pad[index1]=p1[index];
        }
    }

    // =============== [STEP 5] =============== 
    // convolve the polarity inversion line map with a gaussian
    // to identify the region near the plarity inversion line
    // the resultant array is called pmap
    fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0);


   // select out the nx1 x ny1 non-padded array  within pmap
    for (i = 0; i < nx1; i++)
    {
       for (j = 0; j < ny1; j++)
       {
            index  = j * nx1 + i; 
            index1 = (j+20) * nxp + (i+20);
            pmapn[index]=pmap[index1];
        }
    }

    // =============== [STEP 6] =============== 
    // the R parameter is calculated
    for (i = 0; i < nx1; i++)
    {
        for (j = 0; j < ny1; j++)
        {
            index = j * nx1 + i;
            if isnan(pmapn[index]) continue;
            if isnan(rim[index]) continue;
            sum += pmapn[index]*abs(rim[index]);
        }
    }
    
    if (sum < 1.0)
        *Rparam = 0.0;
    else
        *Rparam = log10(sum);

    //printf("R_VALUE=%f\n",*Rparam);

    free_fresize(&fresboxcar);
    free_fresize(&fresgauss);

    return 0;

}

/*===========================================*/

char *smarp_functions_version() // Returns CVS version of smarp_functions.c
{
    return strdup("$Id");
}
Revision:	1.2
Committed:	Mon May 4 21:11:07 2020 UTC (3 years, 4 months ago) by mbobra
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+38 -18 lines
Log Message:	Changed bitmap threshold to pixels > 36 (inside the patch) and added documentation.
#	User	Rev	Content
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			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			// =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			{
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 = sumcdelt1cdelt1(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			{
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			/* 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] = ( (-3bz[j nx + i]) + (4bz[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			derx_bz[j * nx + i] = ( (3bz[j nx + i]) + (-4bz[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] = ( (-3bz[jnx + i]) + (4bz[(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] = ( (3bz[j nx + i]) + (-4bz[(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			//
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			*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			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			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			}
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			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			}
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			}