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

Karen Tian
Powered by
ViewCVS 0.9.4