JSOC/proj/sharp/apps/update_sharp_keys.c - diff

Return to update_sharp_keys.c CVS log

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

Diff for /JSOC/proj/sharp/apps/update_sharp_keys.c between version 1.2 and 1.18

version 1.2, 2013/08/06 23:26:42

version 1.18, 2021/05/26 04:43:52

Line 1

* MODULE NAME: update_sharp_keys.c

* DESCRIPTION: This module recalculates SHARP keywords. This is accomplished by

* DESCRIPTION: This module recalculates SHARP keywords.

* cloning a record, recalculating keywords of choice (i.e., user input),

* This is accomplished by cloning a record, recalculating keywords of choice

* and pointing to the same segments as the old record. Associated error keys

* (i.e., user input), and pointing to the same segments as the old record.

* are computed by default.

* 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

Line 12

* 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

* Version : v0.0 Jun 14 2013

* EXAMPLE :

* update_sharp_keys sharpseries=hmi.sharp_720s sharpceaseries=hmi.sharp_cea_720s //

* update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //

* HARPNUM=1 keylist=USFLUX,TOTPOT

* HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT debug=debug

Line 58

Line 63

#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 = "2013 Jun 14"; /* Version number */

char *version_id = "2021 May 25"; /* Version number */

ModuleArgs_t module_args[] =

{

{ARG_STRING, "sharpseries", NULL, "Input/output Sharp dataseries"},

{ARG_STRING, "sharpseriesin", NULL, "Input Sharp dataseries"},

{ARG_STRING, "sharpceaseries", NULL, "Input/output Sharp CEA 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"},

{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}

};

Line 82 int DoIt(void)

Line 91 int DoIt(void)

/* 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;

Line 116 int DoIt(void)

Line 129 int DoIt(void)

float totaljz_err;

float meanpot_err;

float totpot_err;

float Rparam;

float meanshear_angle_err;

char *sharpseries = (char *) params_get_str(&cmdparams, "sharpseries");

float totfx;

char *sharpceaseries = (char *) params_get_str(&cmdparams, "sharpceaseries");

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(sharpceaquery, "%s[%d]\n", sharpceaseries,harpnum);

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

printf(sharpceaquery, "%s[%d]\n", sharpceaseries,harpnum);

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

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

printf(sharpquery, "%s[%d]\n", sharpseries,harpnum);

DRMS_RecordSet_t *sharpinrecset = drms_open_records(drms_env, sharpquery, &status);

if (status != DRMS_SUCCESS)

DIE("Problem opening sharp recordset.");

if (sharpinrecset->n == 0)

nrecs=sharpinrecset->n;

if (nrecs == 0)

DIE("Sharp recordset contains no records.");

DRMS_RecordSet_t *sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);

printf("nrecs=%d\n",nrecs);

if (status != DRMS_SUCCESS)

DIE("Problem cloning sharp records.");

DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status);

if (status != DRMS_SUCCESS)

DIE("Problem opening sharp cea recordset.");

if (sharpceainrecset->n == 0)

if (sharpceainrecset->n != nrecs)

DIE("Sharp cea recordset contains no records.");

DIE("Sharp cea recordset differs from sharp recordset.");

DRMS_RecordSet_t *sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);

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 meanvfflag = (strstr(keylist,"USFLUX") != NULL); // generalize so that lowercase is acceptable

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);

Line 167 int DoIt(void)

Line 216 int DoIt(void)

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);

DRMS_Record_t *sharpinrec = sharpinrecset->records[0];

DRMS_Record_t *sharpceainrec = sharpceainrecset->records[0];

DRMS_Segment_t *inseg = drms_segment_lookup(sharpceainrec, "Br");

Line 174 int DoIt(void)

Line 233 int DoIt(void)

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)));

Line 197 int DoIt(void)

Line 255 int DoIt(void)

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)));

for (irec=0;irec<nrecs;irec++)

//for (irec=5;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");

// ephemeris variables

float cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2;

double rsun_ref, rsun_obs;

// outrecords

DRMS_Record_t *sharpoutrec, *sharpceaoutrec;

DRMS_Record_t *sharpoutrec, *sharpceaoutrec, *testrec;

nrecs=sharpoutrecset->n;

printf("nrecs=%d\n",nrecs);

// prepare to set CODEVER7 (CVS Version of the SHARP module)

char *cvsinfo0;

Line 224 int DoIt(void)

Line 295 int DoIt(void)

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);

//printf("historyofthemodule=%s\n",historyofthemodule);

history0 = drms_getkey_string(sharpinrec, "HISTORY", &status);

strcat(historyofthemodule,"\n");

strcat(historyofthemodule,history0);

// no longer need inrecs

// grab one record such that we can malloc the R-parameter calculation arrays outside of the loop

drms_close_records(sharpinrecset, DRMS_FREE_RECORD);

// assume scale = round(2./cdelt1) does not change for any given HARP

drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);

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));

for (irec=0;irec<nrecs;irec++)

{

// 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];

cdelt1_orig = drms_getkey_float(sharpceaoutrec, "CDELT1", &status);

if (!sameflag)

dsun_obs = drms_getkey_float(sharpoutrec, "DSUN_OBS", &status);

{

rsun_ref = drms_getkey_double(sharpoutrec, "RSUN_REF", &status);

drms_copykeys(sharpoutrec, sharpinrec, 0, kDRMS_KeyClass_Explicit);

rsun_obs = drms_getkey_double(sharpoutrec, "RSUN_OBS", &status);

drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit);

imcrpix1 = drms_getkey_float(sharpoutrec, "IMCRPIX1", &status);

}

imcrpix2 = drms_getkey_float(sharpoutrec, "IMCRPIX2", &status);

crpix1 = drms_getkey_float(sharpoutrec, "CRPIX1", &status);

TIME trec, tnow; /* 1970.01.01_00:00:00_UTC */

crpix2 = drms_getkey_float(sharpoutrec, "CRPIX2", &status);

tnow = (double)time(NULL);

cdelt1 = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);

tnow += UNIX_epoch;

if (debugflag)

{

printf("i=%d\n",irec);

}

// set CODEVER7 and HISTORY

// 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

// Get bx, by, bz, mask

DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");

// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities

DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceaoutrec, "bitmap");

DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);

int *bitmask = (int *) bitmaskArray->data; // get the previously made mask array

//Use conf_disambig map as a threshold on spaceweather quantities

DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");

DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceaoutrec, "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 *bxSeg = drms_segment_lookup(sharpceaoutrec, "Bp");

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(sharpceaoutrec, "Bt");

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_Segment_t *bzSeg = drms_segment_lookup(sharpceaoutrec, "Br");

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(sharpceaoutrec, "Br_err");

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(sharpceaoutrec, "Bt_err");

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(sharpceaoutrec, "Bp_err");

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 ************************************sdfdsf*/

/***** USFLUX, Example Function 1 *************************************/

if (meanvfflag)

{

// Compute unsigned flux

Line 313 int DoIt(void)

Line 453 int DoIt(void)

drms_setkey_float(sharpceaoutrec, "ERRVF", mean_vf_err);

}

/***** MEANPOT and TOTPOT, Example Function 13 (Requires Keiji's code) **/

/***** 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

Line 346 int DoIt(void)

Line 518 int DoIt(void)

}

/***** MEANGAM, Example Function 3 ************************************/

if (meangamflag)

{

// First compute horizontal field

Line 369 int DoIt(void)

Line 540 int DoIt(void)

}

/***** MEANGBT, Example Function 5 (Requires Function 4) *************/

if (meangbtflag)

{

// First, compute Bt

Line 377 int DoIt(void)

Line 547 int DoIt(void)

// 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))

&mean_derivative_btotal_err, err_term1, err_term2))

{

mean_derivative_btotal = DRMS_MISSING_FLOAT;

mean_derivative_btotal_err = DRMS_MISSING_FLOAT;

Line 393 int DoIt(void)

Line 563 int DoIt(void)

}

/***** 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))

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;

Line 415 int DoIt(void)

Line 584 int DoIt(void)

drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err);

}

/***** 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))

{

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);

}

/***** 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);

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;

Line 470 int DoIt(void)

Line 617 int DoIt(void)

}

/***** 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);

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,

Line 496 int DoIt(void)

Line 642 int DoIt(void)

}

/***** 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);

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,

Line 533 int DoIt(void)

Line 678 int DoIt(void)

}

/***** 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);

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,

Line 558 int DoIt(void)

Line 702 int DoIt(void)

}

/***** MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) **********/

if (meanshrflag || shrgt45flag)

{

// First compute potential field

Line 586 int DoIt(void)

Line 729 int DoIt(void)

drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err);

}

/******************************* END FLAGS **********************************/

/***** 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);

Line 596 int DoIt(void)

Line 818 int DoIt(void)

drms_free_array(bx_errArray);

drms_free_array(by_errArray);

drms_free_array(bz_errArray);

} //endfor

drms_free_array(losArray);

free(cvsinfoall);

free(rim);

free(p1p0);

free(p1n0);

free(p1p);

free(p1n);

free(p1);

free(pmap);

free(p1pad);

free(pmapn);

drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);

} //endfor

drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);

//used for select calculations

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(cvsinfoall);

free(jz_err_squared_smooth);

free(jz_smooth);

free(err_term2);

free(err_term1);

return 1;

} // DoIt

// 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

Legend:

Removed from v.1.2
changed lines
	Added in v.1.18

Karen Tian