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

version 1.5, 2012/10/23 19:35:37

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

Line 20

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

Line 36

* 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 149 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 247 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 294 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 405 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 646 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 730 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 748 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 814 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 866 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 1107 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 1528 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 1589 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 1611 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 1690 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 1702 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, bitmask, cdelt1, rsun_ref, rsun_obs)){

Line 1729 void computeSWIndex(struct swIndex *swKe

Line 1750 void computeSWIndex(struct swIndex *swKe

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

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, bitmask, 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, bitmask, cdelt1, rsun_ref, rsun_obs)) {

swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;

Line 1750 void computeSWIndex(struct swIndex *swKe

Line 1776 void computeSWIndex(struct swIndex *swKe

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, bitmask, cdelt1, rsun_ref, rsun_obs))

swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;

Line 1774 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 1820 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;

int status = DRMS_SUCCESS;

status = 0;

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

drms_setkey_double(outRec, "DATE_B", val);

Line 1838 void setKeys(DRMS_Record_t *outRec, DRMS

Line 1903 void setKeys(DRMS_Record_t *outRec, DRMS

drms_setkey_string(outRec, "DATE", timebuf);

// set cvs commit version into keyword HEADER

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

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

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

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

Legend:

Removed from v.1.5
changed lines
	Added in v.1.11

Karen Tian