JSOC/proj/sharp/apps/smarp_functions.c - annotate

Return to smarp_functions.c CVS log

Up to [Development] / JSOC / proj / sharp / apps

1 mbobra 1.1 2 /*=========================================== 3
4 mbobra 1.3 The following 3 functions calculate the following spaceweather indices on line-of-sight 5 magnetic field data:
6 mbobra 1.1
7 mbobra 1.3 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 mbobra 1.1
11 mbobra 1.2 The indices are calculated on the pixels in which the bitmap segment is greater than 36.
12 mbobra 1.1
13 mbobra 1.2 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 mbobra 1.3 intensity data (quiet, faculae, sunspot).
17 mbobra 1.2 18 Here are all the possible values:
19 mbobra 1.1
20 mbobra 1.2 Value Keyword Definition 21 0 OFFDISK The pixel is located somewhere off the solar disk. 22 1 QUIET The pixel is associated with weak line-of-sight magnetic field. 23 2 ACTIVE The pixel is associated with strong line-of-sight magnetic field. 24 4 SPTQUIET The pixel is associated with no features in the continuum intensity data. 25 8 SPTFACUL The pixel is associated with faculae in the continuum intensity data. 26 16 SPTSPOT The pixel is associated with sunspots in the continuum intensity data. 27 32 ON_PATCH The pixel is inside the patch. 28 29 Here are all the possible permutations: 30 31 Value Definition 32 0 The pixel is located somewhere off the solar disk. 33 5 The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data. 34 9 The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data. 35 17 The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and sunspots in the continuum intensity data. 36 6 The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and no features in the continuum intensity data. 37 10 The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data. 38 18 The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data. 39 37 The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data. 40 41 The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data. 41 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. 42 38 The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and no features in the continuum intensity data. 43 42 The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data. 44 50 The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data. 45 46 Thus pixels with a value greater than 36 are located inside the patch. 47 48 Written by Monica Bobra
49 mbobra 1.1 50 ===========================================/ 51 #include <math.h> 52 #include <mkl.h> 53 54 #define PI (M_PI) 55 #define MUNAUGHT (0.0000012566370614) / magnetic constant / 56 57 /===========================================/ 58 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 mbobra 1.3 // flux = surface integral [(vector LOS) dot (normal vector)], 63 // = surface integral [(magnitude LOS)(magnitude normal)(cos theta)].
64 mbobra 1.1 // 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 67 // To convert G to Gcm^2, simply multiply by the number of square centimeters per pixel. 68 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 69 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2 70 // =Gausscm^2 71
72 mbobra 1.3 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 mbobra 1.1 76 { 77 78 int nx = dims[0]; 79 int ny = dims[1]; 80 int i = 0; 81 int j = 0; 82 int count_mask = 0; 83 double sum = 0.0; 84 absFlux = 0.0; 85 mean_vf_ptr = 0.0; 86 87 88 if (nx <= 0 \|\| ny <= 0) return 1; 89 90 for (i = 0; i < nx; i++) 91 { 92 for (j = 0; j < ny; j++) 93 {
94 mbobra 1.2 if ( bitmask[j * nx + i] < 36 ) continue;
95 mbobra 1.3 if isnan(los[j * nx + i]) continue; 96 sum += (fabs(los[j * nx + i]));
97 mbobra 1.1 count_mask++; 98 } 99 } 100 101 mean_vf_ptr = sumcdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 102 count_mask_ptr = count_mask; 103 //printf("mean_vf_ptr=%f\n",mean_vf_ptr); 104 //printf("count_mask_ptr=%f\n",count_mask_ptr); 105 //printf("cdelt1=%f\n",cdelt1); 106 //printf("rsun_ref=%f\n",rsun_ref); 107 //printf("rsun_obs=%f\n",rsun_obs); 108 109 return 0; 110 } 111 112 /===========================================/
113 mbobra 1.3 /* Example function 2: Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */
114 mbobra 1.1
115 mbobra 1.3 int computeLOSderivative(float los, int dims, float mean_derivative_los_ptr, int bitmask, float derx_los, float dery_los)
116 mbobra 1.1 { 117 118 int nx = dims[0]; 119 int ny = dims[1]; 120 int i = 0; 121 int j = 0; 122 int count_mask = 0; 123 double sum = 0.0;
124 mbobra 1.3 *mean_derivative_los_ptr = 0.0;
125 mbobra 1.1 126 if (nx <= 0 \|\| ny <= 0) return 1; 127 128 /* brute force method of calculating the derivative (no consideration for edges) */ 129 for (i = 1; i <= nx-2; i++) 130 { 131 for (j = 0; j <= ny-1; j++) 132 {
133 mbobra 1.3 derx_los[j * nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])*0.5;
134 mbobra 1.1 } 135 } 136 137 /* brute force method of calculating the derivative (no consideration for edges) */ 138 for (i = 0; i <= nx-1; i++) 139 { 140 for (j = 1; j <= ny-2; j++) 141 {
142 mbobra 1.3 dery_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
143 mbobra 1.1 } 144 } 145 146 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below. 147 ignore the edges for the error terms as those arrays have been initialized to zero. 148 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.*/ 149 i=0; 150 for (j = 0; j <= ny-1; j++) 151 {
152 mbobra 1.3 derx_los[j * nx + i] = ( (-3los[j nx + i]) + (4los[j nx + (i+1)]) - (los[j * nx + (i+2)]) )*0.5;
153 mbobra 1.1 } 154 155 i=nx-1; 156 for (j = 0; j <= ny-1; j++) 157 {
158 mbobra 1.3 derx_los[j * nx + i] = ( (3los[j nx + i]) + (-4los[j nx + (i-1)]) - (-los[j * nx + (i-2)]) )*0.5;
159 mbobra 1.1 } 160 161 j=0; 162 for (i = 0; i <= nx-1; i++) 163 {
164 mbobra 1.3 dery_los[j * nx + i] = ( (-3los[jnx + i]) + (4los[(j+1) nx + i]) - (los[(j+2) * nx + i]) )*0.5;
165 mbobra 1.1 } 166 167 j=ny-1; 168 for (i = 0; i <= nx-1; i++) 169 {
170 mbobra 1.3 dery_los[j * nx + i] = ( (3los[j nx + i]) + (-4los[(j-1) nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
171 mbobra 1.1 } 172
173 mbobra 1.4 for (i = 1; i <= nx-2; i++)
174 mbobra 1.1 {
175 mbobra 1.4 for (j = 1; j <= ny-2; j++)
176 mbobra 1.1 {
177 mbobra 1.2 if ( bitmask[j * nx + i] < 36 ) continue;
178 mbobra 1.4 if isnan(los[j * nx + i]) continue; 179 if isnan(los[(j+1) * nx + i]) continue; 180 if isnan(los[(j-1) * nx + i]) continue; 181 if isnan(los[j * nx + i-1]) continue; 182 if isnan(los[j * nx + i+1]) continue;
183 mbobra 1.3 if isnan(derx_los[j * nx + i]) continue; 184 if isnan(dery_los[j * nx + i]) continue; 185 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 */
186 mbobra 1.1 count_mask++; 187 } 188 } 189
190 mbobra 1.3 mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)(ny-2)) if shape of count_mask = shape of magnetogram 191 //printf("mean_derivative_los_ptr=%f\n",*mean_derivative_los_ptr);
192 mbobra 1.4 //printf("nx=%d\n",nx); 193 //printf("ny=%d\n",ny); 194 //printf("sum=%f\n",sum); 195 return 0;
196 mbobra 1.1 } 197 198 /===========================================/ 199 /* Example function 3: R parameter as defined in Schrijver, 2007 / 200 // 201 // Note that there is a restriction on the function fsample() 202 // If the following occurs: 203 // nx_out > floor((ny_in-1)/scale + 1) 204 // ny_out > floor((ny_in-1)/scale + 1), 205 // where n_out are the dimensions of the output array and n*_in 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 mbobra 1.3 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 mbobra 1.1 214 { 215 int nx = dims[0]; 216 int ny = dims[1]; 217 int i = 0; 218 int j = 0; 219 int index, index1; 220 double sum = 0.0; 221 Rparam = 0.0; 222 struct fresize_struct fresboxcar, fresgauss; 223 struct fint_struct fints; 224 float sigma = 10.0/2.3548; 225 226 // set up convolution kernels 227 init_fresize_boxcar(&fresboxcar,1,1); 228 init_fresize_gaussian(&fresgauss,sigma,20,1); 229 230 // =============== [STEP 1] =============== 231 // bin the line-of-sight magnetogram down by a factor of scale 232 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0); 233 234 mbobra 1.1 // =============== [STEP 2] =============== 235 // identify positive and negative pixels greater than +/- 150 gauss 236 // and label those pixels with a 1.0 in arrays p1p0 and p1n0 237 for (i = 0; i < nx1; i++) 238 { 239 for (j = 0; j < ny1; j++) 240 { 241 index = j nx1 + i; 242 if (rim[index] > 150) 243 p1p0[index]=1.0; 244 else 245 p1p0[index]=0.0; 246 if (rim[index] < -150) 247 p1n0[index]=1.0; 248 else 249 p1n0[index]=0.0; 250 } 251 } 252 253 // =============== [STEP 3] =============== 254 // smooth each of the negative and positive pixel bitmaps 255 mbobra 1.1 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0); 256 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0); 257 258 // =============== [STEP 4] =============== 259 // find the pixels for which p1p and p1n are both equal to 1. 260 // this defines the polarity inversion line 261 for (i = 0; i < nx1; i++) 262 { 263 for (j = 0; j < ny1; j++) 264 { 265 index = j * nx1 + i; 266 if ((p1p[index] > 0.0) && (p1n[index] > 0.0)) 267 p1[index]=1.0; 268 else 269 p1[index]=0.0; 270 } 271 } 272 273 // pad p1 with zeroes so that the gaussian colvolution in step 5 274 // does not cut off data within hwidth of the edge 275 276 mbobra 1.1 // step i: zero p1pad 277 for (i = 0; i < nxp; i++) 278 { 279 for (j = 0; j < nyp; j++) 280 { 281 index = j * nxp + i; 282 p1pad[index]=0.0; 283 } 284 } 285 286 // step ii: place p1 at the center of p1pad 287 for (i = 0; i < nx1; i++) 288 { 289 for (j = 0; j < ny1; j++) 290 { 291 index = j * nx1 + i; 292 index1 = (j+20) * nxp + (i+20); 293 p1pad[index1]=p1[index]; 294 } 295 } 296 297 mbobra 1.1 // =============== [STEP 5] =============== 298 // convolve the polarity inversion line map with a gaussian 299 // to identify the region near the plarity inversion line 300 // the resultant array is called pmap 301 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0); 302 303 304 // select out the nx1 x ny1 non-padded array within pmap 305 for (i = 0; i < nx1; i++) 306 { 307 for (j = 0; j < ny1; j++) 308 { 309 index = j * nx1 + i; 310 index1 = (j+20) * nxp + (i+20); 311 pmapn[index]=pmap[index1]; 312 } 313 } 314 315 // =============== [STEP 6] =============== 316 // the R parameter is calculated 317 for (i = 0; i < nx1; i++) 318 mbobra 1.1 { 319 for (j = 0; j < ny1; j++) 320 { 321 index = j * nx1 + i; 322 if isnan(pmapn[index]) continue; 323 if isnan(rim[index]) continue; 324 sum += pmapn[index]abs(rim[index]); 325 } 326 } 327 328 if (sum < 1.0) 329 Rparam = 0.0; 330 else 331 Rparam = log10(sum); 332 333 //printf("R_VALUE=%f\n",Rparam); 334 335 free_fresize(&fresboxcar); 336 free_fresize(&fresgauss); 337 338 return 0; 339 mbobra 1.1 340 } 341 342 /===========================================/ 343 344 char *smarp_functions_version() // Returns CVS version of smarp_functions.c 345 { 346 return strdup("$Id"); 347 }

Karen Tian