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.21 and 1.22

version 1.21, 2014/02/18 23:35:31

version 1.22, 2014/02/19 14:59:25

Line 1

* sharp.c

* This module creates the pipeline Space Weather Active Region Patches (SHARPs).

* This module creates the pipeline space weather harps

* It is a hard-coded strip-down version of bmap.c.

* It is a hard-coded strip-down version of bmap.c

* It takes the Mharp and Bharp series and create the following quantities:

* It takes the Mharp and Bharp series and crete the following quantities

* Series 1: Sharp_CEA

* CEA remapped magnetogram, bitmap, continuum, doppler (same size in map coordinate, need manual spec?)

* CEA remapped vector field (Br, Bt, Bp) (same as above)

* Space weather indices based on vector cutouts (step 2)

* Series 2: Sharp_cutout:

* cutouts of magnetogram, bitmap, continuum, doppler (HARP defined, various sizes in CCD pixels)

* cutouts of all vector data segments (same as above)

* Series 3: Other remaps

* Author:

* Xudong Sun; Monica Bobra

Line 25

Line 24

* v0.4 Jan 02 2013

* v0.5 Jan 23 2013

* v0.6 Aug 12 2013

* v0.7 Jan 02 2014

* v0.8 Feb 12 2014

* Notes:

* v0.0

Line 48

Line 45

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

* v0.6

* Changes in remapping of bitmap and conf_disambig, now near neighbor without anti-aliasing

* v0.7

* Added full disk as "b"

* Global flag fullDisk is set if "b" is set

* Example:

* Utilize BharpRS and BharpRec all around

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

* Pass mharpRec to set_keys() too in case of full disk

"bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \

* Fixed Bunit (removed from copy_me_keys(), added loops for Bunits in set_keys() here)

"dop=hmi.V_720s[2012.02.20_10:00]" \

* Error for CEA still does account for disambiguation yet

"cont=hmi.Ic_720s[2012.02.20_10:00]" \

* v0.8

"sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"

* Added disambig to azimuth during error propagation

* For comparison:

* Changed usage for disambig: bit 2 (radial acute) for full disk, bit 0 for patch

* bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \

* Fixed disambig cutout for patch: 0 for even, 7 for odd

"out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"

* Example Calls:

* [I] B (full disk disambiguation)

* > sharp "mharp=hmi.Mharp_720s[1832][2012.07.12_15:24]" "b=hmi_test.B_720s[2012.07.12_15:24]" "dop=hmi.V_720s[2012.07.12_15:24]" "cont=hmi.Ic_720s[2012.07.12_15:24]" "sharp_cea=su_xudong.sharp_cea_720s" "sharp_cut=su_xudong.sharp_720s"

* [II] BHARP (patch disambiguation)

* > sharp "mharp=hmi.Mharp_720s[1832][2012.07.12_15:24]" "bharp=hmi.Bharp_720s[1832][2012.07.12_15:24]" "dop=hmi.V_720s[2012.07.12_15:24]" "cont=hmi.Ic_720s[2012.07.12_15:24]" "sharp_cea=su_xudong.sharp_cea_720s" "sharp_cut=su_xudong.sharp_720s"

Line 175 struct swIndex {

Line 165 struct swIndex {

float meanpot_err;

float totpot_err;

float meanshear_angle_err;

float Rparam;

};

// Mapping method

Line 291 void computeSWIndex(struct swIndex *swKe

Line 280 void computeSWIndex(struct swIndex *swKe

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

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

// Changed Dec 30 XS

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

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

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

Line 332 char *CutSegs[] = {"magnetogram", "bitma

Line 320 char *CutSegs[] = {"magnetogram", "bitma

"field_inclination_err", "field_az_err", "inclin_azimuth_err",

"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",

"disambig", "conf_disambig"};

char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",

char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",

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

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

// For Bunits, added Dec 30 XS

char *CutBunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",

"degree", "degree", "Mx/cm^2", "cm/s", "mA", " ",

"length units", "DN/s", "DN/s", " ", " ",

" ",

" ", " ",

"degree", "degree", "Mx/cm^2", "cm/s", " ",

" ", " ", " ",

" ", " "};

char *CEABunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s", " ",

"Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2"};

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

char *module_name = "sharp";

int seed;

int fullDisk; // full disk mode

ModuleArgs_t module_args[] =

{

{ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},

{ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},

{ARG_STRING, "b", kNotSpecified, "Input B series, if set, overrides bharp."},

{ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},

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

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

Line 375 int DoIt(void)

Line 348 int DoIt(void)

int status = DRMS_SUCCESS;

int nrecs, irec;

char *mharpQuery, *bharpQuery, *bQuery;

char *mharpQuery, *bharpQuery;

char *dopQuery, *contQuery;

char *sharpCeaQuery, *sharpCutQuery;

Line 386 int DoIt(void)

Line 359 int DoIt(void)

mharpQuery = (char *) params_get_str(&cmdparams, "mharp");

bharpQuery = (char *) params_get_str(&cmdparams, "bharp");

bQuery = (char *) params_get_str(&cmdparams, "b");

dopQuery = (char *) params_get_str(&cmdparams, "dop");

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

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

Line 396 int DoIt(void)

Line 368 int DoIt(void)

/* Get input data, check everything */

// Full disk mode if "b" is set

if (strcmp(bQuery, kNotSpecified)) {

fullDisk = 1;

bharpQuery = bQuery;

// SHOW(bharpQuery); SHOW("\n");

SHOW("Full disk mode\n");

} else {

fullDisk = 0;

SHOW("Harp mode\n");

}

// Bharp point to B if full disk

if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))

DIE("Input harp data error.");

nrecs = mharpRS->n;

Line 427 int DoIt(void)

Line 387 int DoIt(void)

/* Records in work */

DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;

mharpRec = mharpRS->records[irec];

bharpRec = bharpRS->records[irec];

TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);

