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

File: [Development] / JSOC / proj / sharp / apps / update_sharp_keys.c (download)
Revision: 1.19, Wed May 26 18:25:17 2021 UTC (16 months, 1 week ago) by mbobra
Branch: MAIN
CVS Tags: Ver_LATEST, Ver_9-5, HEAD
Changes since 1.18: +88 -96 lines
Figured out the memory allocation problem.
It turns out that the arrays are not constant per T_REC for any given HARPNUM, e.g.:

> show_info -d key=T_REC ds=hmi.sharp_cea_720s"[3311][]" seg={Bt}
T_REC	Bt	Bt_info
2013.10.24_00:48:00_TAI	Bt.fits	"Gauss, fits, 547x411"
2013.10.24_01:00:00_TAI	Bt.fits	"Gauss, fits, 547x411"
2013.10.24_01:12:00_TAI	Bt.fits	"Gauss, fits, 546x411"

So I moved all of the code that created one fixed array size per HARPNUM into the loop.
Now it creates a new array per T_REC.

/*                                                                      
 *      MODULE NAME: update_sharp_keys.c
 *
 *	    DESCRIPTION: This module recalculates SHARP keywords. 
 *      This is accomplished by cloning a record, recalculating keywords of choice 
 *      (i.e., user input), and pointing to the same segments as the old record. 
 *      Associated error keys are computed by default. 
 * 
 *      This module accounts for versioning by appending to the keyword CODEVER7 and HISTORY:
 *      CODEVER7 will contain multiple lines: the production build of sharp.c, the 
 *      production build of include file sw_functions.c, and the production build
 *      of update_sharp_keys.c. HISTORY will include a human-readable sentence about which
 *      keywords were updated. 
 * 
 *      N.B.: This module calculates the keyword specified in keylist and then does a 
 *      drms_copykeys() to copy over all the remaining keywords. However, if [1] a keyword
 *      does not exist in the .jsd files of the input series, and [2] said keyword X is not
 *      specified in keylist, then [3] drms_copykeys() will copy DRMS_MISSING_VALUE into 
 *      keyword X for all the output series.
 *
 *      This module does not produce segments.
 *
 *      INPUTS     : -- DRMS SHARP series
 *                   -- DRMS SHARP CEA series
 *                   -- HARPNUM
 *                   -- comma separated list of keywords to recalculate
 *                   -- DEBUG flag
 *
 *	AUTHOR     : Monica Bobra
 *
 *	EXAMPLE    :
 *      update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //
 *      HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT debug=debug
 *
 */

/* Include files */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>
#include <string.h>
#include "jsoc_main.h"
#include "astro.h"
#include "fstats.h"
#include "cartography.c"
#include "fresize.h"
#include "finterpolate.h"
#include "img2helioVector.c"
#include "copy_me_keys.c"
#include "errorprop.c"
#include "sw_functions.c"
#include <mkl_blas.h>
#include <mkl_service.h>
#include <mkl_lapack.h>
#include <mkl_vml_functions.h>
#include <omp.h>
 
/* Define values */                     
#define PI              (M_PI)
#define RADSINDEG       (PI/180.)
#define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
#define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
#define SHOW(msg) {printf("%s", msg); fflush(stdout);}
#define DEBUG(true)

char *module_name = "update_sharp_keys";  /* Module name */
char *version_id  = "2021 May 25";        /* Version number */

ModuleArgs_t module_args[] =
{
	{ARG_STRING, "sharpseriesin",      NULL, "Input Sharp dataseries"},
	{ARG_STRING, "sharpceaseriesin",   NULL, "Input Sharp CEA dataseries"},
	{ARG_STRING, "sharpseriesout",     NULL, "Output Sharp dataseries"},
	{ARG_STRING, "sharpceaseriesout",  NULL, "Output Sharp CEA dataseries"},
	{ARG_INT,    "HARPNUM",            NULL, "HARP number"},
	{ARG_STRING, "keylist",            NULL, "comma separated list of keywords to update (for LOS keywords enter MEANGBZLOS or USFLUXLOS)"},
	{ARG_STRING, "debug",              NULL, "debug mode functionality. use like this: debug=debug"},
	{ARG_END}
};


