Return to smarp_functions.c CVS log Up to [Development] / JSOC / proj / sharp / apps

 ``` 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 31 #include 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 } ```