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

Return to sharp.c CVS log

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

Diff for /JSOC/proj/sharp/apps/sharp.c between version 1.1 and 1.14

version 1.1, 2012/08/23 08:38:33

version 1.14, 2013/05/21 17:06:32

Line 19

* Version:

* v0.0 Jul 02 2012

* v0.1 Jul 23 2012

* v0.2 Sep 04 2012

* v0.3 Dec 18 2012

* v0.4 Jan 02 2013

* Notes:

* v0.0

Line 31

Line 34

* Fixed char I/O thanks to Art

* SW indices fixed

* Added doppler and continuum

* Added other keywords: HEADER (populated by cvs build version), DATE_B

* v0.3

* Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"

* v0.4

* Took out convert_inplace(). Was causing all the images to be int

* Example:

* sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \

Line 62

Line 71

#include "errorprop.c"

#include "sw_functions.c"

//#include <mkl.h> // Comment out mkl.h, which can only run on solar3

#include <mkl_blas.h>

#include <mkl_service.h>

#include <mkl_lapack.h>

Line 112

Line 122

// Space weather keywords

struct swIndex {

float mean_vf;

float count_mask;

float absFlux;

float mean_hf;

float mean_gamma;

Line 131 struct swIndex {

Line 142 struct swIndex {

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

};

// Mapping method

Line 147 enum projection {

Line 173 enum projection {

lambert

};

// Ephemeris

// Ephemeris information

struct ephemeris {

double disk_lonc, disk_latc;

double disk_xc, disk_yc;

Line 181 int getInputRS_aux(DRMS_RecordSet_t **in

Line 207 int getInputRS_aux(DRMS_RecordSet_t **in

/* Find record from record set with given T_rec */

int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);

// ===================

/* Create CEA record */

int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,

DRMS_Record_t *dopRec, DRMS_Record_t *contRec,

Line 289 char *CEASegs[] = {"magnetogram", "bitma

Line 313 char *CEASegs[] = {"magnetogram", "bitma

BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};

/* ========================================================================================================== */

char *module_name = "sharp";

char *version_id = "2012 Jul 02"; /* Version number */

char *version_id = "2012 Dec 18"; /* Version number */

int seed;

ModuleArgs_t module_args[] =

{

Line 302 ModuleArgs_t module_args[] =

Line 325 ModuleArgs_t module_args[] =

{ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},

{ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},

{ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},

{ARG_INT, "seed", "987654", "Seed for the random number generator."},

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

Line 326 int DoIt(void)

Line 352 int DoIt(void)

contQuery = (char *) params_get_str(&cmdparams, "cont");

sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");

sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");

seed = params_get_int(&cmdparams, "seed");

/* Get input data, check everything */

Line 403 int DoIt(void)

Line 432 int DoIt(void)

} // irec

drms_close_records(mharpRS, DRMS_FREE_RECORD);

drms_close_records(bharpRS, DRMS_FREE_RECORD);

drms_close_records(dopRS, DRMS_FREE_RECORD); // Dec 18 2012

drms_close_records(contRS, DRMS_FREE_RECORD); // Dec 18 2012

return 0;

Line 564 int createCeaRecord(DRMS_Record_t *mharp

Line 596 int createCeaRecord(DRMS_Record_t *mharp

// Get ephemeris

if (getEphemeris(mharpRec, &(mInfo.ephem))) return 1;

if (getEphemeris(mharpRec, &(mInfo.ephem))) {

SHOW("CEA: get ephemeris error\n");

return 1;

}

// Find position

if (findPosition(mharpRec, &mInfo)) return 1;

if (findPosition(mharpRec, &mInfo)) {

SHOW("CEA: find position error\n");

return 1;

}

// Create xi_out, zeta_out array in mInfo:

// Coordinates to sample in original full disk image

Line 584 int createCeaRecord(DRMS_Record_t *mharp

Line 622 int createCeaRecord(DRMS_Record_t *mharp

// Mapping single segment: Mharp, etc.

if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) return 1;

if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {

if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) return 1;

SHOW("CEA: mapping magnetogram error\n");

return 1;

}

if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {

SHOW("CEA: mapping bitmap error\n");

return 1;

}

printf("Magnetogram mapping done.\n");

if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {

if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) return 1;

SHOW("CEA: mapping dopplergram error\n");

return 1;

}

printf("Dopplergram mapping done.\n");

if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) return 1;

if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {

SHOW("CEA: mapping continuum error\n");

return 1;

}

printf("Intensitygram mapping done.\n");

//if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) return 1;

if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {

//printf("Conf disambig mapping done.\n");

SHOW("CEA: mapping conf_disambig error\n");

return 1;

if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "disambig")) return 1;

}

printf("disambig mapping done.\n");

printf("Conf disambig mapping done.\n");

// Mapping vector B

if (mapVectorB(sharpRec, bharpRec, &mInfo)) return 1;

if (mapVectorB(sharpRec, bharpRec, &mInfo)) {

SHOW("CEA: mapping vector B error\n");

return 1;

}

printf("Vector B mapping done.\n");

// Mapping vector B errors

if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) return 1;

if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) {

SHOW("CEA: mapping vector B uncertainty error\n");

return 1;

}

printf("Vector B error done.\n");

// Keywords & Links

Line 621 int createCeaRecord(DRMS_Record_t *mharp

Line 676 int createCeaRecord(DRMS_Record_t *mharp

if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);

setKeys(sharpRec, bharpRec); // Set all other keywords

drms_copykey(sharpRec, mharpRec, "QUALITY"); // copied from los records

// Space weather

Line 703 int mapScaler(DRMS_Record_t *sharpRec, D

Line 759 int mapScaler(DRMS_Record_t *sharpRec, D

outSeg = drms_segment_lookup(sharpRec, segName);

if (!outSeg) return 1;

DRMS_Type_t arrayType = outSeg->info->type;

// DRMS_Type_t arrayType = outSeg->info->type;

DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);

if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}

// convert to needed data type

drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);

// drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray); // Jan 02 2013

outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];

outArray->parent_segment = outSeg;

// outArray->parent_segment = outSeg;

outArray->israw = 0; // always compressed

outArray->bzero = outSeg->bzero;

outArray->bscale = outSeg->bscale;

Line 721 int mapScaler(DRMS_Record_t *sharpRec, D

Line 777 int mapScaler(DRMS_Record_t *sharpRec, D

if (status) return 0;

if (inArray) drms_free_array(inArray);

if ((xsz != FOURK) || (ysz != FOURK)) free(inData); // Dec 18 2012

if (outArray) drms_free_array(outArray);

return 0;

Line 787 int mapVectorB(DRMS_Record_t *sharpRec,

Line 844 int mapVectorB(DRMS_Record_t *sharpRec,

outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);

outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);

outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];

outArray->parent_segment = outSeg;

// outArray->parent_segment = outSeg;

outArray->israw = 0;

outArray->bzero = outSeg->bzero;

outArray->bscale = outSeg->bscale;

Line 839 int mapVectorBErr(DRMS_Record_t *sharpRe

Line 896 int mapVectorBErr(DRMS_Record_t *sharpRe

outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);

outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);

outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];

outArray->parent_segment = outSeg;

// outArray->parent_segment = outSeg;

outArray->israw = 0;

outArray->bzero = outSeg->bzero;

outArray->bscale = outSeg->bscale;

Line 884 int findPosition(DRMS_Record_t *inRec, s

Line 941 int findPosition(DRMS_Record_t *inRec, s

// We compute minlon & minlat then by

// LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT

float psize = drms_getkey_float(inRec, "SIZE", &status);

// float psize = drms_getkey_float(inRec, "SIZE", &status);

if (psize != psize) {

// if (psize != psize) {

TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1;

if (minlon != minlon || maxlon != maxlon) { // check lons instead of SIZE

TIME t0 = drms_getkey_time(inRec, "T_FRST1", &status); if (status) return 1; // changed from T_FRST to T_FRST1, T_FRST may not exist

double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;

char firstRecQuery[100], t0_str[100];

sprint_time(t0_str, t0, "TAI", 0);

Line 1080 int performSampling(float *outData, floa

Line 1139 int performSampling(float *outData, floa

// Rebinning, smoothing

frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian

free(outData0); // Dec 18 2012

Line 1457 int createCutRecord(DRMS_Record_t *mharp

Line 1517 int createCutRecord(DRMS_Record_t *mharp

break;

}

if (iHarpSeg != nMharpSegs) return 1; // if failed

if (iHarpSeg != nMharpSegs) {

SHOW("Cutout: segment number unmatch\n");

return 1; // if failed

}

printf("Magnetogram cutout done.\n");

// Cutout Doppler

Line 1559 int writeCutout(DRMS_Record_t *outRec, D

Line 1622 int writeCutout(DRMS_Record_t *outRec, D

#if DISAMB_AZI

if (!strcmp(SegName, "azimuth")) {

DRMS_Segment_t *disambSeg = drms_segment_lookup(inRec, "disambig");

DRMS_Segment_t *disambSeg = NULL;

disambSeg = drms_segment_lookup(inRec, "disambig");

if (!disambSeg) {drms_free_array(cutoutArray); return 1;}

DRMS_Array_t *disambArray;

if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {

disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);

if (status) {drms_free_array(cutoutArray); return 1;}

} else {

drms_free_array(cutoutArray);

return 1;

}

double *azimuth = (double *) cutoutArray->data;

Line 1581 int writeCutout(DRMS_Record_t *outRec, D

Line 1646 int writeCutout(DRMS_Record_t *outRec, D

outSeg = drms_segment_lookup(outRec, SegName);

if (!outSeg) return 1;

drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);

// drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray); // Jan 02 2013

outSeg->axis[0] = cutoutArray->axis[0];

outSeg->axis[1] = cutoutArray->axis[1];

cutoutArray->parent_segment = outSeg;

// cutoutArray->parent_segment = outSeg;

cutoutArray->israw = 0; // always compressed

cutoutArray->bzero = outSeg->bzero;

cutoutArray->bscale = outSeg->bscale; // Same as inArray's

Line 1611 void computeSWIndex(struct swIndex *swKe

Line 1676 void computeSWIndex(struct swIndex *swKe

int nx = mInfo->ncol, ny = mInfo->nrow;

int nxny = nx * ny;

int dims[2] = {nx, ny};

// Get bx, by, bz, mask

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

//DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "bitmap");

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

//DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);

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

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

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(inRec, "conf_disambig");

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

DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "disambig");

DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);

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

Line 1637 void computeSWIndex(struct swIndex *swKe

Line 1702 void computeSWIndex(struct swIndex *swKe

DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);

float *bz = (float *) bzArray->data; // bz

DRMS_Segment_t *bz_errSeg = drms_segment_lookup(inRec, BR_ERR_SEG_CEA);

DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);

float *bz_err = (float *) bz_errArray->data; // bz_err

DRMS_Segment_t *by_errSeg = drms_segment_lookup(inRec, BT_ERR_SEG_CEA);

DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);

float *by_err = (float *) by_errArray->data; // by_err

//for (int i = 0; i < nxny; i++) by_err[i] *= -1;

DRMS_Segment_t *bx_errSeg = drms_segment_lookup(inRec, BP_ERR_SEG_CEA);

DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);

float *bx_err = (float *) bx_errArray->data; // bx_err

// Get emphemeris

//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);

Line 1649 void computeSWIndex(struct swIndex *swKe

Line 1727 void computeSWIndex(struct swIndex *swKe

float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status);

float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status);

//float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)

printf("cdelt1=%f\n",cdelt1);

printf("rsun_ref=%f\n",rsun_ref);

printf("rsun_obs=%f\n",rsun_obs);

printf("dsun_obs=%f\n",dsun_obs);

// Temp arrays

float *bh = (float *) (malloc(nxny * sizeof(float)));

float *bt = (float *) (malloc(nxny * sizeof(float)));

float *jz = (float *) (malloc(nxny * sizeof(float)));

float *jz_smooth = (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)));

Line 1672 void computeSWIndex(struct swIndex *swKe

Line 1744 void computeSWIndex(struct swIndex *swKe

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_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));

float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));

//spaceweather quantities computed

// Compute

if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),

if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf), &(swKeys_ptr->mean_vf_err),

mask, cdelt1, rsun_ref, rsun_obs)){

&(swKeys_ptr->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs))

{

swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN

swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;

swKeys_ptr->mean_vf_err = DRMS_MISSING_FLOAT;

swKeys_ptr->count_mask = DRMS_MISSING_INT;

}

for (int i = 0; i < nxny; i++) bpz[i] = bz[i];

greenpot(bpx, bpy, bpz, nx, ny);

computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask);

computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);

if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask))

if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), &(swKeys_ptr->mean_gamma_err),mask, bitmask))

{

swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;

swKeys_ptr->mean_gamma_err = DRMS_MISSING_FLOAT;

}

computeB_total(bx, by, bz, bt, dims, mask);

computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);

if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt))

if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err)))

{

swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;

swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT;

}

if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh))

if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh), &(swKeys_ptr->mean_derivative_bh_err), mask, bitmask, derx_bh, dery_bh))

{

swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;

swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT;

}

if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, derx_bz, dery_bz))

if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err), mask, bitmask, derx_bz, dery_bz))

{

swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN

swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT;

}

computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,

derx, dery);

if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask,

if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &(swKeys_ptr->mean_jz),

cdelt1, rsun_ref, rsun_obs, derx, dery)) {

&(swKeys_ptr->mean_jz_err), &(swKeys_ptr->us_i), &(swKeys_ptr->us_i_err), mask, bitmask, cdelt1,

rsun_ref, rsun_obs, derx, dery))

{

swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;

swKeys_ptr->us_i = DRMS_MISSING_FLOAT;

swKeys_ptr->mean_jz_err = DRMS_MISSING_FLOAT;

swKeys_ptr->us_i_err = DRMS_MISSING_FLOAT;

}

printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);

if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &(swKeys_ptr->mean_alpha), &(swKeys_ptr->mean_alpha_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs))

{

if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, cdelt1, rsun_ref, rsun_obs))

swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;

swKeys_ptr->mean_alpha_err = DRMS_MISSING_FLOAT;

}

if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih),

if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &(swKeys_ptr->mean_ih), &(swKeys_ptr->mean_ih_err), &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),

&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),

&(swKeys_ptr->total_us_ih_err), &(swKeys_ptr->total_abs_ih_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs))

mask, cdelt1, rsun_ref, rsun_obs)) {

{

swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;

swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;

swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;

swKeys_ptr->mean_ih_err = DRMS_MISSING_FLOAT;

swKeys_ptr->total_us_ih_err = DRMS_MISSING_FLOAT;

swKeys_ptr->total_abs_ih_err = DRMS_MISSING_FLOAT;

}

if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),

if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &(swKeys_ptr->totaljz), &(swKeys_ptr->totaljz_err),

mask, cdelt1, rsun_ref, rsun_obs))

mask, bitmask, cdelt1, rsun_ref, rsun_obs))

{

swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;

swKeys_ptr->totaljz_err = DRMS_MISSING_FLOAT;

}

if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims,

if (computeFreeEnergy(bx, by, bpx, bpy, dims,

&(swKeys_ptr->meanpot), &(swKeys_ptr->meanpot_err), &(swKeys_ptr->totpot), &(swKeys_ptr->totpot_err),

&(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),

mask, bitmask, cdelt1, rsun_ref, rsun_obs))