int DoIt(void)                    
{
    int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
    int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);

	int status = DRMS_SUCCESS;
	int nrecs, irec;

    /* Keywords */
	float mean_vf;
    float mean_vf_los; 
	float absFlux;
    float absFlux_los;
    float count_mask;
    float count_mask_los;
	float mean_hf;
	float mean_gamma;
	float mean_derivative_btotal;
	float mean_derivative_bh; 
	float mean_derivative_bz;
	float mean_derivative_los;
	float mean_jz;
	float us_i;
	float mean_alpha;
	float mean_ih;
	float total_us_ih;
	float total_abs_ih;
	float totaljz;
	float totpot;
	float meanpot;
	float area_w_shear_gt_45;
	float meanshear_angle;
	float area_w_shear_gt_45h;
	float meanshear_angleh; 
    float mean_derivative_btotal_err;
    float mean_vf_err;
    float mean_gamma_err;
    float mean_derivative_bh_err;
    float mean_derivative_bz_err;
    float mean_jz_err;
    float us_i_err; 
    float mean_alpha_err;
    float mean_ih_err;
    float total_us_ih_err;
    float total_abs_ih_err;
    float totaljz_err;
    float meanpot_err;
    float totpot_err;
    float Rparam;
    float meanshear_angle_err;
    float totfx;
    float totfy;
    float totfz;
    float totbsq;
    float epsx;
    float epsy;
    float epsz;

	char *sharpseriesin = (char *) params_get_str(&cmdparams, "sharpseriesin");
	char *sharpceaseriesin = (char *) params_get_str(&cmdparams, "sharpceaseriesin");
	char *sharpseriesout = (char *) params_get_str(&cmdparams, "sharpseriesout");
	char *sharpceaseriesout = (char *) params_get_str(&cmdparams, "sharpceaseriesout");

    int sameflag = strcmp(sharpseriesin, sharpseriesout) == 0;
    int testflag = strcmp(sharpceaseriesin, sharpceaseriesout) == 0;
    if (testflag != sameflag)
    DIE("Either both outputs must be the same as their inputs, or both be different.");

	int harpnum = params_get_int(&cmdparams, "HARPNUM");

	char sharpquery[256];
	char sharpceaquery[256];
	sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
	printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);

	sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
	sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);

	DRMS_RecordSet_t *sharpinrecset= drms_open_records(drms_env, sharpquery, &status);
	if (status != DRMS_SUCCESS)
		DIE("Problem opening sharp recordset.");
	nrecs=sharpinrecset->n;
	if (nrecs == 0)
		DIE("Sharp recordset contains no records.");
        printf("nrecs=%d\n",nrecs);

	DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status);
	if (status != DRMS_SUCCESS)
		DIE("Problem opening sharp cea recordset.");
	if (sharpceainrecset->n != nrecs)
		DIE("Sharp cea recordset differs from sharp recordset.");

        DRMS_RecordSet_t *sharpoutrecset, *sharpceaoutrecset;
        if (sameflag)
        {
	sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
	if (status != DRMS_SUCCESS)
		DIE("Problem cloning sharp records.");
	sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
	if (status != DRMS_SUCCESS)
		DIE("Problem cloning sharp cea records.");
        }
        else
        {
	sharpoutrecset = drms_create_records(drms_env, nrecs, sharpseriesout, DRMS_PERMANENT, &status);
	if (status != DRMS_SUCCESS)
		DIE("Problem creating sharp records.");
	sharpceaoutrecset = drms_create_records(drms_env, nrecs, sharpceaseriesout, DRMS_PERMANENT, &status);
	if (status != DRMS_SUCCESS)
		DIE("Problem creating sharp cea records.");
        }

	char *keylist = (char *) params_get_str(&cmdparams, "keylist");
	char *debug   = (char *) params_get_str(&cmdparams, "debug");

    // Flags to indicate which keyword will be recalculated
    int meanvflosflag = (strstr(keylist,"USFLUXL")  != NULL);
    int meanglosflag  = (strstr(keylist,"MEANGBL") != NULL);
    int meanvfflag    = (strstr(keylist,"USFLUX")  != NULL); 
    int totpotflag    = (strstr(keylist,"TOTPOT")  != NULL);  
    int meangamflag   = (strstr(keylist,"MEANGAM") != NULL);
    int meangbtflag   = (strstr(keylist,"MEANGBT") != NULL);
    int meangbzflag   = (strstr(keylist,"MEANGBZ") != NULL);
    int meangbhflag   = (strstr(keylist,"MEANGBH") != NULL);
    int meanjzdflag   = (strstr(keylist,"MEANJZD") != NULL);
    int totusjzflag   = (strstr(keylist,"TOTUSJZ") != NULL);
    int meanalpflag   = (strstr(keylist,"MEANALP") != NULL);
    int meanjzhflag   = (strstr(keylist,"MEANJZH") != NULL);
    int totusjhflag   = (strstr(keylist,"TOTUSJH") != NULL);
    int absnjzhflag   = (strstr(keylist,"ABSNJZH") != NULL);
    int savncppflag   = (strstr(keylist,"SAVNCPP") != NULL);
    int meanpotflag   = (strstr(keylist,"MEANPOT") != NULL);
    int meanshrflag   = (strstr(keylist,"MEANSHR") != NULL);
    int shrgt45flag   = (strstr(keylist,"SHRGT45") != NULL);
    int r_valueflag   = (strstr(keylist,"R_VALUE") != NULL);
    int totbsqflag    = (strstr(keylist,"TOTBSQ")  != NULL);
    int totfxflag     = (strstr(keylist,"TOTFX")   != NULL);
    int totfyflag     = (strstr(keylist,"TOTFY")   != NULL);
    int totfzflag     = (strstr(keylist,"TOTFZ")   != NULL);
    int epsxflag      = (strstr(keylist,"EPSX")    != NULL);
    int epsyflag      = (strstr(keylist,"EPSY")    != NULL);
    int epszflag      = (strstr(keylist,"EPSZ")    != NULL);
    int debugflag     = (strstr(debug,"debug")     != NULL);
    
	for (irec=0;irec<nrecs;irec++)
	{
        DRMS_Record_t *sharpinrec = sharpinrecset->records[irec];
	    DRMS_Record_t *sharpceainrec = sharpceainrecset->records[irec];
	    DRMS_Segment_t *inseg = drms_segment_lookup(sharpceainrec, "Br");
	    int nx = inseg->axis[0];
	    int ny = inseg->axis[1];
	    int nxny = nx * ny;
	    int dims[2] = {nx, ny};

	    // Temp arrays 	
	    float *bh      = (float *) (malloc(nxny * sizeof(float)));
	    float *bt      = (float *) (malloc(nxny * sizeof(float)));
	    float *jz      = (float *) (malloc(nxny * sizeof(float)));
	    float *bpx     = (float *) (malloc(nxny * sizeof(float)));
	    float *bpy     = (float *) (malloc(nxny * sizeof(float)));
	    float *bpz     = (float *) (malloc(nxny * sizeof(float)));
	    float *derx    = (float *) (malloc(nxny * sizeof(float)));
	    float *dery    = (float *) (malloc(nxny * sizeof(float)));
	    float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
	    float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
	    float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
	    float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
	    float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
	    float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
	    float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
	    float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
        float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
        float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
        float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
        float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
	    float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
	    float *err_term1 = (float *) (calloc(nxny, sizeof(float)));
	    float *err_term2 = (float *) (calloc(nxny, sizeof(float)));
        float *fx        = (float *) (malloc(nxny * sizeof(float)));
        float *fy        = (float *) (malloc(nxny * sizeof(float)));
        float *fz        = (float *) (malloc(nxny * sizeof(float)));
        float *derx_los  = (float *) (malloc(nxny * sizeof(float)));
	    float *dery_los  = (float *) (malloc(nxny * sizeof(float)));
 
        // ephemeris variables
	    float  cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2;
	    double rsun_ref, rsun_obs;
         
        // outrecords 
	    DRMS_Record_t *sharpoutrec, *sharpceaoutrec, *testrec;

        // prepare to set CODEVER7 (CVS Version of the SHARP module)
	    char *cvsinfo0;
	    char *history0;
	    char *cvsinfo1 = strdup("$Id: update_sharp_keys.c,v 1.19 2021/05/26 19:25:17 mbobra Exp $");
	    char *cvsinfo2 = sw_functions_version();
	    char *cvsinfoall = (char *)malloc(2048);
            char historyofthemodule[2048];
            cvsinfo0    = drms_getkey_string(sharpinrec, "CODEVER7", &status);
            strcpy(cvsinfoall,cvsinfo0);
            strcat(cvsinfoall,"\n");
            strcat(cvsinfoall,cvsinfo2);

        // prepare to set HISTORY (History of the data)
	    char timebuf[1024];
	    float UNIX_epoch = -220924792.000; 	
	    sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
            sprintf(historyofthemodule,"The following keywords were re-computed on %s: %s.",timebuf,keylist); 
	    //printf("historyofthemodule=%s\n",historyofthemodule);
	    history0    = drms_getkey_string(sharpinrec, "HISTORY", &status);
	    strcat(historyofthemodule,"\n");
            strcat(historyofthemodule,history0);
 
        // grab one record such that we can malloc the R-parameter calculation arrays outside of the loop
        // assume scale = round(2./cdelt1) does not change for any given HARP
        testrec  = sharpceainrecset->records[nrecs/2];
        dsun_obs = drms_getkey_float(testrec, "DSUN_OBS",   &status); 
        rsun_ref = drms_getkey_double(testrec, "RSUN_REF", &status);
        cdelt1_orig = drms_getkey_float(testrec, "CDELT1",   &status);
        int scale = round(2.0/cdelt1);
        int nx1 = nx/scale;
        int ny1 = ny/scale;
        int nxp = nx1+40;
        int nyp = ny1+40;

        // check to make sure that the dimensions passed into fsample() are adequate
	    if (nx1 > floor((nx-1)/scale + 1) )
		DIE("X-dimension of output array in fsample() is too large.");
	    if (ny1 > floor((ny-1)/scale + 1) )
        	DIE("Y-dimension of output array in fsample() is too large.");
        // malloc some arrays for the R parameter calculation
        float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
        float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
        float *p1n0    = (float *)malloc(nx1*ny1*sizeof(float));
        float *p1p     = (float *)malloc(nx1*ny1*sizeof(float));
        float *p1n     = (float *)malloc(nx1*ny1*sizeof(float));
        float *p1      = (float *)malloc(nx1*ny1*sizeof(float));
        float *pmap    = (float *)malloc(nxp*nyp*sizeof(float));
        float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
        float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
    
        // Get emphemeris
        sharpinrec    = sharpinrecset->records[irec];
        sharpceainrec = sharpceainrecset->records[irec];
        cdelt1_orig   = drms_getkey_float(sharpceainrec, "CDELT1",   &status);
        dsun_obs      = drms_getkey_float(sharpinrec, "DSUN_OBS",   &status);
        rsun_ref      = drms_getkey_double(sharpinrec, "RSUN_REF", &status);
        rsun_obs      = drms_getkey_double(sharpinrec, "RSUN_OBS", &status);
        imcrpix1      = drms_getkey_float(sharpinrec, "IMCRPIX1", &status);
        imcrpix2      = drms_getkey_float(sharpinrec, "IMCRPIX2", &status);
        crpix1        = drms_getkey_float(sharpinrec, "CRPIX1", &status);
        crpix2        = drms_getkey_float(sharpinrec, "CRPIX2", &status);
        cdelt1      = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);

        sharpoutrec = sharpoutrecset->records[irec];
        sharpceaoutrec = sharpceaoutrecset->records[irec];
           if (!sameflag)
           {
             drms_copykeys(sharpoutrec, sharpinrec, 0, kDRMS_KeyClass_Explicit);
             drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit);
           }

           TIME trec, tnow; /* 1970.01.01_00:00:00_UTC */
           tnow = (double)time(NULL);
           tnow += UNIX_epoch;

           if (debugflag)
           {
              printf("i=%d\n",irec);
           }

            // set CODEVER7 and HISTORY and DATE
	    drms_setkey_string(sharpoutrec, "CODEVER7", cvsinfoall);
            drms_setkey_string(sharpoutrec,"HISTORY",historyofthemodule);           
	    drms_setkey_string(sharpceaoutrec, "CODEVER7", cvsinfoall);
            drms_setkey_string(sharpceaoutrec,"HISTORY",historyofthemodule);
            drms_setkey_time(sharpoutrec,"DATE",tnow);  
            drms_setkey_time(sharpceaoutrec,"DATE",tnow);

	    // Get all the segments if computing vector magnetic field keywords

	    DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");
	    DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
	    int *bitmask = (int *) bitmaskArray->data;	 // get the previously made mask array	
 
	    DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");       
	    DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);	
            if (status != DRMS_SUCCESS)
            DIE("No CONF_DISAMBIG segment.");
	    int *mask = (int *) maskArray->data;		// get the previously made mask array	

	    DRMS_Segment_t *losSeg = drms_segment_lookup(sharpceainrec, "magnetogram");
	    DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
            if (status != DRMS_SUCCESS)
            DIE("No magnetogram segment.");
	    float *los = (float *) losArray->data;		// los
            //printf("Line 387");

	    DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceainrec, "Bp");
	    DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
            if (status != DRMS_SUCCESS)
            DIE("No Bp segment.");
	    float *bx = (float *) bxArray->data;		// bx
            //printf("Line 393");	

	    DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceainrec, "Bt");
	    DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);
	    if (status != DRMS_SUCCESS)
            DIE("No Bt segment.");
            float *by = (float *) byArray->data;		// by
	    for (int i = 0; i < nxny; i++) by[i] *= -1;
            //printf("Line 399");
	
	    DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceainrec, "Br");
	    //printf("Line 408\n");
            if (bzSeg == NULL)
            DIE("somekind of pointer problem.");            

            DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
	    //printf("status=%d\n",status);
            if (status != DRMS_SUCCESS)
            DIE("No Br segment.");
            float *bz = (float *) bzArray->data;		// bz
            //printf("Line 414");
   	    
            DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceainrec, "Br_err");
	    DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
	    if (status != DRMS_SUCCESS)
            DIE("No Br_err segment.");
            float *bz_err = (float *) bz_errArray->data;		// bz_err
            //printf("Line 408");
	    
            DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceainrec, "Bt_err");
	    DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
	    if (status != DRMS_SUCCESS)
            DIE("No Bt_err segment.");
            float *by_err = (float *) by_errArray->data;		// by_err
            //printf("Line 412");
	    
            DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceainrec, "Bp_err");
	    DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
	    if (status != DRMS_SUCCESS)
            DIE("No Bp_err segment.");
            float *bx_err = (float *) bx_errArray->data;		// bx_err

        /***** USFLUX, Example Function 1 *************************************/ 
	    if (meanvfflag)
	    {
            // Compute unsigned flux 
            if (computeAbsFlux(bz_err, bz , dims, &absFlux, &mean_vf,  &mean_vf_err, 
                               &count_mask, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
            {
                mean_vf = DRMS_MISSING_FLOAT;
                mean_vf_err = DRMS_MISSING_FLOAT;
                count_mask = DRMS_MISSING_INT;
            }

            drms_setkey_float(sharpoutrec, "USFLUX",  mean_vf);
            drms_setkey_float(sharpoutrec, "CMASK",   count_mask); 	
            drms_setkey_float(sharpoutrec, "ERRVF",   mean_vf_err);

            drms_setkey_float(sharpceaoutrec, "USFLUX",  mean_vf);
            drms_setkey_float(sharpceaoutrec, "CMASK",   count_mask); 	
            drms_setkey_float(sharpceaoutrec, "ERRVF",   mean_vf_err);
        }

        /***** RPARAM, Example Function 14 *************************************/ 
	    if (r_valueflag)
	    {
            if (computeR(bz_err, los , dims, &Rparam, cdelt1, rim, p1p0, p1n0,
                        p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))

            {
		        Rparam = DRMS_MISSING_FLOAT;
	        }

            drms_setkey_float(sharpoutrec, "R_VALUE",  Rparam);
            drms_setkey_float(sharpceaoutrec, "R_VALUE", Rparam);
        }

        /***** MEANGBZ, Example Function 7 ************************************/
        if (meangbzflag)
        {
            // Compute Bz derivative
            if (computeBzderivative(bz, bz_err, dims, &mean_derivative_bz, &mean_derivative_bz_err, mask, bitmask, derx_bz, dery_bz, err_term1, err_term2))
            {
                mean_derivative_bz     = DRMS_MISSING_FLOAT; 
                mean_derivative_bz_err = DRMS_MISSING_FLOAT; 
            }
	
            // Set sharp keys   
            drms_setkey_float(sharpoutrec, "MEANGBZ",  mean_derivative_bz); 
            drms_setkey_float(sharpoutrec, "ERRBZ",  mean_derivative_bz_err); 
	
            // Set sharp cea keys   
            drms_setkey_float(sharpceaoutrec, "MEANGBZ", mean_derivative_bz); 
            drms_setkey_float(sharpceaoutrec, "ERRBZ", mean_derivative_bz_err); 	
        }
  
        /***** MEANPOT and TOTPOT, Example Function 13  (Requires Keiji's code) **/
   	    if (totpotflag || meanpotflag)
        {
            // First compute potential field	
            for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
                greenpot(bpx, bpy, bpz, nx, ny);
 
                // Then compute energy   
                if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims, 
                                      &meanpot, &meanpot_err, &totpot, &totpot_err, 
			                          mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 
                {
                    meanpot = DRMS_MISSING_FLOAT; 
                    totpot = DRMS_MISSING_FLOAT;
                    meanpot_err = DRMS_MISSING_FLOAT;
                    totpot_err = DRMS_MISSING_FLOAT;
	            }
	   
            // Set sharp keys          
            drms_setkey_float(sharpoutrec, "MEANPOT",  meanpot); 
            drms_setkey_float(sharpoutrec, "TOTPOT",   totpot); 	
            drms_setkey_float(sharpoutrec, "ERRMPOT",  meanpot_err); 
            drms_setkey_float(sharpoutrec, "ERRTPOT",  totpot_err);  

            // Set sharp cea keys
            drms_setkey_float(sharpceaoutrec, "MEANPOT",  meanpot); 
            drms_setkey_float(sharpceaoutrec, "TOTPOT",   totpot); 	
            drms_setkey_float(sharpceaoutrec, "ERRMPOT",  meanpot_err); 
            drms_setkey_float(sharpceaoutrec, "ERRTPOT",  totpot_err);
        }
       
        /***** MEANGAM, Example Function 3 ************************************/
        if (meangamflag)
        {
            // First compute horizontal field
            computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);

            // Then compute inclination angle, gamma
            if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &mean_gamma, &mean_gamma_err,mask, bitmask))
            {	
                mean_gamma = DRMS_MISSING_FLOAT;
                mean_gamma_err = DRMS_MISSING_FLOAT;
            }

            // Set sharp keys    
            drms_setkey_float(sharpoutrec, "MEANGAM",  mean_gamma); 
            drms_setkey_float(sharpoutrec, "ERRGAM",   mean_gamma_err); 	

            // Set sharp cea keys
            drms_setkey_float(sharpceaoutrec, "MEANGAM", mean_gamma); 
            drms_setkey_float(sharpceaoutrec, "ERRGAM",  mean_gamma_err); 	
        }

        /***** MEANGBT, Example Function 5 (Requires Function 4) *************/
        if (meangbtflag)
        {
            // First, compute Bt 
            computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);

            // Then take the derivative of Bt
	        if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask, bitmask, derx_bt, dery_bt, bt_err, 
                                        &mean_derivative_btotal_err, err_term1, err_term2))
            {
	            mean_derivative_btotal     = DRMS_MISSING_FLOAT;
	            mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
            }

            // Set sharp keys    
            drms_setkey_float(sharpoutrec, "MEANGBT",  mean_derivative_btotal); 
            drms_setkey_float(sharpoutrec, "ERRBT",  mean_derivative_btotal_err); 
	
            // Set sharp cea keys  
            drms_setkey_float(sharpceaoutrec, "MEANGBT", mean_derivative_btotal); 
            drms_setkey_float(sharpceaoutrec, "ERRBT", mean_derivative_btotal_err); 	
        }
           
        /***** MEANGBH, Example Function 6 (Requires Function 2) *************/
        if (meangbhflag)
        {
            // First, compute Bh
            computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask);

            // Then take the derivative of Bh
	        if (computeBhderivative(bh, bh_err, dims, &mean_derivative_bh, &mean_derivative_bh_err, mask, bitmask, derx_bh, dery_bh, err_term1, err_term2))
                {
	        mean_derivative_bh       = DRMS_MISSING_FLOAT;
                mean_derivative_bh_err   = DRMS_MISSING_FLOAT;
	        }

            // Set sharp keys   
            drms_setkey_float(sharpoutrec, "MEANGBH", mean_derivative_bh); 
            drms_setkey_float(sharpoutrec, "ERRBH", mean_derivative_bh_err); 
	
            // Set sharp cea keys  
            drms_setkey_float(sharpceaoutrec, "MEANGBH", mean_derivative_bh); 
            drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err); 	
        }

        /***** MEANJZD and TOTUSJZ, Example Function 9 (Requires Function 8) ***/
        if (meanjzdflag || totusjzflag)
        {
            // First, compute Jz on all pixels
	        computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
                      derx, dery, err_term1, err_term2);
   
            // Then take the sums of the appropriate pixels
            if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &mean_jz,
                               &mean_jz_err, &us_i, &us_i_err, mask, bitmask, cdelt1,
                               rsun_ref, rsun_obs, derx, dery))

                {
                mean_jz                = DRMS_MISSING_FLOAT;
	        us_i                   = DRMS_MISSING_FLOAT;
                mean_jz_err            = DRMS_MISSING_FLOAT;
                us_i_err               = DRMS_MISSING_FLOAT;
	        }

            // Set sharp keys   
            drms_setkey_float(sharpoutrec, "MEANJZD",  mean_jz); 
            drms_setkey_float(sharpoutrec, "TOTUSJZ",  us_i); 	
            drms_setkey_float(sharpoutrec, "ERRJZ",  mean_jz_err); 
            drms_setkey_float(sharpoutrec, "ERRUSI",  us_i_err); 	

            // Set sharp cea keys   
            drms_setkey_float(sharpceaoutrec, "MEANJZD", mean_jz); 
            drms_setkey_float(sharpceaoutrec, "TOTUSJZ", us_i); 	
            drms_setkey_float(sharpceaoutrec, "ERRJZ", mean_jz_err); 
            drms_setkey_float(sharpceaoutrec, "ERRUSI", us_i_err); 	
        }

        /***** MEANALP, Example Function 10 (Requires Function 8)*********/
        if (meanalpflag)
        {
            // First, compute Jz on all pixels
   	        computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
                      derx, dery, err_term1, err_term2);

            // Then compute alpha quantities on select pixels
	        if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &mean_alpha, &mean_alpha_err, 
                             mask, bitmask, cdelt1, rsun_ref, rsun_obs))
                {
	        mean_alpha             = DRMS_MISSING_FLOAT;
                mean_alpha_err         = DRMS_MISSING_FLOAT;
                }

            // Set sharp keys   	
            drms_setkey_float(sharpoutrec, "MEANALP", mean_alpha); 
            drms_setkey_float(sharpoutrec, "ERRALP", mean_alpha_err); 	

            // Set sharp cea keys   
            drms_setkey_float(sharpceaoutrec, "MEANALP", mean_alpha); 
            drms_setkey_float(sharpceaoutrec, "ERRALP", mean_alpha_err); 	

        }

        /***** MEANJZH, TOTUSJH, ABSNJZH, Example Function 11 (Requires Function 8) ***/
        if (meanjzhflag || totusjhflag || absnjzhflag)
        {
            // First, compute Jz on all pixels
    	    computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
                      derx, dery, err_term1, err_term2);

            // Then compute helicity quantities on select pixels
            if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &mean_ih, &mean_ih_err, &total_us_ih, 
                                &total_abs_ih, &total_us_ih_err, &total_abs_ih_err, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
                {  
	        mean_ih                = DRMS_MISSING_FLOAT; 
	        total_us_ih            = DRMS_MISSING_FLOAT;
  	        total_abs_ih           = DRMS_MISSING_FLOAT;
                mean_ih_err            = DRMS_MISSING_FLOAT;
                total_us_ih_err        = DRMS_MISSING_FLOAT;
                total_abs_ih_err       = DRMS_MISSING_FLOAT;
	        }

            // Set sharp keys 
            drms_setkey_float(sharpoutrec, "MEANJZH",  mean_ih); 
            drms_setkey_float(sharpoutrec, "TOTUSJH",  total_us_ih); 	
            drms_setkey_float(sharpoutrec, "ABSNJZH",  total_abs_ih); 
            drms_setkey_float(sharpoutrec, "ERRMIH",  mean_ih_err); 	
            drms_setkey_float(sharpoutrec, "ERRTUI",  total_us_ih_err); 
            drms_setkey_float(sharpoutrec, "ERRTAI",  total_abs_ih_err); 	

            // Set sharp cea keys 
            drms_setkey_float(sharpceaoutrec, "MEANJZH", mean_ih); 
            drms_setkey_float(sharpceaoutrec, "TOTUSJH", total_us_ih); 	
            drms_setkey_float(sharpceaoutrec, "ABSNJZH", total_abs_ih); 
            drms_setkey_float(sharpceaoutrec, "ERRMIH", mean_ih_err); 	
            drms_setkey_float(sharpceaoutrec, "ERRTUI", total_us_ih_err); 
            drms_setkey_float(sharpceaoutrec, "ERRTAI", total_abs_ih_err); 	
        }

        /***** SAVNCPP, Example Function 12 (Requires Function 8) *******************/
        if (savncppflag)
        {
            // First, compute Jz on all pixels
  	        computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 
                      derx, dery, err_term1, err_term2);

            // Then compute sums of Jz on select pixels
  	        if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &totaljz, &totaljz_err, 
				                         mask, bitmask, cdelt1, rsun_ref, rsun_obs))
                {
	        totaljz = DRMS_MISSING_FLOAT;
                totaljz_err = DRMS_MISSING_FLOAT;
	        }

            // Set sharp keys 	
            drms_setkey_float(sharpoutrec, "SAVNCPP",  totaljz); 
            drms_setkey_float(sharpoutrec, "ERRJHT",  totaljz_err); 	

            // Set sharp cea keys 
            drms_setkey_float(sharpceaoutrec, "SAVNCPP", totaljz); 
            drms_setkey_float(sharpceaoutrec, "ERRJHT", totaljz_err); 	
        }

        /***** MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) **********/
        if (meanshrflag || shrgt45flag)
        {
            // First compute potential field
	        for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
	        greenpot(bpx, bpy, bpz, nx, ny); 

            // Then compute shear angles
	        if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims, 
		   		                  &meanshear_angle, &meanshear_angle_err, &area_w_shear_gt_45,  
				                  mask, bitmask)) 
                {
	        meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
		area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
                meanshear_angle_err= DRMS_MISSING_FLOAT;
	        }

            // Set sharp keys 	
            drms_setkey_float(sharpoutrec, "MEANSHR",  meanshear_angle); 
            drms_setkey_float(sharpoutrec, "SHRGT45",  area_w_shear_gt_45); 	
            drms_setkey_float(sharpoutrec, "ERRMSHA",  meanshear_angle_err); 	

            // Set sharp cea keys 
            drms_setkey_float(sharpceaoutrec, "MEANSHR", meanshear_angle); 
            drms_setkey_float(sharpceaoutrec, "SHRGT45", area_w_shear_gt_45); 	
            drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err); 	  
        }

        /***** TOTFX, TOTFY, TOTFX, TOTBSQ, EPSX, EPSY, EPSZ , Example Function 16 (Lorentz forces) **********/
        if (totfxflag || totfyflag || totfzflag || totbsqflag || epsxflag || epsyflag || epszflag)
        {
            // Compute Lorentz forces
        	if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &totfx, &totfy, &totfz, &totbsq,
                               &epsx, &epsy, &epsz, mask, bitmask, cdelt1, rsun_ref, rsun_obs))
            {  
		totfx = DRMS_MISSING_FLOAT;
                totfy = DRMS_MISSING_FLOAT;
		totfz = DRMS_MISSING_FLOAT;
		totbsq = DRMS_MISSING_FLOAT;
                epsx = DRMS_MISSING_FLOAT;
		epsy = DRMS_MISSING_FLOAT;
                epsz = DRMS_MISSING_FLOAT;
            }

            // Set sharp keys 	
            drms_setkey_float(sharpoutrec, "TOTFX", totfx); 
            drms_setkey_float(sharpoutrec, "TOTFY", totfy); 	
            drms_setkey_float(sharpoutrec, "TOTFZ", totfz); 	
            drms_setkey_float(sharpoutrec, "TOTBSQ",totbsq); 
            drms_setkey_float(sharpoutrec, "EPSX",  epsx); 	
            drms_setkey_float(sharpoutrec, "EPSY",  epsy); 	
            drms_setkey_float(sharpoutrec, "EPSZ",  epsz); 	

            // Set sharp cea keys 
            drms_setkey_float(sharpceaoutrec, "TOTFX", totfx); 
            drms_setkey_float(sharpceaoutrec, "TOTFY", totfy); 	
            drms_setkey_float(sharpceaoutrec, "TOTFZ", totfz); 	
            drms_setkey_float(sharpceaoutrec, "TOTBSQ",totbsq); 
            drms_setkey_float(sharpceaoutrec, "EPSX",  epsx); 	
            drms_setkey_float(sharpceaoutrec, "EPSY",  epsy); 	
            drms_setkey_float(sharpceaoutrec, "EPSZ",  epsz); 	
              
        }

        /***** USFLUX, Example Function 17 *************************************/ 
	    if (meanvflosflag)
	    {
               //printf("772");
               if (computeAbsFlux_los(los, dims, &absFlux_los, &mean_vf_los,
                           &count_mask_los, bitmask, cdelt1, rsun_ref, rsun_obs))
               {
                absFlux_los = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
                mean_vf_los = DRMS_MISSING_FLOAT;
                count_mask_los = DRMS_MISSING_INT;
               }
            //printf("line757\n");
            //printf("mean_vf_los=%f\n",mean_vf_los);
            //printf("cdelt1=%f\n",cdelt1); 
            //printf("rsun_ref=%f\n",rsun_ref);
            //printf("rsun_obs=%f\n",rsun_obs);
            //printf("Key #%d done.\n", irec); 
            drms_setkey_float(sharpoutrec, "USFLUXL",  mean_vf_los);
            drms_setkey_float(sharpoutrec, "CMASKL",   count_mask_los); 	

            drms_setkey_float(sharpceaoutrec, "USFLUXL",  mean_vf_los);
            drms_setkey_float(sharpceaoutrec, "CMASKL",   count_mask_los); 	
        }

        /***** MEANGBZ, Example Function 18 ************************************/
        if (meanglosflag)
        {
        //printf("line 794\n");
            // Compute Bz derivative
            if (computeLOSderivative(los, dims, &mean_derivative_los, bitmask, derx_los, dery_los))
            {
                //printf("line 799\n");
                mean_derivative_los     = DRMS_MISSING_FLOAT; 
            }
	
            // Set sharp keys
            //printf("mean_derivative_los=%f\n",mean_derivative_los);
            drms_setkey_float(sharpoutrec, "MEANGBL",  mean_derivative_los); 
	
            // Set sharp cea keys   
            drms_setkey_float(sharpceaoutrec, "MEANGBL", mean_derivative_los); 
        }

        /******************************* END FLAGS **********************************/
        //printf("line809\n");
	drms_free_array(bitmaskArray);		
	drms_free_array(maskArray);
	drms_free_array(bxArray);           
	drms_free_array(byArray);
	drms_free_array(bzArray);
	drms_free_array(bx_errArray);           
	drms_free_array(by_errArray);
	drms_free_array(bz_errArray);
	drms_free_array(losArray);
	free(cvsinfoall);
    free(rim);
    free(p1p0);
    free(p1n0);
    free(p1p);
    free(p1n);
    free(p1);
    free(pmap);
    free(p1pad); 
    free(pmapn);
    free(fx); free(fy); free(fz);
    free(bh); free(bt); free(jz);
    free(bpx); free(bpy); free(bpz);
    free(derx); free(dery);
    free(derx_bt); free(dery_bt);
    free(derx_los); free(dery_los);
    free(derx_bz); free(dery_bz);
    free(derx_bh); free(dery_bh);
    free(bt_err); free(bh_err);  free(jz_err);
    free(jz_err_squared); free(jz_rms_err);
    free(jz_err_squared_smooth);
    free(jz_smooth);
    free(err_term2);
    free(err_term1);

	} //endfor

    // Close all the records
 	drms_close_records(sharpinrecset, DRMS_FREE_RECORD);
	drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);

	drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);
	drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);

	printf("Success=%d\n",sameflag);

return 0;    

}	// DoIt

Karen Tian
Powered by
ViewCVS 0.9.4