(file) Return to update_sharp_keys.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

  1 mbobra 1.1 /*                                                                      
  2             *      MODULE NAME: update_sharp_keys.c
  3             *
  4             *	DESCRIPTION: This module recalculates SHARP keywords. This is accomplished by 
  5             *      cloning a record, recalculating keywords of choice (i.e., user input), 
  6             *      and pointing to the same segments as the old record. Associated error keys
  7             *      are computed by default. 
  8             * 
  9             *      This module accounts for versioning by appending to the keyword CODEVER7 and HISTORY:
 10             *      CODEVER7 will contain multiple lines: the production build of sharp.c, the 
 11             *      production build of include file sw_functions.c, and the production build
 12             *      of update_sharp_keys.c. HISTORY will include a human-readable sentence about which
 13             *      keywords were updated. 
 14             * 
 15             *      This module does not produce segments.
 16             *
 17             *      INPUTS     : -- DRMS SHARP series
 18             *                   -- DRMS SHARP CEA series
 19             *                   -- HARPNUM
 20             *                   -- comma separated list of keywords to recalculate
 21             *
 22 mbobra 1.1  *	AUTHOR     : Monica Bobra
 23             *
 24             *	Version    :   v0.0	Jun 14 2013
 25 mbobra 1.6  *                     v0.1     Feb 11 2014 Added different input and output series
 26 mbobra 1.1  *
 27             *	EXAMPLE    :
 28 mbobra 1.6  *      update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //
 29             *      HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT
 30 mbobra 1.1  *
 31             */
 32            
 33            /* Include files */
 34            #include <stdio.h>
 35            #include <stdlib.h>
 36            #include <time.h>
 37            #include <sys/time.h>
 38            #include <math.h>
 39            #include <string.h>
 40            #include "jsoc_main.h"
 41            #include "astro.h"
 42            #include "fstats.h"
 43            #include "cartography.c"
 44            #include "fresize.h"
 45            #include "finterpolate.h"
 46            #include "img2helioVector.c"
 47            #include "copy_me_keys.c"
 48            #include "errorprop.c"
 49            #include "sw_functions.c"
 50            #include <mkl_blas.h>
 51 mbobra 1.1 #include <mkl_service.h>
 52            #include <mkl_lapack.h>
 53            #include <mkl_vml_functions.h>
 54            #include <omp.h>
 55             
 56            /* Define values */                     
 57            #define PI              (M_PI)
 58            #define RADSINDEG       (PI/180.)
 59            #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
 60            #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
 61            #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
 62            
 63            char *module_name = "update_sharp_keys";  /* Module name */
 64 mbobra 1.6 char *version_id  = "2014 Feb 12";        /* Version number */
 65 mbobra 1.1 
 66            ModuleArgs_t module_args[] =
 67            {
 68 mbobra 1.6 	{ARG_STRING, "sharpseriesin",     NULL, "Input Sharp dataseries"},
 69            	{ARG_STRING, "sharpceaseriesin",  NULL, "Input Sharp CEA dataseries"},
 70            	{ARG_STRING, "sharpseriesout",     NULL, "Output Sharp dataseries"},
 71            	{ARG_STRING, "sharpceaseriesout",  NULL, "Output Sharp CEA dataseries"},
 72 mbobra 1.1 	{ARG_INT,    "HARPNUM",         NULL, "HARP number"},
 73            	{ARG_STRING, "keylist",         NULL, "comma separated list of keywords to update"},
 74            	{ARG_END}
 75            };
 76            
 77            
 78            int DoIt(void)                    
 79            {
 80                    int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
 81                    int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
 82            
 83            	int status = DRMS_SUCCESS;
 84            	int nrecs, irec;
 85            
 86                    /* Keywords */
 87            	float mean_vf; 
 88            	float absFlux;
 89                    float count_mask;
 90            	float mean_hf;
 91            	float mean_gamma;
 92            	float mean_derivative_btotal;
 93 mbobra 1.1 	float mean_derivative_bh; 
 94            	float mean_derivative_bz;
 95            	float mean_jz;
 96            	float us_i;
 97            	float mean_alpha;
 98            	float mean_ih;
 99            	float total_us_ih;
100            	float total_abs_ih;
101            	float totaljz;
102            	float totpot;
103            	float meanpot;
104            	float area_w_shear_gt_45;
105            	float meanshear_angle;
106            	float area_w_shear_gt_45h;
107            	float meanshear_angleh; 
108                    float mean_derivative_btotal_err;
109                    float mean_vf_err;
110                    float mean_gamma_err;
111                    float mean_derivative_bh_err;
112                    float mean_derivative_bz_err;
113                    float mean_jz_err;
114 mbobra 1.1         float us_i_err; 
115                    float mean_alpha_err;
116                    float mean_ih_err;
117                    float total_us_ih_err;
118                    float total_abs_ih_err;
119                    float totaljz_err;
120                    float meanpot_err;
121                    float totpot_err;
122 mbobra 1.6         float Rparam;
123 mbobra 1.1         float meanshear_angle_err;
124 mbobra 1.6 	char *sharpseriesin = (char *) params_get_str(&cmdparams, "sharpseriesin");
125            	char *sharpceaseriesin = (char *) params_get_str(&cmdparams, "sharpceaseriesin");
126            	char *sharpseriesout = (char *) params_get_str(&cmdparams, "sharpseriesout");
127            	char *sharpceaseriesout = (char *) params_get_str(&cmdparams, "sharpceaseriesout");
128            
129                    int sameflag = strcmp(sharpseriesin, sharpseriesout) == 0;
130                    int testflag = strcmp(sharpceaseriesin, sharpceaseriesout) == 0;
131                    if (testflag != sameflag)
132                      DIE("Either both outputs must be the same as their inputs, or both be different.");
133 mbobra 1.1 
134            	int harpnum = params_get_int(&cmdparams, "HARPNUM");
135 mbobra 1.6 
136 mbobra 1.1 	char sharpquery[256];
137            	char sharpceaquery[256];
138 mbobra 1.6 	printf("sameflag = %d\n", sameflag);
139            	printf("testflag = %d\n", testflag);
140            	sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
141            	printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
142 mbobra 1.1 
143 mbobra 1.6 	sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
144            	sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
145 mbobra 1.1 
146 mbobra 1.6 	DRMS_RecordSet_t *sharpinrecset= drms_open_records(drms_env, sharpquery, &status);
147 mbobra 1.1 	if (status != DRMS_SUCCESS)
148            		DIE("Problem opening sharp recordset.");
149 mbobra 1.6 	nrecs=sharpinrecset->n;
150            	if (nrecs == 0)
151 mbobra 1.1 		DIE("Sharp recordset contains no records.");
152 mbobra 1.6         printf("nrecs=%d\n",nrecs);
153 mbobra 1.1 
154            	DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status);
155            	if (status != DRMS_SUCCESS)
156            		DIE("Problem opening sharp cea recordset.");
157 mbobra 1.6 	if (sharpceainrecset->n != nrecs)
158            		DIE("Sharp cea recordset differs from sharp recordset.");
159            
160                    DRMS_RecordSet_t *sharpoutrecset, *sharpceaoutrecset;
161                    if (sameflag)
162                    {
163            	sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
164            	if (status != DRMS_SUCCESS)
165            		DIE("Problem cloning sharp records.");
166            	sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
167 mbobra 1.1 	if (status != DRMS_SUCCESS)
168            		DIE("Problem cloning sharp cea records.");
169 mbobra 1.6         }
170                    else
171                    {
172            	sharpoutrecset = drms_create_records(drms_env, nrecs, sharpseriesout, DRMS_PERMANENT, &status);
173            	if (status != DRMS_SUCCESS)
174            		DIE("Problem creating sharp records.");
175            	sharpceaoutrecset = drms_create_records(drms_env, nrecs, sharpceaseriesout, DRMS_PERMANENT, &status);
176            	if (status != DRMS_SUCCESS)
177            		DIE("Problem creating sharp cea records.");
178                    }
179            
180 mbobra 1.1 
181            	char *keylist = (char *) params_get_str(&cmdparams, "keylist");
182            
183                    // Flags to indicate which keyword will be recalculated
184            	int meanvfflag  = (strstr(keylist,"USFLUX")  != NULL);  // generalize so that lowercase is acceptable
185            	int totpotflag  = (strstr(keylist,"TOTPOT")  != NULL);  
186                    int meangamflag = (strstr(keylist,"MEANGAM") != NULL);
187                    int meangbtflag = (strstr(keylist,"MEANGBT") != NULL);
188                    int meangbzflag = (strstr(keylist,"MEANGBZ") != NULL);
189                    int meangbhflag = (strstr(keylist,"MEANGBH") != NULL);
190                    int meanjzdflag = (strstr(keylist,"MEANJZD") != NULL);
191                    int totusjzflag = (strstr(keylist,"TOTUSJZ") != NULL);
192                    int meanalpflag = (strstr(keylist,"MEANALP") != NULL);
193                    int meanjzhflag = (strstr(keylist,"MEANJZH") != NULL);
194                    int totusjhflag = (strstr(keylist,"TOTUSJH") != NULL);
195                    int absnjzhflag = (strstr(keylist,"ABSNJZH") != NULL);
196                    int savncppflag = (strstr(keylist,"SAVNCPP") != NULL);
197                    int meanpotflag = (strstr(keylist,"MEANPOT") != NULL);
198                    int meanshrflag = (strstr(keylist,"MEANSHR") != NULL);
199                    int shrgt45flag = (strstr(keylist,"SHRGT45") != NULL);
200 mbobra 1.6         int r_valueflag = (strstr(keylist,"R_VALUE") != NULL);
201 mbobra 1.1 	DRMS_Record_t *sharpinrec = sharpinrecset->records[0];
202            	DRMS_Record_t *sharpceainrec = sharpceainrecset->records[0];
203            	DRMS_Segment_t *inseg = drms_segment_lookup(sharpceainrec, "Br");
204            	int nx = inseg->axis[0];
205            	int ny = inseg->axis[1];
206            	int nxny = nx * ny;
207            	int dims[2] = {nx, ny};
208 mbobra 1.6  
209 mbobra 1.1 	// Temp arrays 	
210            	float *bh      = (float *) (malloc(nxny * sizeof(float)));
211            	float *bt      = (float *) (malloc(nxny * sizeof(float)));
212            	float *jz      = (float *) (malloc(nxny * sizeof(float)));
213            	float *bpx     = (float *) (malloc(nxny * sizeof(float)));
214            	float *bpy     = (float *) (malloc(nxny * sizeof(float)));
215            	float *bpz     = (float *) (malloc(nxny * sizeof(float)));
216            	float *derx    = (float *) (malloc(nxny * sizeof(float)));
217            	float *dery    = (float *) (malloc(nxny * sizeof(float)));
218            	float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
219            	float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
220            	float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
221            	float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
222            	float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
223            	float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
224            	float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
225            	float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
226                    float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
227                    float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
228                    float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
229                    float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
230 mbobra 1.1 	float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
231            
232                    // ephemeris variables
233            	float  cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2;
234            	double rsun_ref, rsun_obs;
235            
236                    // outrecords 
237 mbobra 1.6 	DRMS_Record_t *sharpoutrec, *sharpceaoutrec, *testrec;
238 mbobra 1.1 
239                    // prepare to set CODEVER7 (CVS Version of the SHARP module)
240            	char *cvsinfo0;
241            	char *history0;
242 mbobra 1.7 	char *cvsinfo1 = strdup("$Id: update_sharp_keys.c,v 1.6 2014/05/16 21:56:33 mbobra Exp $");
243 mbobra 1.1 	char *cvsinfo2 = sw_functions_version();
244            	char *cvsinfoall = (char *)malloc(2048);
245                    char historyofthemodule[2048];
246                    cvsinfo0    = drms_getkey_string(sharpinrec, "CODEVER7", &status);
247                    strcpy(cvsinfoall,cvsinfo0);
248                    strcat(cvsinfoall,"\n");
249                    strcat(cvsinfoall,cvsinfo2);
250            
251                    // prepare to set HISTORY (History of the data)
252            	char timebuf[1024];
253            	float UNIX_epoch = -220924792.000; 	
254            	sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
255                    sprintf(historyofthemodule,"The following keywords were re-computed on %s: %s.",timebuf,keylist); 
256            	printf("historyofthemodule=%s\n",historyofthemodule);
257            	history0    = drms_getkey_string(sharpinrec, "HISTORY", &status);
258            	strcat(historyofthemodule,"\n");
259                    strcat(historyofthemodule,history0);
260             
261 mbobra 1.6         // grab one record such that we can malloc the R-parameter calculation arrays outside of the loop
262                    // assume scale = round(2./cdelt1) does not change for any given HARP
263                    testrec  = sharpceainrecset->records[nrecs/2];
264                    dsun_obs = drms_getkey_float(testrec, "DSUN_OBS",   &status); 
265                    rsun_ref = drms_getkey_double(testrec, "RSUN_REF", &status);
266                    cdelt1_orig = drms_getkey_float(testrec, "CDELT1",   &status);
267                    cdelt1      = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
268                    int scale = round(2.0/cdelt1);
269                    int nx1 = nx/scale;
270                    int ny1 = ny/scale;
271 mbobra 1.7         int nxp = nx1+40;
272                    int nyp = ny1+40;
273 mbobra 1.6 
274            	if (nx1 > floor((nx-1)/scale + 1) )
275            		DIE("X-dimension of output array in fsample() is too large.");
276            	if (ny1 > floor((ny-1)/scale + 1) )
277            		DIE("Y-dimension of output array in fsample() is too large.");
278             
279                    // malloc some arrays for the R parameter calculation
280 mbobra 1.7         float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
281                    float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
282                    float *p1n0    = (float *)malloc(nx1*ny1*sizeof(float));
283                    float *p1p     = (float *)malloc(nx1*ny1*sizeof(float));
284                    float *p1n     = (float *)malloc(nx1*ny1*sizeof(float));
285                    float *p1      = (float *)malloc(nx1*ny1*sizeof(float));
286                    float *pmap    = (float *)malloc(nxp*nyp*sizeof(float));
287                    float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
288                    float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
289                
290 mbobra 1.1 	for (irec=0;irec<nrecs;irec++)
291            	{
292               	   // Get emphemeris
293 mbobra 1.7   	   sharpinrec    = sharpinrecset->records[irec];
294 mbobra 1.6   	   sharpceainrec = sharpceainrecset->records[irec];
295 mbobra 1.7 	   cdelt1_orig   = drms_getkey_float(sharpceainrec, "CDELT1",   &status);
296            	   dsun_obs      = drms_getkey_float(sharpinrec, "DSUN_OBS",   &status);
297            	   rsun_ref      = drms_getkey_double(sharpinrec, "RSUN_REF", &status);
298            	   rsun_obs      = drms_getkey_double(sharpinrec, "RSUN_OBS", &status);
299            	   imcrpix1      = drms_getkey_float(sharpinrec, "IMCRPIX1", &status);
300            	   imcrpix2      = drms_getkey_float(sharpinrec, "IMCRPIX2", &status);
301            	   crpix1        = drms_getkey_float(sharpinrec, "CRPIX1", &status);
302            	   crpix2        = drms_getkey_float(sharpinrec, "CRPIX2", &status);
303 mbobra 1.6 
304 mbobra 1.1   	   sharpoutrec = sharpoutrecset->records[irec];
305              	   sharpceaoutrec = sharpceaoutrecset->records[irec];
306 mbobra 1.6            if (!sameflag)
307                       {
308                         drms_copykeys(sharpoutrec, sharpinrec, 0, kDRMS_KeyClass_Explicit);
309                         drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit);
310                       }
311            
312                       TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
313                       tnow = (double)time(NULL);
314                       tnow += UNIX_epoch;
315 mbobra 1.1 
316 mbobra 1.6 	   // set CODEVER7 and HISTORY and DATE
317 mbobra 1.1 	   drms_setkey_string(sharpoutrec, "CODEVER7", cvsinfoall);
318                       drms_setkey_string(sharpoutrec,"HISTORY",historyofthemodule);           
319            	   drms_setkey_string(sharpceaoutrec, "CODEVER7", cvsinfoall);
320                       drms_setkey_string(sharpceaoutrec,"HISTORY",historyofthemodule);
321 mbobra 1.6            drms_setkey_time(sharpoutrec,"DATE",tnow);  
322                       drms_setkey_time(sharpceaoutrec,"DATE",tnow);
323 mbobra 1.1 
324            	   // Get bx, by, bz, mask	
325            	   // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
326 mbobra 1.6 	   DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");
327 mbobra 1.1 	   DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
328            	   int *bitmask = (int *) bitmaskArray->data;	 // get the previously made mask array	
329            
330            	   //Use conf_disambig map as a threshold on spaceweather quantities 
331 mbobra 1.6 	   DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");       
332 mbobra 1.1 	   DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);	
333                       if (status != DRMS_SUCCESS)
334                          DIE("No CONF_DISAMBIG segment.");
335            	   int *mask = (int *) maskArray->data;		// get the previously made mask array	
336            
337 mbobra 1.6  	   //Use magnetogram map to compute R
338            	   DRMS_Segment_t *losSeg = drms_segment_lookup(sharpceainrec, "magnetogram");
339            	   DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
340            	   float *los = (float *) losArray->data;		// los
341            
342            	   DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceainrec, "Bp");
343 mbobra 1.1 	   DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
344            	   float *bx = (float *) bxArray->data;		// bx
345            	
346 mbobra 1.6 	   DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceainrec, "Bt");
347 mbobra 1.1 	   DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);
348            	   float *by = (float *) byArray->data;		// by
349            	   for (int i = 0; i < nxny; i++) by[i] *= -1;
350            	
351 mbobra 1.6 	   DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceainrec, "Br");
352 mbobra 1.1 	   DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
353            	   float *bz = (float *) bzArray->data;		// bz
354            
355 mbobra 1.6    	   DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceainrec, "Br_err");
356 mbobra 1.1 	   DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
357            	   float *bz_err = (float *) bz_errArray->data;		// bz_err
358            
359 mbobra 1.6 	   DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceainrec, "Bt_err");
360 mbobra 1.1 	   DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
361            	   float *by_err = (float *) by_errArray->data;		// by_err
362            
363 mbobra 1.6 	   DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceainrec, "Bp_err");
364 mbobra 1.1 	   DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
365            	   float *bx_err = (float *) bx_errArray->data;		// bx_err
366            
367 mbobra 1.6            /***** USFLUX, Example Function 1 *************************************/ 
368 mbobra 1.1 	   if (meanvfflag)
369            	   {
370                          // Compute unsigned flux 
371                          if (computeAbsFlux(bz_err, bz , dims, &absFlux, &mean_vf,  &mean_vf_err, 
372                                             &count_mask, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
373                          {
374            	         mean_vf           = DRMS_MISSING_FLOAT;
375                             mean_vf_err       = DRMS_MISSING_FLOAT;
376                             count_mask        = DRMS_MISSING_INT;
377            	      }
378            
379                       drms_setkey_float(sharpoutrec, "USFLUX",  mean_vf);
380                       drms_setkey_float(sharpoutrec, "CMASK",   count_mask); 	
381                       drms_setkey_float(sharpoutrec, "ERRVF",   mean_vf_err);
382            
383                       drms_setkey_float(sharpceaoutrec, "USFLUX",  mean_vf);
384                       drms_setkey_float(sharpceaoutrec, "CMASK",   count_mask); 	
385                       drms_setkey_float(sharpceaoutrec, "ERRVF",   mean_vf_err);
386                       }
387 mbobra 1.6 
388                        /***** RPARAM, Example Function 14 *************************************/ 
389            	   if (r_valueflag)
390            	   {
391                          if (computeR(bz_err, los , dims, &Rparam, cdelt1, rim, p1p0, p1n0,
392 mbobra 1.7                            p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
393            
394 mbobra 1.6               {
395            		 Rparam            = DRMS_MISSING_FLOAT;
396            	      }
397            
398                       drms_setkey_float(sharpoutrec, "R_VALUE",  Rparam);
399                       drms_setkey_float(sharpceaoutrec, "R_VALUE", Rparam);
400            
401                       }
402              
403 mbobra 1.1            /***** MEANPOT and TOTPOT, Example Function 13  (Requires Keiji's code) **/
404                   
405               	   if (totpotflag || meanpotflag)
406            	   {
407                          // First compute potential field	
408                          for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
409            	      greenpot(bpx, bpy, bpz, nx, ny);
410             
411                          // Then compute energy   
412            	      if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims, 
413                                                &meanpot, &meanpot_err, &totpot, &totpot_err, 
414            			            mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 
415                          {
416            	         meanpot           = DRMS_MISSING_FLOAT; 
417            	         totpot            = DRMS_MISSING_FLOAT;
418                             meanpot_err       = DRMS_MISSING_FLOAT;
419                             totpot_err        = DRMS_MISSING_FLOAT;
420            	      }
421            	   
422                          // Set sharp keys          
423                          drms_setkey_float(sharpoutrec, "MEANPOT",  meanpot); 
424 mbobra 1.1               drms_setkey_float(sharpoutrec, "TOTPOT",   totpot); 	
425                          drms_setkey_float(sharpoutrec, "ERRMPOT",  meanpot_err); 
426                          drms_setkey_float(sharpoutrec, "ERRTPOT",  totpot_err);  
427            
428                          // Set sharp cea keys
429                          drms_setkey_float(sharpceaoutrec, "MEANPOT",  meanpot); 
430                          drms_setkey_float(sharpceaoutrec, "TOTPOT",   totpot); 	
431                          drms_setkey_float(sharpceaoutrec, "ERRMPOT",  meanpot_err); 
432                          drms_setkey_float(sharpceaoutrec, "ERRTPOT",  totpot_err);
433                       }
434                   
435                       /***** MEANGAM, Example Function 3 ************************************/
436               
437                       if (meangamflag)
438                       {
439                          // First compute horizontal field
440                          computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);
441            
442                          // Then compute inclination angle, gamma
443                          if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &mean_gamma, &mean_gamma_err,mask, bitmask))
444            	      {	
445 mbobra 1.1                  mean_gamma                 =  DRMS_MISSING_FLOAT;
446                             mean_gamma_err             =  DRMS_MISSING_FLOAT;
447                          }
448            
449                          // Set sharp keys    
450                          drms_setkey_float(sharpoutrec, "MEANGAM",  mean_gamma); 
451                          drms_setkey_float(sharpoutrec, "ERRGAM",   mean_gamma_err); 	
452            
453                          // Set sharp cea keys
454                          drms_setkey_float(sharpceaoutrec, "MEANGAM", mean_gamma); 
455                          drms_setkey_float(sharpceaoutrec, "ERRGAM",  mean_gamma_err); 	
456                       }
457            
458                       /***** MEANGBT, Example Function 5 (Requires Function 4) *************/
459            
460                       if (meangbtflag)
461                       {
462                          // First, compute Bt 
463                          computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);
464            
465                          // Then take the derivative of Bt
466 mbobra 1.1 	      if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask, bitmask, derx_bt, dery_bt, bt_err, 
467                                                      &mean_derivative_btotal_err))
468                          {
469            	         mean_derivative_btotal     = DRMS_MISSING_FLOAT;
470            	         mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
471                          }
472            
473                          // Set sharp keys    
474                          drms_setkey_float(sharpoutrec, "MEANGBT",  mean_derivative_btotal); 
475                          drms_setkey_float(sharpoutrec, "ERRBT",  mean_derivative_btotal_err); 
476            	
477                          // Set sharp cea keys  
478                          drms_setkey_float(sharpceaoutrec, "MEANGBT", mean_derivative_btotal); 
479                          drms_setkey_float(sharpceaoutrec, "ERRBT", mean_derivative_btotal_err); 	
480                       }
481               
482                       /***** MEANGBH, Example Function 6 (Requires Function 2) *************/
483            
484                       if (meangbhflag)
485                       {
486                          // First, compute Bh
487 mbobra 1.1    	      computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);
488            
489                          // Then take the derivative of Bh
490            	      if (computeBhderivative(bh, bh_err, dims, &mean_derivative_bh, &mean_derivative_bh_err, mask, bitmask, derx_bh, dery_bh))
491                          {
492            	         mean_derivative_bh       = DRMS_MISSING_FLOAT;
493                             mean_derivative_bh_err   = DRMS_MISSING_FLOAT;
494            	      }
495            
496                          // Set sharp keys   
497                          drms_setkey_float(sharpoutrec, "MEANGBH", mean_derivative_bh); 
498                          drms_setkey_float(sharpoutrec, "ERRBH", mean_derivative_bh_err); 
499            	
500                          // Set sharp cea keys  
501                          drms_setkey_float(sharpceaoutrec, "MEANGBH", mean_derivative_bh); 
502                          drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err); 	
503                       }
504            
505                       /***** MEANGBZ, Example Function 7 ************************************/
506            
507                       if (meangbzflag)
508 mbobra 1.1            {
509                          // Compute Bz derivative
510                          if (computeBzderivative(bz, bz_err, dims, &mean_derivative_bz, &mean_derivative_bz_err, mask, bitmask, derx_bz, dery_bz))
511                          {
512            	         mean_derivative_bz     = DRMS_MISSING_FLOAT; 
513                             mean_derivative_bz_err = DRMS_MISSING_FLOAT; 
514                          }
515            	
516                          // Set sharp keys   
517                          drms_setkey_float(sharpoutrec, "MEANGBZ",  mean_derivative_bz); 
518                          drms_setkey_float(sharpoutrec, "ERRBZ",  mean_derivative_bz_err); 
519            	
520                          // Set sharp cea keys   
521                          drms_setkey_float(sharpceaoutrec, "MEANGBZ", mean_derivative_bz); 
522                          drms_setkey_float(sharpceaoutrec, "ERRBZ", mean_derivative_bz_err); 	
523                       }
524            
525                       /***** MEANJZD and TOTUSJZ, Example Function 9 (Requires Function 8) ***/
526            
527                       if (meanjzdflag || totusjzflag)
528                       {
529 mbobra 1.1               // First, compute Jz on all pixels
530            	      computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
531                                 derx, dery);
532               
533                          // Then take the sums of the appropriate pixels
534                          if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &mean_jz,
535                                   &mean_jz_err, &us_i, &us_i_err, mask, bitmask, cdelt1,
536                                   rsun_ref, rsun_obs, derx, dery))
537            
538            
539                          {
540                             mean_jz                = DRMS_MISSING_FLOAT;
541            	         us_i                   = DRMS_MISSING_FLOAT;
542                             mean_jz_err            = DRMS_MISSING_FLOAT;
543                             us_i_err               = DRMS_MISSING_FLOAT;
544            	      }
545            
546                          // Set sharp keys   
547                          drms_setkey_float(sharpoutrec, "MEANJZD",  mean_jz); 
548                          drms_setkey_float(sharpoutrec, "TOTUSJZ",  us_i); 	
549                          drms_setkey_float(sharpoutrec, "ERRJZ",  mean_jz_err); 
550 mbobra 1.1               drms_setkey_float(sharpoutrec, "ERRUSI",  us_i_err); 	
551            
552                          // Set sharp cea keys   
553                          drms_setkey_float(sharpceaoutrec, "MEANJZD", mean_jz); 
554                          drms_setkey_float(sharpceaoutrec, "TOTUSJZ", us_i); 	
555                          drms_setkey_float(sharpceaoutrec, "ERRJZ", mean_jz_err); 
556                          drms_setkey_float(sharpceaoutrec, "ERRUSI", us_i_err); 	
557                       }
558            
559                       /***** MEANALP, Example Function 10 (Requires Function 8)*********/
560            
561                       if (meanalpflag)
562                       {
563                          // First, compute Jz on all pixels
564               	      computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
565                                    derx, dery);
566            
567                          // Then compute alpha quantities on select pixels
568            	      if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &mean_alpha, &mean_alpha_err, 
569                                           mask, bitmask, cdelt1, rsun_ref, rsun_obs))
570                          {
571 mbobra 1.1 	         mean_alpha             = DRMS_MISSING_FLOAT;
572                             mean_alpha_err         = DRMS_MISSING_FLOAT;
573                          }
574            
575                          // Set sharp keys   	
576                          drms_setkey_float(sharpoutrec, "MEANALP", mean_alpha); 
577                          drms_setkey_float(sharpoutrec, "ERRALP", mean_alpha_err); 	
578            
579                          // Set sharp cea keys   
580                          drms_setkey_float(sharpceaoutrec, "MEANALP", mean_alpha); 
581                          drms_setkey_float(sharpceaoutrec, "ERRALP", mean_alpha_err); 	
582            
583                       }
584            
585                       /***** MEANJZH, TOTUSJH, ABSNJZH, Example Function 11 (Requires Function 8) ***/
586            
587                       if (meanjzhflag || totusjhflag || absnjzhflag)
588                       {
589                          // First, compute Jz on all pixels
590                	      computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
591                                    derx, dery);
592 mbobra 1.1 
593                          // Then compute helicity quantities on select pixels
594                       if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &mean_ih, &mean_ih_err, &total_us_ih, 
595                                           &total_abs_ih, &total_us_ih_err, &total_abs_ih_err, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
596                          {  
597            	         mean_ih                = DRMS_MISSING_FLOAT; 
598            	         total_us_ih            = DRMS_MISSING_FLOAT;
599              	         total_abs_ih           = DRMS_MISSING_FLOAT;
600                             mean_ih_err            = DRMS_MISSING_FLOAT;
601                             total_us_ih_err        = DRMS_MISSING_FLOAT;
602                             total_abs_ih_err       = DRMS_MISSING_FLOAT;
603            	      }
604            
605                          // Set sharp keys 
606                          drms_setkey_float(sharpoutrec, "MEANJZH",  mean_ih); 
607                          drms_setkey_float(sharpoutrec, "TOTUSJH",  total_us_ih); 	
608                          drms_setkey_float(sharpoutrec, "ABSNJZH",  total_abs_ih); 
609                          drms_setkey_float(sharpoutrec, "ERRMIH",  mean_ih_err); 	
610                          drms_setkey_float(sharpoutrec, "ERRTUI",  total_us_ih_err); 
611                          drms_setkey_float(sharpoutrec, "ERRTAI",  total_abs_ih_err); 	
612            
613 mbobra 1.1               // Set sharp cea keys 
614                          drms_setkey_float(sharpceaoutrec, "MEANJZH", mean_ih); 
615                          drms_setkey_float(sharpceaoutrec, "TOTUSJH", total_us_ih); 	
616                          drms_setkey_float(sharpceaoutrec, "ABSNJZH", total_abs_ih); 
617                          drms_setkey_float(sharpceaoutrec, "ERRMIH", mean_ih_err); 	
618                          drms_setkey_float(sharpceaoutrec, "ERRTUI", total_us_ih_err); 
619                          drms_setkey_float(sharpceaoutrec, "ERRTAI", total_abs_ih_err); 	
620                       }
621            
622                       /***** SAVNCPP, Example Function 12 (Requires Function 8) *******************/
623               
624                       if (savncppflag)
625                       {
626                          // First, compute Jz on all pixels
627              	      computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
628                                    derx, dery);
629            
630                          // Then compute sums of Jz on select pixels
631              	      if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &totaljz, &totaljz_err, 
632            				           mask, bitmask, cdelt1, rsun_ref, rsun_obs))
633                          {
634 mbobra 1.1 	         totaljz                = DRMS_MISSING_FLOAT;
635                             totaljz_err            = DRMS_MISSING_FLOAT;
636            	      }
637            
638                          // Set sharp keys 	
639                          drms_setkey_float(sharpoutrec, "SAVNCPP",  totaljz); 
640                          drms_setkey_float(sharpoutrec, "ERRJHT",  totaljz_err); 	
641            
642                          // Set sharp cea keys 
643                          drms_setkey_float(sharpceaoutrec, "SAVNCPP", totaljz); 
644                          drms_setkey_float(sharpceaoutrec, "ERRJHT", totaljz_err); 	
645                       }
646            
647                       /***** MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) **********/
648            
649                       if (meanshrflag || shrgt45flag)
650                       {
651                          // First compute potential field
652            	      for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
653            	      greenpot(bpx, bpy, bpz, nx, ny); 
654            
655 mbobra 1.1               // Then compute shear angles
656            	      if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims, 
657            		   		    &meanshear_angle, &meanshear_angle_err, &area_w_shear_gt_45,  
658            				    mask, bitmask)) 
659                          {
660            	          meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
661            		  area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
662                              meanshear_angle_err= DRMS_MISSING_FLOAT;
663            	      }
664            
665                          // Set sharp keys 	
666                          drms_setkey_float(sharpoutrec, "MEANSHR",  meanshear_angle); 
667                          drms_setkey_float(sharpoutrec, "SHRGT45",  area_w_shear_gt_45); 	
668                          drms_setkey_float(sharpoutrec, "ERRMSHA",  meanshear_angle_err); 	
669            
670                          // Set sharp cea keys 
671                          drms_setkey_float(sharpceaoutrec, "MEANSHR", meanshear_angle); 
672                          drms_setkey_float(sharpceaoutrec, "SHRGT45", area_w_shear_gt_45); 	
673                          drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err); 	  
674                       }
675            
676 mbobra 1.1            /******************************* END FLAGS **********************************/
677            
678            	   drms_free_array(bitmaskArray);		
679            	   drms_free_array(maskArray);
680            	   drms_free_array(bxArray);           
681            	   drms_free_array(byArray);
682            	   drms_free_array(bzArray);
683            	   drms_free_array(bx_errArray);           
684            	   drms_free_array(by_errArray);
685            	   drms_free_array(bz_errArray);
686 mbobra 1.6 	   drms_free_array(losArray);
687 mbobra 1.1 	} //endfor
688            
689 mbobra 1.6  	drms_close_records(sharpinrecset, DRMS_FREE_RECORD);
690            	drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);
691            
692 mbobra 1.1 	drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);
693            	drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);
694               
695                    //used for select calculations
696            	free(bh); free(bt); free(jz); 
697            	free(bpx); free(bpy); free(bpz); 
698            	free(derx); free(dery);	 
699            	free(derx_bt); free(dery_bt); 	
700            	free(derx_bz); free(dery_bz);	 
701            	free(derx_bh); free(dery_bh); 
702            	free(bt_err); free(bh_err);  free(jz_err); 
703                    free(jz_err_squared); free(jz_rms_err);
704            	free(cvsinfoall);
705 mbobra 1.7 
706 mbobra 1.6         free(rim);
707                    free(p1p0);
708                    free(p1n0);
709                    free(p1p);
710                    free(p1n);
711                    free(p1);
712                    free(pmap);
713 mbobra 1.7         free(p1pad);
714                    free(pmapn);
715 mbobra 1.6 
716 arta   1.3 return 0;    
717 mbobra 1.1 }	// DoIt

Karen Tian
Powered by
ViewCVS 0.9.4