if (!fullDisk) {

bharpRec = bharpRS->records[irec];

} else {

if (getInputRec_aux(&bharpRec, bharpRS, trec)) { // Bharp point to full disk B

printf("Fetching B failed, image #%d skipped.\n", irec);

continue;

}

struct swIndex swKeys;

DRMS_Record_t *dopRec = NULL, *contRec = NULL;

if (getInputRec_aux(&dopRec, dopRS, trec)) {

printf("Fetching Doppler failed, image #%d skipped.\n", irec);

continue;

Line 523 int getInputRS(DRMS_RecordSet_t **mharpR

Line 473 int getInputRS(DRMS_RecordSet_t **mharpR

*mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);

if (status || (*mharpRS_ptr)->n == 0) return 1;

if (fullDisk) {

if (getInputRS_aux(bharpRS_ptr, bharpQuery, *mharpRS_ptr)) return 1;

} else {

*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);

if (status || (*bharpRS_ptr)->n == 0) return 1;

if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;

}

return 0;

Line 756 int createCeaRecord(DRMS_Record_t *mharp

Line 703 int createCeaRecord(DRMS_Record_t *mharp

drms_copykey(sharpRec, mharpRec, "T_REC");

drms_copykey(sharpRec, mharpRec, "HARPNUM");

if (fullDisk) {

DRMS_Link_t *bLink = hcon_lookup_lower(&sharpRec->links, "B");

if (bLink) drms_link_set("B", sharpRec, bharpRec);

} else {

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

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

}

DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");

if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);

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

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

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

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

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

// Space weather

Line 820 int mapScaler(DRMS_Record_t *sharpRec, D

Line 762 int mapScaler(DRMS_Record_t *sharpRec, D

inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);

if (!inArray) return 1;

if (!strcmp(segName, "conf_disambig") || !strcmp(segName, "bitmap")) {

// Moved out so it works for FD conf_disambig as well

// Jan 2 2014 XS

interpOpt = 3; // Aug 12 XS, near neighbor

}

float *inData;

int xsz = inArray->axis[0], ysz = inArray->axis[1];

if ((xsz != FOURK) || (ysz != FOURK)) { // for bitmap, make tmp full disk

interpOpt = 3; // Aug 12 XS, near neighbor

float *inData0 = (float *) inArray->data;

inData = (float *) (calloc(FOURK2, sizeof(float)));

int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;

Line 1498 int readVectorB(DRMS_Record_t *inRec, fl

Line 1435 int readVectorB(DRMS_Record_t *inRec, fl

int llx, lly; // lower-left corner

int bmx, bmy; // bitmap size

if (fullDisk) {

llx = lly = 0;

bmx = bmy = FOURK;

} else {

llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;

lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;

bmx = inArray_ambig->axis[0];

bmy = inArray_ambig->axis[1];

}

int kx, ky, kOff;

int ix = 0, jy = 0, yOff = 0, iData = 0;

int xDim = FOURK, yDim = FOURK;

int amb = 0;

for (jy = 0; jy < yDim; jy++)

{

Line 1534 int readVectorB(DRMS_Record_t *inRec, fl

Line 1466 int readVectorB(DRMS_Record_t *inRec, fl

continue;

} else {

kOff = ky * bmx + kx;

// if (ambig[kOff] % 2) { // 180

if (ambig[kOff] % 2) { // 180

// Feb 12 2014, use bit #2 for full disk, lowest bit for patch

if (fullDisk) { amb = (ambig[kOff] / 4) % 2; } else { amb = ambig[kOff] % 2; }

if (amb) { // Feb 12 2014, use bit #2

bx_img[iData] *= -1.; by_img[iData] *= -1.;

}

Line 1577 int readVectorBErr(DRMS_Record_t *inRec,

Line 1506 int readVectorBErr(DRMS_Record_t *inRec,

DRMS_Array_t *inArrays[9];

// Read full disk images

// Do we need disambig? Dec 30 XS

for (int iSeg = 0; iSeg < 9; iSeg++) {

Line 1591 int readVectorBErr(DRMS_Record_t *inRec,

Line 1519 int readVectorBErr(DRMS_Record_t *inRec,

float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];

float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];

// Add disambig, Feb 12 2014

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

DRMS_Array_t *inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);

if (status) return 1;

char *ambig = (char *)inArray_ambig->data;

int llx, lly; // lower-left corner

int bmx, bmy; // bitmap size

if (fullDisk) {

llx = lly = 0;

bmx = bmy = FOURK;

} else {

llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;

lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;

bmx = inArray_ambig->axis[0];

bmy = inArray_ambig->axis[1];

}

int idx, idx_a;

int amb;

for (int j = 0; j < bmy; j++) {

for (int i = 0; i < bmx; i++) {

idx_a = j * bmx + i;

idx = (j + lly) * FOURK + (i + llx);

// Feb 12 2014, use bit #2 for full disk, lowest bit for patch

if (fullDisk) { amb = (ambig[idx_a] / 4) % 2; } else { amb = ambig[idx_a] % 2; }

if (amb) { bA0[idx] += 180.; }

}

// Convert errors to variances, correlation coefficients to covariances

for (int i = 0; i < FOURK2; i++) {

Line 1632 int readVectorBErr(DRMS_Record_t *inRec,

Line 1527 int readVectorBErr(DRMS_Record_t *inRec,

if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;

bT[i] = bT0[i];

bI[i] = bI0[i]; // in deg, coverted in errorprop

bI[i] = bI0[i];

bA[i] = bA0[i];

errbT[i] = errbT0[i] * errbT0[i];

Line 1648 int readVectorBErr(DRMS_Record_t *inRec,

Line 1543 int readVectorBErr(DRMS_Record_t *inRec,

for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);

drms_free_array(inArray_ambig); // Feb 12 2014

return 0;

Line 1724 int createCutRecord(DRMS_Record_t *mharp

Line 1618 int createCutRecord(DRMS_Record_t *mharp

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

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

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

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

// Stats

Line 1781 int writeCutout(DRMS_Record_t *outRec, D

Line 1675 int writeCutout(DRMS_Record_t *outRec, D

return 1;

}

// Feb 12 2014, fool-proof, for patch, change everything to 0 or 7!!!

// This is a fix for disambiguation before Aug 2013

if (!strcmp(SegName, "disambig") && !fullDisk) {

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

for (int i = 0; i < nxny; i++) {

if (((int)disamb[i]) % 2) { disamb[i] = 7; } else { disamb[i] = 0; }

}

/* Adding disambiguation resolution to cutout azimuth? */

#if DISAMB_AZI

int amb;

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

DRMS_Segment_t *disambSeg = NULL;

disambSeg = drms_segment_lookup(inRec, "disambig");

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

DRMS_Array_t *disambArray;

if (fullDisk) { // Jan 2 2014 XS

disambArray = drms_segment_readslice(disambSeg, DRMS_TYPE_CHAR, ll, ur, &status);

if (status) return 1;

} else {

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

Line 1811 int writeCutout(DRMS_Record_t *outRec, D

Line 1690 int writeCutout(DRMS_Record_t *outRec, D

drms_free_array(cutoutArray);

return 1;

}

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

char *disamb = (char *) disambArray->data;

for (int n = 0; n < nxny; n++) {

// if (disamb[n] % 2) azimuth[n] += 180.; // Nov 12 2013 Fixed!!!

if (disamb[n] % 2) azimuth[n] += 180.; // Nov 12 2013 Fixed!!!

// Feb 12 2014, use bit #2 for full disk, lowest bit for patch

if (fullDisk) { amb = (disamb[n] / 4) % 2; } else { amb = disamb[n] % 2; }

if (amb) azimuth[n] += 180.;

}

drms_free_array(disambArray);

}

Line 1884 void computeSWIndex(struct swIndex *swKe

Line 1759 void computeSWIndex(struct swIndex *swKe

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

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

//Use magnetogram map to compute R

DRMS_Segment_t *losSeg = drms_segment_lookup(inRec, "magnetogram");

DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);

float *los = (float *) losArray->data; // los

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

Line 1914 void computeSWIndex(struct swIndex *swKe

Line 1784 void computeSWIndex(struct swIndex *swKe

// convert cdelt1_orig from degrees to arcsec

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

int nx1 = nx*cdelt1/2;

int ny1 = ny*cdelt1/2;

// Temp arrays

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

Line 1939 void computeSWIndex(struct swIndex *swKe

Line 1807 void computeSWIndex(struct swIndex *swKe

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

// 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(nx1*ny1*sizeof(float));

//spaceweather quantities computed

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

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

{

Line 2046 void computeSWIndex(struct swIndex *swKe

Line 1906 void computeSWIndex(struct swIndex *swKe

swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;

}

if (computeR(bz_err, los , dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,

p1p, p1n, p1, pmap, nx1, ny1))

{

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

}

// Clean up the arrays

drms_free_array(bitmaskArray); // Dec 18 2012 Xudong

Line 2070 void computeSWIndex(struct swIndex *swKe

Line 1923 void computeSWIndex(struct swIndex *swKe

free(bt_err); free(bh_err); free(jz_err);

free(jz_err_squared); free(jz_rms_err);

free(jz_err_squared_smooth);

free(rim);

free(p1p0);

free(p1n0);

free(p1p);

free(p1n);

free(p1);

free(pmap);

}

Line 2120 void setSWIndex(DRMS_Record_t *outRec, s

Line 1964 void setSWIndex(DRMS_Record_t *outRec, s

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

drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);

};

Line 2128 void setSWIndex(DRMS_Record_t *outRec, s

Line 1971 void setSWIndex(DRMS_Record_t *outRec, s

void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)

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

{

copy_me_keys(inRec, outRec);

copy_me_keys(bharpRec, outRec);

copy_patch_keys(inRec, outRec);

copy_patch_keys(mharpRec, outRec); // Dec 30

copy_geo_keys(inRec, outRec);

copy_geo_keys(mharpRec, outRec); // Dec 30

copy_ambig_keys(inRec, outRec);

copy_ambig_keys(bharpRec, outRec);

int status = 0;

// Change a few geometry keywords for CEA & cutout records

// Change a few geometry keywords for CEA records

if (mInfo != NULL) { // CEA

if (mInfo != NULL) {

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

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

Line 2158 void setKeys(DRMS_Record_t *outRec, DRMS

Line 2000 void setKeys(DRMS_Record_t *outRec, DRMS

drms_setkey_string(outRec, "CTYPE2", key);

drms_setkey_float(outRec, "CROTA2", 0.0);

// Jan 2 2014 XS

int nSeg = ARRLENGTH(CEASegs);

for (int iSeg = 0; iSeg < nSeg; iSeg++) {

DRMS_Segment_t *outSeg = NULL;

outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);

if (!outSeg) continue;

// Set Bunit

char bunit_xxx[20];

sprintf(bunit_xxx, "BUNIT_%03d", iSeg);

//printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);

drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);

}

} else { // Cutout

float disk_xc, disk_yc;

if (fullDisk) {

disk_xc = drms_getkey_float(bharpRec, "CRPIX1", &status);

disk_yc = drms_getkey_float(bharpRec, "CRPIX2", &status);

} else {

disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);

disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);

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

}

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

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

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

float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &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.);

Line 2190 void setKeys(DRMS_Record_t *outRec, DRMS

Line 2013 void setKeys(DRMS_Record_t *outRec, DRMS

drms_setkey_float(outRec, "CRVAL1", 0);

drms_setkey_float(outRec, "CRVAL2", 0);

// Jan 2 2014 XS

int nSeg = ARRLENGTH(CutSegs);

for (int iSeg = 0; iSeg < nSeg; iSeg++) {

DRMS_Segment_t *outSeg = NULL;

outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);

if (!outSeg) continue;

// Set Bunit

char bunit_xxx[20];

sprintf(bunit_xxx, "BUNIT_%03d", iSeg);

//printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);

drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);

}

char timebuf[1024];

}

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

double val;

TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */

tnow = (double)time(NULL);

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

tnow += UNIX_epoch;

drms_setkey_double(outRec, "DATE_B", val);

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

val = drms_getkey_time(bharpRec, "DATE", &status);

drms_setkey_string(outRec, "DATE", timebuf);

drms_setkey_time(outRec, "DATE_B", val);

drms_setkey_time(outRec, "DATE", tnow);

// set cvs commit version into keyword HEADER

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

Line 2225 void setKeys(DRMS_Record_t *outRec, DRMS

Line 2035 void setKeys(DRMS_Record_t *outRec, DRMS

};

/* ############# Nearest neighbour interpolation ############### */

Legend:

Removed from v.1.21
changed lines
	Added in v.1.22

Karen Tian