Return to smarp_functions.c CVS log | Up to [Development] / JSOC / proj / sharp / apps |
File: [Development] / JSOC / proj / sharp / apps / smarp_functions.c
(download)
Revision: 1.4, Tue May 25 00:15:39 2021 UTC (11 months, 4 weeks ago) by mbobra Branch: MAIN CVS Tags: HEAD Changes since 1.3: +11 -11 lines removed check for ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) |
/*=========================================== The following 3 functions calculate the following spaceweather indices on line-of-sight magnetic field data: USFLUX Total unsigned flux in Maxwells MEANGBZ Mean value of the line-of-sight (approximately vertical in remapped and reprojected data) field gradient, in Gauss/Mm R_VALUE R parameter (Schrijver, 2007) 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, faculae, 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 LOS) dot (normal vector)], // = surface integral [(magnitude LOS)*(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_los(float *los, 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(los[j * nx + i]) continue; sum += (fabs(los[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_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */ int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los) { int nx = dims[0]; int ny = dims[1]; int i = 0; int j = 0; int count_mask = 0; double sum = 0.0; *mean_derivative_los_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_los[j * nx + i] = (los[j * nx + i+1] - los[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_los[j * nx + i] = (los[(j+1) * nx + i] - los[(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_los[j * nx + i] = ( (-3*los[j * nx + i]) + (4*los[j * nx + (i+1)]) - (los[j * nx + (i+2)]) )*0.5; } i=nx-1; for (j = 0; j <= ny-1; j++) { derx_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[j * nx + (i-1)]) - (-los[j * nx + (i-2)]) )*0.5; } j=0; for (i = 0; i <= nx-1; i++) { dery_los[j * nx + i] = ( (-3*los[j*nx + i]) + (4*los[(j+1) * nx + i]) - (los[(j+2) * nx + i]) )*0.5; } j=ny-1; for (i = 0; i <= nx-1; i++) { dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5; } for (i = 1; i <= nx-2; i++) { for (j = 1; j <= ny-2; j++) { if ( bitmask[j * nx + i] < 36 ) continue; if isnan(los[j * nx + i]) continue; if isnan(los[(j+1) * nx + i]) continue; if isnan(los[(j-1) * nx + i]) continue; if isnan(los[j * nx + i-1]) continue; if isnan(los[j * nx + i+1]) continue; if isnan(derx_los[j * nx + i]) continue; if isnan(dery_los[j * nx + i]) continue; 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 */ count_mask++; } } *mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram //printf("mean_derivative_los_ptr=%f\n",*mean_derivative_los_ptr); //printf("nx=%d\n",nx); //printf("ny=%d\n",ny); //printf("sum=%f\n",sum); 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_los(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"); }
Karen Tian |
Powered by ViewCVS 0.9.4 |