mask, cdelt1, rsun_ref, rsun_obs)) {

{

swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN

swKeys_ptr->totpot = DRMS_MISSING_FLOAT;

swKeys_ptr->meanpot_err = DRMS_MISSING_FLOAT;

swKeys_ptr->totpot_err = DRMS_MISSING_FLOAT;

}

if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,

if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims,

&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),

&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),

&(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),

mask, bitmask)) {

mask)) {

swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN

swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;

swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN

swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;

swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;

}

// Clean up

// Clean up the arrays

drms_free_array(bitmaskArray); // Dec 18 2012 Xudong

drms_free_array(maskArray);

drms_free_array(bxArray);

drms_free_array(byArray);

drms_free_array(bzArray);

free(bh); free(bt); free(jz);

free(bh); free(bt); free(jz); free(jz_smooth);

free(bpx); free(bpy); free(bpz);

free(derx); free(dery);

free(derx_bt); free(dery_bt);

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

}

* Set space weather indices, no error checking for now

Line 1783 void setSWIndex(DRMS_Record_t *outRec, s

Line 1891 void setSWIndex(DRMS_Record_t *outRec, s

drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);

drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);

drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);

drms_setkey_float(outRec, "CMASK", swKeys_ptr->count_mask);

drms_setkey_float(outRec, "ERRBT", swKeys_ptr->mean_derivative_btotal_err);

drms_setkey_float(outRec, "ERRVF", swKeys_ptr->mean_vf_err);

drms_setkey_float(outRec, "ERRGAM", swKeys_ptr->mean_gamma_err);

drms_setkey_float(outRec, "ERRBH", swKeys_ptr->mean_derivative_bh_err);

drms_setkey_float(outRec, "ERRBZ", swKeys_ptr->mean_derivative_bz_err);

drms_setkey_float(outRec, "ERRJZ", swKeys_ptr->mean_jz_err);

drms_setkey_float(outRec, "ERRUSI", swKeys_ptr->us_i_err);

drms_setkey_float(outRec, "ERRALP", swKeys_ptr->mean_alpha_err);

drms_setkey_float(outRec, "ERRMIH", swKeys_ptr->mean_ih_err);

drms_setkey_float(outRec, "ERRTUI", swKeys_ptr->total_us_ih_err);

drms_setkey_float(outRec, "ERRTAI", swKeys_ptr->total_abs_ih_err);

drms_setkey_float(outRec, "ERRJHT", swKeys_ptr->totaljz_err);

drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);

drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);

drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);

};

* Set all keywords, no error checking for now

Line 1797 void setKeys(DRMS_Record_t *outRec, DRMS

Line 1920 void setKeys(DRMS_Record_t *outRec, DRMS

copy_patch_keys(inRec, outRec);

copy_geo_keys(inRec, outRec);

copy_ambig_keys(inRec, outRec);

char timebuf[1024];

float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */

double val;

int status = DRMS_SUCCESS;

val = drms_getkey_double(inRec, "DATE",&status);

drms_setkey_double(outRec, "DATE_B", val);

sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);

drms_setkey_string(outRec, "DATE", timebuf);

// set cvs commit version into keyword HEADER

char *cvsinfo = strdup("$Id$");

// status = drms_setkey_string(outRec, "HEADER", cvsinfo);

status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);

};

Legend:

Removed from v.1.1
changed lines
	Added in v.1.14

Karen Tian