ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/JSOC/proj/sharp/apps/sharp.c

Comparing proj/sharp/apps/sharp.c (file contents):
Revision 1.4 by xudong, Fri Sep 7 22:08:36 2012 UTC vs.
Revision 1.10 by xudong, Wed Jan 23 21:48:26 2013 UTC

#	Line 19 | Line 19
19		*      Version:
20		*              v0.0    Jul 02 2012
21		*              v0.1    Jul 23 2012
22	<	*              v0.2    Sep 04 2012
22	>	*              v0.2    Sep 04 2012
23	>	*              v0.3    Dec 18 2012
24	>	*              v0.4    Jan 02 2013
25	>	*              v0.5    Jam 23 2012
26		*
27		*      Notes:
28		*              v0.0
#	Line 33 | Line 36
36		*              SW indices fixed
37		*              Added doppler and continuum
38		*              Added other keywords: HEADER (populated by cvs build version), DATE_B
39	+	*              v0.3
40	+	*              Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"
41	+	*              v0.4
42	+	*              Took out convert_inplace(). Was causing all the images to be int
43	+	*              v0.5
44	+	*              Corrected ephemeris keywords, added argument mInfo for setKeys()
45		*
46		*      Example:
47		*      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
48	<	"bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
49	<	"dop=hmi.V_720s[2012.02.20_10:00]" \
50	<	"cont=hmi.Ic_720s[2012.02.20_10:00]" \
51	<	"sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
48	>	"bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
49	>	"dop=hmi.V_720s[2012.02.20_10:00]" \
50	>	"cont=hmi.Ic_720s[2012.02.20_10:00]" \
51	>	"sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
52		*      For comparison:
53		*      bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
54	<	"out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
54	>	"out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
55		*
56		*
57		*/
58	<
58	>
59		#include <stdio.h>
60		#include <stdlib.h>
61		#include <time.h>
#	Line 69 | Line 78
78		#include <mkl_lapack.h>
79		#include <mkl_vml_functions.h>
80		#include <omp.h>
81	<
81	>
82		#define PI              (M_PI)
83		#define RADSINDEG               (PI/180.)
84		#define RAD2ARCSEC              (648000./M_PI)
#	Line 93 | Line 102
102		#define SHOW(msg) {printf("%s", msg); fflush(stdout);}
103
104		#define kNotSpecified "Not Specified"
105	<
105	>
106		// Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
107		// and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
108		// PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun from disk center.
#	Line 149 | Line 158 | enum projection {
158		lambert
159		};
160
161	+	// WSC code
162	+	char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
163	+	"SIN", "ZEA"};
164	+
165		// Ephemeris
166		struct ephemeris {
167		double disk_lonc, disk_latc;
#	Line 217 | Line 230 | void vectorTransform(float *bx_map, floa
230
231		/* Map and propogate errors */
232		int getBErr(float *bx_err, float *by_err, float *bz_err,
233	<	DRMS_Record_t *inRec, struct mapInfo *mInfo);
233	>	DRMS_Record_t *inRec, struct mapInfo *mInfo);
234
235		/* Read full disk vector magnetogram */
236		int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
#	Line 247 | Line 260 | void computeSWIndex(struct swIndex *swKe
260		void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
261
262		/* Set all keywords, no error checking for now */
263	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);
263	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo);
264
265		// ===================
266
#	Line 288 | Line 301 | char *CutSegs[] = {"magnetogram", "bitma
301		"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
302		"disambig", "conf_disambig"};
303		char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",
304	<	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
304	>	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
305		/* ========================================================================================================== */
306
307
308
309		char *module_name = "sharp";
310	<	char *version_id = "2012 Jul 02";  /* Version number */
310	>	char *version_id = "2012 Dec 18";  /* Version number */
311
312		ModuleArgs_t module_args[] =
313		{
#	Line 344 | Line 357 | int DoIt(void)
357		/* Start */
358
359		printf("==============\nStart. %d image(s) in total.\n", nrecs);
360	<
360	>
361		for (irec = 0; irec < nrecs; irec++) {
362
363		/* Records in work */
#	Line 354 | Line 367 | int DoIt(void)
367		bharpRec = bharpRS->records[irec];
368
369		TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
370	<
370	>
371		struct swIndex swKeys;
372
373		DRMS_Record_t *dopRec = NULL, *contRec = NULL;
#	Line 366 | Line 379 | int DoIt(void)
379		printf("Fetching continuum failed, image #%d skipped.\n", irec);
380		continue;
381		}
382	<
382	>
383		/* Create CEA record */
384
385		DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
#	Line 404 | Line 417 | int DoIt(void)
417		printf("Image #%d done.\n", irec);
418
419		} // irec
420	+
421
422		drms_close_records(mharpRS, DRMS_FREE_RECORD);
423		drms_close_records(bharpRS, DRMS_FREE_RECORD);
424	+	drms_close_records(dopRS, DRMS_FREE_RECORD);                            // Dec 18 2012
425	+	drms_close_records(contRS, DRMS_FREE_RECORD);                           // Dec 18 2012
426
427		return 0;
428	<
428	>
429		}       // DoIt
430
431
#	Line 435 | Line 451 | int getInputRS(DRMS_RecordSet_t **mharpR
451
452		*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
453		if (status || (*bharpRS_ptr)->n == 0) return 1;
454	<
454	>
455		if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
456
457		return 0;
#	Line 490 | Line 506 | int getInputRS_aux(DRMS_RecordSet_t **in
506
507		// Check if all T_rec are available, need to match both ways
508		int n = harpRS->n, n0 = (*inRS_ptr)->n;
509	<
509	>
510		for (int i0 = 0; i0 < n0; i0++) {
511		DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
512		TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
#	Line 565 | Line 581 | int createCeaRecord(DRMS_Record_t *mharp
581		mInfo.nbin = NBIN;
582
583		// Get ephemeris
584	<
584	>
585		if (getEphemeris(mharpRec, &(mInfo.ephem))) {
586		SHOW("CEA: get ephemeris error\n");
587		return 1;
#	Line 591 | Line 607 | int createCeaRecord(DRMS_Record_t *mharp
607		findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
608
609		// Mapping single segment: Mharp, etc.
610	<
610	>
611		if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
612		SHOW("CEA: mapping magnetogram error\n");
613		return 1;
#	Line 613 | Line 629 | int createCeaRecord(DRMS_Record_t *mharp
629		return 1;
630		}
631		printf("Intensitygram mapping done.\n");
632	<
632	>
633		if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
634		SHOW("CEA: mapping conf_disambig error\n");
635		return 1;
#	Line 646 | Line 662 | int createCeaRecord(DRMS_Record_t *mharp
662		DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
663		if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
664
665	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
665	>	setKeys(sharpRec, bharpRec, &mInfo);            // Set all other keywords
666		drms_copykey(sharpRec, mharpRec, "QUALITY");            // copied from los records
667	<
667	>
668		// Space weather
669
670		computeSWIndex(swKeys_ptr, sharpRec, &mInfo);           // compute it!
671		printf("Space weather indices done.\n");
672
673		setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
674	<
674	>
675		// Stats
676
677		int nCEASegs = ARRLENGTH(CEASegs);
#	Line 663 | Line 679 | int createCeaRecord(DRMS_Record_t *mharp
679		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
680		DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
681		int stat = set_statistics(outSeg, outArray, 1);
682	<	//              printf("%d => %d\n", iSeg, stat);
682	>	//              printf("%d => %d\n", iSeg, stat);
683		drms_free_array(outArray);
684		}
685
#	Line 722 | Line 738 | int mapScaler(DRMS_Record_t *sharpRec, D
738
739		float *map = (float *) (malloc(nxny * sizeof(float)));
740		if (performSampling(map, inData, mInfo))
741	<	{if (inArray) drms_free_array(inArray); free(map); return 1;}
741	>	{if (inArray) drms_free_array(inArray); free(map); return 1;}
742
743		// Write out
744
#	Line 730 | Line 746 | int mapScaler(DRMS_Record_t *sharpRec, D
746		outSeg = drms_segment_lookup(sharpRec, segName);
747		if (!outSeg) return 1;
748
749	<	DRMS_Type_t arrayType = outSeg->info->type;
749	>	//      DRMS_Type_t arrayType = outSeg->info->type;
750		DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
751		if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
752
753		// convert to needed data type
754
755	<	drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);
755	>	//      drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);         // Jan 02 2013
756
757		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
758	<	outArray->parent_segment = outSeg;
758	>	//      outArray->parent_segment = outSeg;
759		outArray->israw = 0;            // always compressed
760		outArray->bzero = outSeg->bzero;
761		outArray->bscale = outSeg->bscale;
762	<
762	>
763		status = drms_segment_write(outSeg, outArray, 0);
764		if (status) return 0;
765
766		if (inArray) drms_free_array(inArray);
767	+	if ((xsz != FOURK) || (ysz != FOURK)) free(inData);                     // Dec 18 2012
768		if (outArray) drms_free_array(outArray);
769		return 0;
770
#	Line 781 | Line 798 | int mapVectorB(DRMS_Record_t *sharpRec,
798		// Mapping
799
800		float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation
801	<
801	>
802		bx_map = (float *) (malloc(nxny * sizeof(float)));
803		if (performSampling(bx_map, bx_img, mInfo))
804	<	{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
805	<
804	>	{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
805	>
806		by_map = (float *) (malloc(nxny * sizeof(float)));
807		if (performSampling(by_map, by_img, mInfo))
808	<	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
809	<
808	>	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
809	>
810		bz_map = (float *) (malloc(nxny * sizeof(float)));
811		if (performSampling(bz_map, bz_img, mInfo))
812	<	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
812	>	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
813
814		free(bx_img); free(by_img); free(bz_img);
815
#	Line 814 | Line 831 | int mapVectorB(DRMS_Record_t *sharpRec,
831		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
832		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
833		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
834	<	outArray->parent_segment = outSeg;
834	>	//              outArray->parent_segment = outSeg;
835		outArray->israw = 0;
836		outArray->bzero = outSeg->bzero;
837		outArray->bscale = outSeg->bscale;
#	Line 853 | Line 870 | int mapVectorBErr(DRMS_Record_t *sharpRe
870		free(bx_err); free(by_err); free(bz_err);
871		return 1;
872		}
873	<
873	>
874		// Write out
875
876		DRMS_Segment_t *outSeg;
#	Line 866 | Line 883 | int mapVectorBErr(DRMS_Record_t *sharpRe
883		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
884		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
885		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
886	<	outArray->parent_segment = outSeg;
886	>	//              outArray->parent_segment = outSeg;
887		outArray->israw = 0;
888		outArray->bzero = outSeg->bzero;
889		outArray->bscale = outSeg->bscale;
#	Line 965 | Line 982 | int getEphemeris(DRMS_Record_t *inRec, s
982		float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
983		float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
984		float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
985	<	printf("cdelt=%f\n",cdelt);
985	>	printf("cdelt=%f\n",cdelt);
986		ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0
987		ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
988
#	Line 1081 | Line 1098 | int performSampling(float *outData, floa
1098
1099		float *xi_out = mInfo->xi_out;
1100		float *zeta_out = mInfo->zeta_out;
1101	<
1101	>
1102		// Interpolation
1103
1104		struct fint_struct pars;
#	Line 1107 | Line 1124 | int performSampling(float *outData, floa
1124		// Rebinning, smoothing
1125
1126		frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
1127	+	free(outData0);         // Dec 18 2012
1128
1129		//
1130
#	Line 1176 | Line 1194 | void vectorTransform(float *bx_map, floa
1194
1195		}
1196		}
1197	<
1197	>
1198		}
1199
1200
#	Line 1187 | Line 1205 | void vectorTransform(float *bx_map, floa
1205		*/
1206
1207		int getBErr(float *bx_err, float *by_err, float *bz_err,
1208	<	DRMS_Record_t *inRec, struct mapInfo *mInfo)
1208	>	DRMS_Record_t *inRec, struct mapInfo *mInfo)
1209		{
1210
1211		int status = 0;
#	Line 1245 | Line 1263 | int getBErr(float *bx_err, float *by_err
1263		int ind_map, ind_img;
1264
1265		double bpSigma2, btSigma2, brSigma2;            // variances after prop
1266	<
1266	>
1267		for (int row = 0; row < mInfo->nrow; row++) {
1268		for (int col = 0; col < mInfo->ncol; col++) {
1269
#	Line 1315 | Line 1333 | int readVectorB(DRMS_Record_t *inRec, fl
1333
1334		DRMS_Segment_t *inSeg;
1335		DRMS_Array_t *inArray_ambig;
1336	<	DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1336	>	DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1337
1338		char *ambig;
1339		float *bTotal, *bAzim, *bIncl;
#	Line 1410 | Line 1428 | int readVectorBErr(DRMS_Record_t *inRec,
1428
1429		float *data_ptr[9];
1430		char *segName[9] = {"field", "inclination", "azimuth",
1431	<	"field_err", "inclination_err", "azimuth_err",
1432	<	"field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1431	>	"field_err", "inclination_err", "azimuth_err",
1432	>	"field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1433		DRMS_Segment_t *inSegs[9];
1434		DRMS_Array_t *inArrays[9];
1435
#	Line 1447 | Line 1465 | int readVectorBErr(DRMS_Record_t *inRec,
1465		errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG;
1466		errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG;
1467		errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1468	<
1468	>
1469		}
1470
1471		//
1472
1473		for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
1474	<
1474	>
1475		return 0;
1476
1477		}
#	Line 1528 | Line 1546 | int createCutRecord(DRMS_Record_t *mharp
1546		if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
1547
1548		setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
1549	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
1550	<
1549	>	setKeys(sharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout
1550	>
1551		// Stats
1552	<
1552	>
1553		int nCutSegs = ARRLENGTH(CutSegs);
1554		for (int iSeg = 0; iSeg < nCutSegs; iSeg++) {
1555		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
#	Line 1539 | Line 1557 | int createCutRecord(DRMS_Record_t *mharp
1557		set_statistics(outSeg, outArray, 1);
1558		drms_free_array(outArray);
1559		}
1560	<
1560	>
1561		return 0;
1562
1563		}
#	Line 1584 | Line 1602 | int writeCutout(DRMS_Record_t *outRec, D
1602		} else {
1603		return 1;
1604		}
1605	<
1605	>
1606		/* Adding disambiguation resolution to cutout azimuth? */
1607	<
1607	>
1608		#if DISAMB_AZI
1609		if (!strcmp(SegName, "azimuth")) {
1610	<	DRMS_Segment_t *disambSeg = drms_segment_lookup(inRec, "disambig");
1610	>	DRMS_Segment_t *disambSeg = NULL;
1611	>	disambSeg = drms_segment_lookup(inRec, "disambig");
1612		if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
1613		DRMS_Array_t *disambArray;
1614		if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1615		disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1616		if (status) {drms_free_array(cutoutArray); return 1;}
1617		} else {
1618	+	drms_free_array(cutoutArray);
1619		return 1;
1620		}
1621		double *azimuth = (double *) cutoutArray->data;
#	Line 1606 | Line 1626 | int writeCutout(DRMS_Record_t *outRec, D
1626		drms_free_array(disambArray);
1627		}
1628		#endif
1629	<
1629	>
1630		/* Write out */
1631
1632		outSeg = drms_segment_lookup(outRec, SegName);
1633		if (!outSeg) return 1;
1634	<	drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);
1634	>	//      drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);      // Jan 02 2013
1635		outSeg->axis[0] = cutoutArray->axis[0];
1636		outSeg->axis[1] = cutoutArray->axis[1];
1637	<	cutoutArray->parent_segment = outSeg;
1637	>	//      cutoutArray->parent_segment = outSeg;
1638		cutoutArray->israw = 0;         // always compressed
1639		cutoutArray->bzero = outSeg->bzero;
1640		cutoutArray->bscale = outSeg->bscale;               // Same as inArray's
#	Line 1643 | Line 1663 | void computeSWIndex(struct swIndex *swKe
1663		int dims[2] = {nx, ny};
1664		// Get bx, by, bz, mask
1665
1666	<	// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1667	<	//DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "bitmap");
1668	<	//DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1669	<	//int *mask = (int *) maskArray->data;          // get the previously made mask array
1670	<
1671	<	//Use conf_disambig map as a threshold on spaceweather quantities
1666	>	// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1667	>	DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
1668	>	DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
1669	>	int *bitmask = (int *) bitmaskArray->data;              // get the previously made mask array
1670	>
1671	>	//Use conf_disambig map as a threshold on spaceweather quantities
1672		DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1673		DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1674		int *mask = (int *) maskArray->data;            // get the previously made mask array
1675	<
1675	>
1676		DRMS_Segment_t *bxSeg = drms_segment_lookup(inRec, BP_SEG_CEA);
1677		DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
1678		float *bx = (float *) bxArray->data;            // bx
#	Line 1669 | Line 1689 | void computeSWIndex(struct swIndex *swKe
1689		// Get emphemeris
1690
1691		//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1692	<	float cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);
1693	<	float dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1692	>	float  cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);
1693	>	float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1694		double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);
1695		double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);
1696		float imcrpix1    = drms_getkey_float(inRec, "IMCRPIX1", &status);
#	Line 1679 | Line 1699 | void computeSWIndex(struct swIndex *swKe
1699		float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);
1700
1701		//float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)
1702	<
1703	<	printf("cdelt1=%f\n",cdelt1);
1704	<	printf("rsun_ref=%f\n",rsun_ref);
1705	<	printf("rsun_obs=%f\n",rsun_obs);
1706	<	printf("dsun_obs=%f\n",dsun_obs);
1707	<
1702	>
1703	>	printf("cdelt1=%f\n",cdelt1);
1704	>	printf("rsun_ref=%f\n",rsun_ref);
1705	>	printf("rsun_obs=%f\n",rsun_obs);
1706	>	printf("dsun_obs=%f\n",dsun_obs);
1707	>
1708		// Temp arrays
1709
1710		float *bh = (float *) (malloc(nxny * sizeof(float)));
1711		float *bt = (float *) (malloc(nxny * sizeof(float)));
1712		float *jz = (float *) (malloc(nxny * sizeof(float)));
1713	+	float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
1714		float *bpx = (float *) (malloc(nxny * sizeof(float)));
1715		float *bpy = (float *) (malloc(nxny * sizeof(float)));
1716		float *bpz = (float *) (malloc(nxny * sizeof(float)));
#	Line 1701 | Line 1722 | void computeSWIndex(struct swIndex *swKe
1722		float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
1723		float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1724		float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1725	<
1726	<	// Compute
1725	>
1726	>	// The Computations
1727
1728		if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1729	<	mask, cdelt1, rsun_ref, rsun_obs)){
1729	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs)){
1730		swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
1731		swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1732		}
1733
1734		for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
1735	<	greenpot(bpx, bpy, bpz, nx, ny);
1735	>	greenpot(bpx, bpy, bpz, nx, ny);
1736
1737	<	computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask);
1737	>	computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
1738
1739	<	if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask))
1739	>	if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask, bitmask))
1740		swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
1741
1742	<	computeB_total(bx, by, bz, bt, dims, mask);
1742	>	computeB_total(bx, by, bz, bt, dims, mask, bitmask);
1743
1744	<	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt))
1744	>	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt))
1745		swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
1746
1747	<	if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh))
1747	>	if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, bitmask, derx_bh, dery_bh))
1748		swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
1749
1750	<	if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, derx_bz, dery_bz))
1750	>	if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, bitmask, derx_bz, dery_bz))
1751		swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1752
1753	+	computeJz(bx, by, dims, jz, mask, bitmask, cdelt1, rsun_ref, rsun_obs, derx, dery);
1754
1755	<
1756	<	if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask,
1757	<	cdelt1, rsun_ref, rsun_obs, derx, dery)) {
1758	<	swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
1755	>	struct fresize_struct convolution_array;
1756	>	init_fresize_gaussian(&convolution_array, 4.0f, 12, 1);
1757	>	fresize(&convolution_array, jz, jz_smooth, nx, ny, nx, nx, ny, nx, 0, 0, 0.0f);
1758	>
1759	>	if(computeJzsmooth(bx, by, dims, jz_smooth, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask,
1760	>	cdelt1, rsun_ref, rsun_obs, derx, dery))
1761	>	{
1762	>	swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
1763		swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
1764		}
1765	<
1740	<	printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
1765	>	printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
1766
1767	<	if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, cdelt1, rsun_ref, rsun_obs))
1767	>
1768	>	if (computeAlpha(bz, dims, jz_smooth, &(swKeys_ptr->mean_alpha), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1769		swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
1770
1771	<	if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih),
1771	>	if (computeHelicity(bz, dims, jz_smooth, &(swKeys_ptr->mean_ih),
1772		&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
1773	<	mask, cdelt1, rsun_ref, rsun_obs)) {
1773	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
1774		swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
1775		swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
1776		swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
1777		}
1778
1779	<	if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),
1780	<	mask, cdelt1, rsun_ref, rsun_obs))
1779	>	if (computeSumAbsPerPolarity(bz, jz_smooth, dims, &(swKeys_ptr->totaljz),
1780	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1781		swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
1782	<
1782	>
1783
1784		if (computeFreeEnergy(bx, by, bpx, bpy, dims,
1785		&(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
1786	<	mask, cdelt1, rsun_ref, rsun_obs)) {
1786	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
1787		swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1788		swKeys_ptr->totpot = DRMS_MISSING_FLOAT;
1789		}
1790	<
1790	>
1791		if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,
1792		&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),
1793		&(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
1794	<	mask)) {
1794	>	mask, bitmask)) {
1795		swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1796		swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
1797		swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1798		swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;
1799		}
1800	<
1800	>
1801		// Clean up
1802
1803	+	drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong
1804		drms_free_array(maskArray);
1805	<	drms_free_array(bxArray);
1805	>	drms_free_array(bxArray);
1806		drms_free_array(byArray);
1807		drms_free_array(bzArray);
1808
1809	<	free(bh); free(bt); free(jz);
1809	>	free(bh); free(bt); free(jz); free(jz_smooth);
1810		free(bpx); free(bpy); free(bpz);
1811		free(derx); free(dery);
1812		free(derx_bt); free(dery_bt);
1813		free(derx_bz); free(dery_bz);
1814		free(derx_bh); free(dery_bh);
1815
1816	+	free_fresize(&convolution_array);
1817		}
1818
1819
#	Line 1820 | Line 1848 | void setSWIndex(DRMS_Record_t *outRec, s
1848		*
1849		*/
1850
1851	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)
1851	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo)
1852		{
1853	<	copy_me_keys(inRec, outRec);
1854	<	copy_patch_keys(inRec, outRec);
1855	<	copy_geo_keys(inRec, outRec);
1856	<	copy_ambig_keys(inRec, outRec);
1857	<
1858	<	char timebuf[1024];
1859	<	float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1860	<	double val;
1861	<	int status = DRMS_SUCCESS;
1862	<
1863	<	val = drms_getkey_double(inRec, "DATE",&status);
1864	<	drms_setkey_double(outRec, "DATE_B", val);
1865	<	sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
1866	<	drms_setkey_string(outRec, "DATE", timebuf);
1867	<
1868	<	// set cvs commit version into keyword HEADER
1869	<	char *cvsinfo = strdup("$Header$");
1870	<	//   status = drms_setkey_string(outRec, "HEADER", cvsinfo);
1853	>	copy_me_keys(inRec, outRec);
1854	>	copy_patch_keys(inRec, outRec);
1855	>	copy_geo_keys(inRec, outRec);
1856	>	copy_ambig_keys(inRec, outRec);
1857	>
1858	>	int status = 0;
1859	>
1860	>	// Change a few geometry keywords for CEA records
1861	>	if (mInfo != NULL) {
1862	>
1863	>	drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
1864	>	drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
1865	>
1866	>	drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
1867	>	drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
1868	>	drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
1869	>	drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
1870	>	drms_setkey_string(outRec, "CUNIT1", "degree");
1871	>	drms_setkey_string(outRec, "CUNIT2", "degree");
1872	>
1873	>	char key[64];
1874	>	snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
1875	>	drms_setkey_string(outRec, "CTYPE1", key);
1876	>	snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
1877	>	drms_setkey_string(outRec, "CTYPE2", key);
1878	>	drms_setkey_float(outRec, "CROTA2", 0.0);
1879	>
1880	>	} else {
1881	>
1882	>	float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status);
1883	>	float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status);
1884	>	float x_ll = drms_getkey_float(inRec, "CRPIX1", &status);
1885	>	float y_ll = drms_getkey_float(inRec, "CRPIX2", &status);
1886	>	// Defined as disk center's pixel address wrt lower-left of cutout
1887	>	drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
1888	>	drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
1889	>	// Always 0.
1890	>	drms_setkey_float(outRec, "CRVAL1", 0);
1891	>	drms_setkey_float(outRec, "CRVAL2", 0);
1892	>
1893	>	}
1894	>
1895	>	char timebuf[1024];
1896	>	float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1897	>	double val;
1898	>	status = 0;
1899	>
1900	>	val = drms_getkey_double(inRec, "DATE",&status);
1901	>	drms_setkey_double(outRec, "DATE_B", val);
1902	>	sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
1903	>	drms_setkey_string(outRec, "DATE", timebuf);
1904	>
1905	>	// set cvs commit version into keyword HEADER
1906	>	char *cvsinfo = strdup("$Header$");
1907	>	//   status = drms_setkey_string(outRec, "HEADER", cvsinfo);
1908		status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);
1909	<
1909	>
1910		};
1911
1912		//

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines