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.2 and 1.11

version 1.2, 2012/08/27 19:56:00

version 1.11, 2013/01/24 00:05:58

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

* v0.5 Jam 23 2012

* Notes:

* v0.0

Line 31

Line 35

* 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

* v0.5

* Corrected ephemeris keywords, added argument mInfo for setKeys()

* Example:

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

Line 147 enum projection {

Line 158 enum projection {

lambert

};

// WSC code

char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",

"SIN", "ZEA"};

// Ephemeris

struct ephemeris {

double disk_lonc, disk_latc;

Line 245 void computeSWIndex(struct swIndex *swKe

Line 260 void computeSWIndex(struct swIndex *swKe

void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);

/* Set all keywords, no error checking for now */

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo);

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

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

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

char *module_name = "sharp";

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

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

ModuleArgs_t module_args[] =

{

Line 403 int DoIt(void)

Line 418 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 582 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 608 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")) return 1;

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

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")) {

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

return 1;

}

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 617 int createCeaRecord(DRMS_Record_t *mharp

Line 662 int createCeaRecord(DRMS_Record_t *mharp

DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");

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

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

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

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

// Space weather

Line 700 int mapScaler(DRMS_Record_t *sharpRec, D

Line 746 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 718 int mapScaler(DRMS_Record_t *sharpRec, D

Line 764 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 784 int mapVectorB(DRMS_Record_t *sharpRec,

Line 831 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 836 int mapVectorBErr(DRMS_Record_t *sharpRe

Line 883 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 1077 int performSampling(float *outData, floa

Line 1124 int performSampling(float *outData, floa

// Rebinning, smoothing

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

free(outData0); // Dec 18 2012

Line 1454 int createCutRecord(DRMS_Record_t *mharp

Line 1502 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 1495 int createCutRecord(DRMS_Record_t *mharp

Line 1546 int createCutRecord(DRMS_Record_t *mharp

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

setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices

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

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

// Stats

Line 1556 int writeCutout(DRMS_Record_t *outRec, D

Line 1607 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 1578 int writeCutout(DRMS_Record_t *outRec, D

Line 1631 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 1664 void computeSWIndex(struct swIndex *swKe

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

Line 1657 void computeSWIndex(struct swIndex *swKe

Line 1710 void computeSWIndex(struct swIndex *swKe

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 1669 void computeSWIndex(struct swIndex *swKe

Line 1723 void computeSWIndex(struct swIndex *swKe

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

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

// Compute

// The Computations

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

mask, cdelt1, rsun_ref, rsun_obs)){

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;

}

Line 1680 void computeSWIndex(struct swIndex *swKe

Line 1734 void computeSWIndex(struct swIndex *swKe

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, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);

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

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

swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;

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

computeB_total(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))

swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;

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

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

swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;

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

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

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

computeJz(bx, by, dims, jz, mask, bitmask, cdelt1, rsun_ref, rsun_obs, derx, dery);

struct fresize_struct convolution_array;

init_fresize_gaussian(&convolution_array, 4.0f, 12, 1);

fresize(&convolution_array, jz, jz_smooth, nx, ny, nx, nx, ny, nx, 0, 0, 0.0f);

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

if(computeJzsmooth(bx, by, dims, jz_smooth, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask,

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

cdelt1, rsun_ref, rsun_obs, derx, dery))

{

swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;

swKeys_ptr->us_i = DRMS_MISSING_FLOAT;

}

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

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

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

swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;

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

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

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

mask, cdelt1, rsun_ref, rsun_obs)) {

mask, bitmask, 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;

}

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

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

mask, cdelt1, rsun_ref, rsun_obs))

mask, bitmask, cdelt1, rsun_ref, rsun_obs))

swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;

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

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

mask, cdelt1, rsun_ref, rsun_obs)) {

mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {

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

swKeys_ptr->totpot = DRMS_MISSING_FLOAT;

}

Line 1732 void computeSWIndex(struct swIndex *swKe

Line 1791 void computeSWIndex(struct swIndex *swKe

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

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

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

mask)) {

mask, bitmask)) {

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

Line 1741 void computeSWIndex(struct swIndex *swKe

Line 1800 void computeSWIndex(struct swIndex *swKe

// Clean up

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_fresize(&convolution_array);

}

Line 1787 void setSWIndex(DRMS_Record_t *outRec, s

Line 1848 void setSWIndex(DRMS_Record_t *outRec, s

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo)

{

copy_me_keys(inRec, outRec);

copy_patch_keys(inRec, outRec);

copy_geo_keys(inRec, outRec);

copy_ambig_keys(inRec, outRec);

int status = 0;

// Change a few geometry keywords for CEA records

if (mInfo != NULL) {

drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);

drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);

drms_setkey_float(outRec, "CRVAL1", mInfo->xc);

drms_setkey_float(outRec, "CRVAL2", mInfo->yc);

drms_setkey_float(outRec, "CDELT1", mInfo->xscale);

drms_setkey_float(outRec, "CDELT2", mInfo->yscale);

drms_setkey_string(outRec, "CUNIT1", "degree");

drms_setkey_string(outRec, "CUNIT2", "degree");

char key[64];

snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);

drms_setkey_string(outRec, "CTYPE1", key);

snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);

drms_setkey_string(outRec, "CTYPE2", key);

drms_setkey_float(outRec, "CROTA2", 0.0);

} else {

float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status);

float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status);

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

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

// Defined as disk center's pixel address wrt lower-left of cutout

drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);

drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);

// Always 0.

drms_setkey_float(outRec, "CRVAL1", 0);

drms_setkey_float(outRec, "CRVAL2", 0);

}

char timebuf[1024];

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

double val;

status = 0;

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.2
changed lines
	Added in v.1.11

Karen Tian