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.2 by xudong, Mon Aug 27 19:56:00 2012 UTC vs.
Revision 1.18 by mbobra, Sat Nov 2 20:31:08 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
23	+	*              v0.3    Dec 18 2012
24	+	*              v0.4    Jan 02 2013
25	+	*    v0.5  Jan 23 2013
26	+	*              v0.6    Aug 12 2013
27		*
28		*      Notes:
29		*              v0.0
#	Line 31 | Line 36
36		*              Fixed char I/O thanks to Art
37		*              SW indices fixed
38		*              Added doppler and continuum
39	+	*              Added other keywords: HEADER (populated by cvs build version), DATE_B
40	+	*              v0.3
41	+	*              Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"
42	+	*              v0.4
43	+	*              Took out convert_inplace(). Was causing all the images to be int
44	+	*    v0.5
45	+	*    Corrected ephemeris keywords, added argument mInfo for setKeys()
46	+	*              v0.6
47	+	*              Changes in remapping of bitmap and conf_disambig, now near neighbor without anti-aliasing
48	+	*
49		*
50		*      Example:
51		*      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
52	<	"bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
53	<	"dop=hmi.V_720s[2012.02.20_10:00]" \
54	<	"cont=hmi.Ic_720s[2012.02.20_10:00]" \
55	<	"sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
52	>	"bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
53	>	"dop=hmi.V_720s[2012.02.20_10:00]" \
54	>	"cont=hmi.Ic_720s[2012.02.20_10:00]" \
55	>	"sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
56		*      For comparison:
57		*      bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
58	<	"out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
58	>	"out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
59		*
60		*
61		*/
62	<
62	>
63		#include <stdio.h>
64		#include <stdlib.h>
65		#include <time.h>
#	Line 62 | Line 77
77		#include "errorprop.c"
78		#include "sw_functions.c"
79
80	+	//#include <mkl.h> // Comment out mkl.h, which can only run on solar3
81		#include <mkl_blas.h>
82		#include <mkl_service.h>
83		#include <mkl_lapack.h>
#	Line 76 | Line 92
92		#define FOURK2    (16777216)
93
94		#define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
95	<
95	>
96		// Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
97		#define NYQVIST         (0.015)
98	<
98	>
99	>	// Some other things
100		#ifndef MIN
101		#define MIN(a,b) (((a)<(b)) ? (a) : (b))
102		#endif
#	Line 91 | Line 108
108		#define SHOW(msg) {printf("%s", msg); fflush(stdout);}
109
110		#define kNotSpecified "Not Specified"
111	<
111	>
112		// Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
113		// and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
114		// 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 107 | Line 124
124		#define INTERP                  0
125		#define dpath    "/home/jsoc/cvs/Development/JSOC"
126
127	+
128		/* ========================================================================================================== */
129
130		// Space weather keywords
131		struct swIndex {
132	<	float mean_vf;
132	>	float mean_vf;
133	>	float count_mask;
134		float absFlux;
135		float mean_hf;
136		float mean_gamma;
137		float mean_derivative_btotal;
138	<	float mean_derivative_bh;
138	>	float mean_derivative_bh;
139		float mean_derivative_bz;
140		float mean_jz;
141		float us_i;
#	Line 130 | Line 149 | struct swIndex {
149		float area_w_shear_gt_45;
150		float meanshear_angle;
151		float area_w_shear_gt_45h;
152	<	float meanshear_angleh;
152	>	float meanshear_angleh;
153	>	float mean_derivative_btotal_err;
154	>	float mean_vf_err;
155	>	float mean_gamma_err;
156	>	float mean_derivative_bh_err;
157	>	float mean_derivative_bz_err;
158	>	float mean_jz_err;
159	>	float us_i_err;
160	>	float mean_alpha_err;
161	>	float mean_ih_err;
162	>	float total_us_ih_err;
163	>	float total_abs_ih_err;
164	>	float totaljz_err;
165	>	float meanpot_err;
166	>	float totpot_err;
167	>	float meanshear_angle_err;
168		};
169
170		// Mapping method
#	Line 147 | Line 181 | enum projection {
181		lambert
182		};
183
184	<	// Ephemeris
184	>	// WSC code
185	>	char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
186	>	"SIN", "ZEA"};
187	>
188	>	// Ephemeris information
189		struct ephemeris {
190		double disk_lonc, disk_latc;
191		double disk_xc, disk_yc;
#	Line 165 | Line 203 | struct mapInfo {
203		float *xi_out, *zeta_out;       // coordinate on full disk image to sample at
204		};
205
168	–
206		/* ========================================================================================================== */
207
208		/* Get all input data series */
#	Line 181 | Line 218 | int getInputRS_aux(DRMS_RecordSet_t **in
218		/* Find record from record set with given T_rec */
219		int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
220
184	–	// ===================
185	–
221		/* Create CEA record */
222		int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
223	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
223	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
224		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
225
226		/* Mapping single segment, wrapper */
#	Line 208 | Line 243 | int getEphemeris(DRMS_Record_t *inRec, s
243		void findCoord(struct mapInfo *mInfo);
244
245		/* Mapping function */
246	<	int performSampling(float *outData, float *inData, struct mapInfo *mInfo);
246	>	int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
247
248		/* Performing local vector transformation */
249		void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
250
251		/* Map and propogate errors */
252		int getBErr(float *bx_err, float *by_err, float *bz_err,
253	<	DRMS_Record_t *inRec, struct mapInfo *mInfo);
253	>	DRMS_Record_t *inRec, struct mapInfo *mInfo);
254
255		/* Read full disk vector magnetogram */
256		int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
257
258		/* Read variances and covariances of vector magnetograms */
259	<	int readVectorBErr(DRMS_Record_t *bharpRec,
259	>	int readVectorBErr(DRMS_Record_t *bharpRec,
260		float *bT, float *bI, float *bA,
261	<	float *errbT, float *errbI, float *errbA,
261	>	float *errbT, float *errbI, float *errbA,
262		float *errbTbI, float *errbTbA, float *errbIbA);
263
264		// ===================
#	Line 245 | Line 280 | void computeSWIndex(struct swIndex *swKe
280		void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
281
282		/* Set all keywords, no error checking for now */
283	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);
283	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo);
284
285		// ===================
286
#	Line 276 | Line 311 | char *BharpSegs[] = {"inclination", "azi
311		"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
312		"disambig", "conf_disambig"};
313		// For stats
314	<	char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
314	>	char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
315		"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
316		"damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
317		"conv_flag", // fixed
#	Line 285 | Line 320 | char *CutSegs[] = {"magnetogram", "bitma
320		"field_inclination_err", "field_az_err", "inclin_azimuth_err",
321		"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
322		"disambig", "conf_disambig"};
323	<	char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",
324	<	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
290	<	/* ========================================================================================================== */
291	<
323	>	char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",
324	>	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
325
326	+	/* ========================================================================================================== */
327
328		char *module_name = "sharp";
329	<	char *version_id = "2012 Jul 02";  /* Version number */
329	>	int seed;
330
331		ModuleArgs_t module_args[] =
332		{
#	Line 302 | Line 336 | ModuleArgs_t module_args[] =
336		{ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
337		{ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
338		{ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
339	+	{ARG_INT,    "seed", "987654", "Seed for the random number generator."},
340		{ARG_END}
341		};
342
343	<	int DoIt(void)
343	>	int DoIt(void)
344		{
345	<
345	>	int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
346	>	int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
347	>
348		int status = DRMS_SUCCESS;
349		int nrecs, irec;
350
#	Line 326 | Line 363 | int DoIt(void)
363		contQuery = (char *) params_get_str(&cmdparams, "cont");
364		sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
365		sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
366	+
367	+	seed = params_get_int(&cmdparams, "seed");
368
369		/* Get input data, check everything */
370
#	Line 333 | Line 372 | int DoIt(void)
372		DIE("Input harp data error.");
373		nrecs = mharpRS->n;
374
375	<	if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
375	>	if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
376		DIE("Input doppler data error.");
377
378	<	if (getInputRS_aux(&contRS, contQuery, mharpRS))
378	>	if (getInputRS_aux(&contRS, contQuery, mharpRS))
379		DIE("Input continuum data error.");
380
381		/* Start */
382
383		printf("==============\nStart. %d image(s) in total.\n", nrecs);
384	<
384	>
385		for (irec = 0; irec < nrecs; irec++) {
386
387		/* Records in work */
#	Line 352 | Line 391 | int DoIt(void)
391		bharpRec = bharpRS->records[irec];
392
393		TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
394	<
394	>
395		struct swIndex swKeys;
396
397		DRMS_Record_t *dopRec = NULL, *contRec = NULL;
#	Line 364 | Line 403 | int DoIt(void)
403		printf("Fetching continuum failed, image #%d skipped.\n", irec);
404		continue;
405		}
406	<
406	>
407		/* Create CEA record */
408
409		DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
#	Line 402 | Line 441 | int DoIt(void)
441		printf("Image #%d done.\n", irec);
442
443		} // irec
444	+
445
446		drms_close_records(mharpRS, DRMS_FREE_RECORD);
447		drms_close_records(bharpRS, DRMS_FREE_RECORD);
448	+	drms_close_records(dopRS, DRMS_FREE_RECORD);                            // Dec 18 2012
449	+	drms_close_records(contRS, DRMS_FREE_RECORD);                           // Dec 18 2012
450
451		return 0;
452	<
452	>
453		}       // DoIt
454
455
#	Line 433 | Line 475 | int getInputRS(DRMS_RecordSet_t **mharpR
475
476		*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
477		if (status || (*bharpRS_ptr)->n == 0) return 1;
478	<
478	>
479		if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
480
481		return 0;
#	Line 460 | Line 502 | int compareHarp(DRMS_RecordSet_t *mharpR
502		for (int i = 0; i < nrecs; i++) {
503		mharpRec_t = mharpRS->records[i];
504		bharpRec_t = bharpRS->records[i];
505	<	if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
505	>	if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
506		drms_getkey_int(bharpRec_t, "HARPNUM", &status)) ||
507	<	(drms_getkey_time(mharpRec_t, "T_REC", &status) !=
507	>	(drms_getkey_time(mharpRec_t, "T_REC", &status) !=
508		drms_getkey_time(bharpRec_t, "T_REC", &status)))
509		{
510		return 1;
#	Line 473 | Line 515 | int compareHarp(DRMS_RecordSet_t *mharpR
515
516		}
517
518	<	/*
518	>	/*
519		* Get other data series, check all T_REC are available
520		*
521		*/
#	Line 488 | Line 530 | int getInputRS_aux(DRMS_RecordSet_t **in
530
531		// Check if all T_rec are available, need to match both ways
532		int n = harpRS->n, n0 = (*inRS_ptr)->n;
533	<
533	>
534		for (int i0 = 0; i0 < n0; i0++) {
535		DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
536		TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
#	Line 517 | Line 559 | int getInputRS_aux(DRMS_RecordSet_t **in
559
560		}
561
562	<	/*
562	>	/*
563		* Find record from record set with given T_rec
564		*
565		*/
#	Line 547 | Line 589 | int getInputRec_aux(DRMS_Record_t **inRe
589		*
590		*/
591
592	<	int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
593	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
592	>	int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
593	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
594		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
595		{
596
#	Line 560 | Line 602 | int createCeaRecord(DRMS_Record_t *mharp
602		mInfo.proj = (enum projection) cyleqa;          // projection method
603		mInfo.xscale = XSCALE;
604		mInfo.yscale = YSCALE;
605	<	mInfo.nbin = NBIN;
605	>
606	>	int ncol0, nrow0;             // oversampled map size
607
608		// Get ephemeris
609	<
610	<	if (getEphemeris(mharpRec, &(mInfo.ephem))) return 1;
609	>
610	>	if (getEphemeris(mharpRec, &(mInfo.ephem))) {
611	>	SHOW("CEA: get ephemeris error\n");
612	>	return 1;
613	>	}
614
615		// Find position
616
617	<	if (findPosition(mharpRec, &mInfo)) return 1;
617	>	if (findPosition(mharpRec, &mInfo)) {
618	>	SHOW("CEA: find position error\n");
619	>	return 1;
620	>	}
621	>
622	>	// ========================================
623	>	// Do this for all bitmaps, Aug 12 2013 XS
624	>	// ========================================
625	>
626	>	mInfo.nbin = 1;                       // for bitmaps. suppress anti-aliasing
627	>	ncol0 = mInfo.ncol;
628	>	nrow0 = mInfo.nrow;
629	>
630	>	mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
631	>	mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
632	>
633	>	findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
634	>
635	>	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
636	>	SHOW("CEA: mapping bitmap error\n");
637	>	return 1;
638	>	}
639	>	printf("Bitmap mapping done.\n");
640	>
641	>	if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
642	>	SHOW("CEA: mapping conf_disambig error\n");
643	>	return 1;
644	>	}
645	>	printf("Conf disambig mapping done.\n");
646	>
647	>	free(mInfo.xi_out);
648	>	free(mInfo.zeta_out);
649
650	+	// ========================================
651	+	// Do this again for floats, Aug 12 2013 XS
652	+	// ========================================
653		// Create xi_out, zeta_out array in mInfo:
654		// Coordinates to sample in original full disk image
655
656	<	int ncol0, nrow0;               // oversampled map size
656	>	mInfo.nbin = NBIN;
657		ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
658		nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
659
#	Line 583 | Line 663 | int createCeaRecord(DRMS_Record_t *mharp
663		findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
664
665		// Mapping single segment: Mharp, etc.
666	<
667	<	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) return 1;
668	<	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) return 1;
666	>
667	>	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
668	>	SHOW("CEA: mapping magnetogram error\n");
669	>	return 1;
670	>	}
671		printf("Magnetogram mapping done.\n");
672	<
673	<	if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) return 1;
672	>
673	>	if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
674	>	SHOW("CEA: mapping dopplergram error\n");
675	>	return 1;
676	>	}
677		printf("Dopplergram mapping done.\n");
678
679	<	if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) return 1;
679	>	if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
680	>	SHOW("CEA: mapping continuum error\n");
681	>	return 1;
682	>	}
683		printf("Intensitygram mapping done.\n");
684	<
597	<	if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) return 1;
598	<	printf("Conf disambig mapping done.\n");
599	<
684	>
685		// Mapping vector B
686
687	<	if (mapVectorB(sharpRec, bharpRec, &mInfo)) return 1;
687	>	if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
688	>	SHOW("CEA: mapping vector B error\n");
689	>	return 1;
690	>	}
691		printf("Vector B mapping done.\n");
692
693		// Mapping vector B errors
694
695	<	if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) return 1;
695	>	if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) {
696	>	SHOW("CEA: mapping vector B uncertainty error\n");
697	>	return 1;
698	>	}
699		printf("Vector B error done.\n");
700
701		// Keywords & Links
#	Line 617 | Line 708 | int createCeaRecord(DRMS_Record_t *mharp
708		DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
709		if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
710
711	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
712	<
711	>	setKeys(sharpRec, bharpRec, &mInfo);            // Set all other keywords
712	>	drms_copykey(sharpRec, mharpRec, "QUALITY");            // copied from los records
713	>
714		// Space weather
715
716		computeSWIndex(swKeys_ptr, sharpRec, &mInfo);           // compute it!
717		printf("Space weather indices done.\n");
718
719		setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
720	<
720	>
721		// Stats
722
723		int nCEASegs = ARRLENGTH(CEASegs);
#	Line 633 | Line 725 | int createCeaRecord(DRMS_Record_t *mharp
725		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
726		DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
727		int stat = set_statistics(outSeg, outArray, 1);
728	<	//              printf("%d => %d\n", iSeg, stat);
728	>	//              printf("%d => %d\n", iSeg, stat);
729		drms_free_array(outArray);
730		}
731
#	Line 644 | Line 736 | int createCeaRecord(DRMS_Record_t *mharp
736		}
737
738
739	<	/*
739	>	/*
740		* Mapping a single segment
741		* Read in full disk image, utilize mapImage for mapping
742		* then write the segment out, segName same in in/out Rec
#	Line 658 | Line 750 | int mapScaler(DRMS_Record_t *sharpRec, D
750		int status = 0;
751		int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
752		int dims[2] = {nx, ny};
753	+	int interpOpt = INTERP;         // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
754
755		// Input full disk array
756
#	Line 672 | Line 765 | int mapScaler(DRMS_Record_t *sharpRec, D
765		float *inData;
766		int xsz = inArray->axis[0], ysz = inArray->axis[1];
767		if ((xsz != FOURK) || (ysz != FOURK)) {         // for bitmap, make tmp full disk
768	+	interpOpt = 3;          // Aug 12 XS, near neighbor
769		float *inData0 = (float *) inArray->data;
770		inData = (float *) (calloc(FOURK2, sizeof(float)));
771		int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
#	Line 691 | Line 785 | int mapScaler(DRMS_Record_t *sharpRec, D
785		// Mapping
786
787		float *map = (float *) (malloc(nxny * sizeof(float)));
788	<	if (performSampling(map, inData, mInfo))
789	<	{if (inArray) drms_free_array(inArray); free(map); return 1;}
788	>	if (performSampling(map, inData, mInfo, interpOpt))             // Add interpOpt for different types, Aug 12 XS
789	>	{if (inArray) drms_free_array(inArray); free(map); return 1;}
790
791		// Write out
792
793	<	DRMS_Segment_t *outSeg = NULL;
793	>	DRMS_Segment_t *outSeg = NULL;
794		outSeg = drms_segment_lookup(sharpRec, segName);
795		if (!outSeg) return 1;
796
797	<	DRMS_Type_t arrayType = outSeg->info->type;
797	>	//  DRMS_Type_t arrayType = outSeg->info->type;
798		DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
799		if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
800
801		// convert to needed data type
802
803	<	drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);
803	>	//  drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);         // Jan 02 2013
804
805		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
806	<	outArray->parent_segment = outSeg;
806	>	//  outArray->parent_segment = outSeg;
807		outArray->israw = 0;            // always compressed
808		outArray->bzero = outSeg->bzero;
809		outArray->bscale = outSeg->bscale;
810	<
810	>
811		status = drms_segment_write(outSeg, outArray, 0);
812		if (status) return 0;
813
814		if (inArray) drms_free_array(inArray);
815	+	if ((xsz != FOURK) || (ysz != FOURK)) free(inData);                     // Dec 18 2012
816		if (outArray) drms_free_array(outArray);
817		return 0;
818
819		}
820
821
822	<	/*
822	>	/*
823		* Mapping vector magnetogram
824		*
825		*/
#	Line 751 | Line 846 | int mapVectorB(DRMS_Record_t *sharpRec,
846		// Mapping
847
848		float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation
849	<
849	>
850		bx_map = (float *) (malloc(nxny * sizeof(float)));
851	<	if (performSampling(bx_map, bx_img, mInfo))
852	<	{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
853	<
851	>	if (performSampling(bx_map, bx_img, mInfo, INTERP))
852	>	{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
853	>
854		by_map = (float *) (malloc(nxny * sizeof(float)));
855	<	if (performSampling(by_map, by_img, mInfo))
856	<	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
857	<
855	>	if (performSampling(by_map, by_img, mInfo, INTERP))
856	>	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
857	>
858		bz_map = (float *) (malloc(nxny * sizeof(float)));
859	<	if (performSampling(bz_map, bz_img, mInfo))
860	<	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
859	>	if (performSampling(bz_map, bz_img, mInfo, INTERP))
860	>	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
861
862		free(bx_img); free(by_img); free(bz_img);
863
#	Line 784 | Line 879 | int mapVectorB(DRMS_Record_t *sharpRec,
879		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
880		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
881		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
882	<	outArray->parent_segment = outSeg;
882	>	//              outArray->parent_segment = outSeg;
883		outArray->israw = 0;
884		outArray->bzero = outSeg->bzero;
885		outArray->bscale = outSeg->bscale;
#	Line 800 | Line 895 | int mapVectorB(DRMS_Record_t *sharpRec,
895		}
896
897
898	<	/*
898	>	/*
899		* Mapping vector magnetogram errors
900		*
901		*/
#	Line 823 | Line 918 | int mapVectorBErr(DRMS_Record_t *sharpRe
918		free(bx_err); free(by_err); free(bz_err);
919		return 1;
920		}
921	<
921	>
922		// Write out
923
924		DRMS_Segment_t *outSeg;
#	Line 836 | Line 931 | int mapVectorBErr(DRMS_Record_t *sharpRe
931		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
932		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
933		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
934	<	outArray->parent_segment = outSeg;
934	>	//              outArray->parent_segment = outSeg;
935		outArray->israw = 0;
936		outArray->bzero = outSeg->bzero;
937		outArray->bscale = outSeg->bscale;
#	Line 853 | Line 948 | int mapVectorBErr(DRMS_Record_t *sharpRe
948
949
950
951	<	/*
951	>	/*
952		* Determine reference point coordinate and patch size according to keywords
953		* xc, yc are the coordinate of patch center, in degrees
954		* ncol and nrow are the final size
#	Line 881 | Line 976 | int findPosition(DRMS_Record_t *inRec, s
976		// We compute minlon & minlat then by
977		// LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
978
979	<	float psize = drms_getkey_float(inRec, "SIZE", &status);
980	<	if (psize != psize) {
981	<	TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1;
979	>	//  float psize = drms_getkey_float(inRec, "SIZE", &status);
980	>	//  if (psize != psize) {
981	>
982	>	if (minlon != minlon || maxlon != maxlon) {         // check lons instead of SIZE
983	>	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
984		double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
985		char firstRecQuery[100], t0_str[100];
986		sprint_time(t0_str, t0, "TAI", 0);
#	Line 923 | Line 1020 | int getEphemeris(DRMS_Record_t *inRec, s
1020		int status = 0;
1021
1022		float crota2 = drms_getkey_float(inRec, "CROTA2", &status);     // rotation
1023	<	double sina = sin(crota2 * RADSINDEG);
1023	>	double sina = sin(crota2 * RADSINDEG);
1024		double cosa = cos(crota2 * RADSINDEG);
1025
1026		ephem->pa = - crota2 * RADSINDEG;
#	Line 935 | Line 1032 | int getEphemeris(DRMS_Record_t *inRec, s
1032		float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
1033		float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
1034		float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
938	–	printf("cdelt=%f\n",cdelt);
1035		ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0
1036		ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
1037
#	Line 1000 | Line 1096 | void findCoord(struct mapInfo *mInfo)
1096		x = (col0 + 0.5 - ncol0/2.) * xscale0;          // in rad
1097		y = (row0 + 0.5 - nrow0/2.) * yscale0;
1098
1099	<	/* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1099	>	/* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1100		* the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
1101		* projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
1102	<	* is the heliographic longitude and latitude of the map center. Both are in degree.
1102	>	* is the heliographic longitude and latitude of the map center. Both are in degree.
1103		*/
1104
1105		if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
#	Line 1016 | Line 1112 | void findCoord(struct mapInfo *mInfo)
1112		* image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
1113		*/
1114
1115	<	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1115	>	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1116		disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1117		xi_out[ind_map] = -1;
1118		zeta_out[ind_map] = -1;
#	Line 1032 | Line 1128 | void findCoord(struct mapInfo *mInfo)
1128		}
1129
1130
1131	<	/*
1131	>	/*
1132		* Sampling function
1133		* oversampling by nbin, then binning using a Gaussian
1134		* save results in outData, always of float type
1135		*
1136		*/
1137
1138	<	int performSampling(float *outData, float *inData, struct mapInfo *mInfo)
1138	>	int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
1139		{
1140
1141		int status = 0;
1142	+	int ind_map;
1143
1144		int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN
1145		int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
1146
1147	<	float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1147	>	// Changed Aug 12 2013, XS, for bitmaps
1148	>	float *outData0;
1149	>	if (interpOpt == 3 && mInfo->nbin == 1) {
1150	>	outData0 = outData;
1151	>	} else {
1152	>	outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1153	>	}
1154
1155		float *xi_out = mInfo->xi_out;
1156		float *zeta_out = mInfo->zeta_out;
1157	<
1157	>
1158		// Interpolation
1159
1160		struct fint_struct pars;
1161	<	int interpOpt = INTERP;         // Use Wiener by default, 6 order, 1 constraint
1161	>	// Aug 12 2013, passed in as argument now
1162
1163		switch (interpOpt) {
1164		case 0:                 // Wiener, 6 order, 1 constraint
#	Line 1067 | Line 1170 | int performSampling(float *outData, floa
1170		case 2:                 // Bilinear
1171		init_finterpolate_linear(&pars, 1.);
1172		break;
1173	+	case 3:                 // Near neighbor
1174	+	break;
1175		default:
1176		return 1;
1177		}
1178
1179	<	finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1179	>	printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
1180	>	if (interpOpt == 3) {                   // Aug 6 2013, Xudong
1181	>	for (int row0 = 0; row0 < nrow0; row0++) {
1182	>	for (int col0 = 0; col0 < ncol0; col0++) {
1183	>	ind_map = row0 * ncol0 + col0;
1184	>	outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
1185	>	}
1186	>	}
1187	>	} else {
1188	>	finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1189		FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1190	+	}
1191
1192		// Rebinning, smoothing
1193
1194	<	frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
1194	>	if (interpOpt == 3 && mInfo->nbin == 1) {
1195	>	return 0;
1196	>	} else {
1197	>	frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
1198	>	free(outData0);               // Dec 18 2012
1199	>	}
1200
1201		//
1202
#	Line 1085 | Line 1205 | int performSampling(float *outData, floa
1205		}
1206
1207
1208	<	/*
1208	>	/*
1209		* Performing local vector transformation
1210		*  xyz: z refers to vertical (radial) component, x EW (phi), y NS
1211		*
#	Line 1146 | Line 1266 | void vectorTransform(float *bx_map, floa
1266
1267		}
1268		}
1269	<
1269	>
1270		}
1271
1272
1273
1274	<	/*
1274	>	/*
1275		* Map and propogate vector field errors
1276		*
1277		*/
1278
1279		int getBErr(float *bx_err, float *by_err, float *bz_err,
1280	<	DRMS_Record_t *inRec, struct mapInfo *mInfo)
1280	>	DRMS_Record_t *inRec, struct mapInfo *mInfo)
1281		{
1282
1283		int status = 0;
#	Line 1176 | Line 1296 | int getBErr(float *bx_err, float *by_err
1296		float *errbTbA = (float *) (malloc(FOURK2 * sizeof(float)));
1297		float *errbIbA = (float *) (malloc(FOURK2 * sizeof(float)));
1298
1299	<	if (readVectorBErr(inRec,
1299	>	if (readVectorBErr(inRec,
1300		bT, bI, bA,
1301	<	errbT, errbI, errbA,
1301	>	errbT, errbI, errbA,
1302		errbTbI, errbTbA, errbIbA)) {
1303		printf("Read full disk variances & covariances error\n");
1304		free(bT); free(bI); free(bA);
#	Line 1215 | Line 1335 | int getBErr(float *bx_err, float *by_err
1335		int ind_map, ind_img;
1336
1337		double bpSigma2, btSigma2, brSigma2;            // variances after prop
1338	<
1338	>
1339		for (int row = 0; row < mInfo->nrow; row++) {
1340		for (int col = 0; col < mInfo->ncol; col++) {
1341
#	Line 1231 | Line 1351 | int getBErr(float *bx_err, float *by_err
1351		continue;
1352		}
1353
1354	<	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1354	>	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1355		disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1356		bx_err[ind_map] = DRMS_MISSING_FLOAT;
1357		bx_err[ind_map] = DRMS_MISSING_FLOAT;
#	Line 1244 | Line 1364 | int getBErr(float *bx_err, float *by_err
1364
1365		ind_img = round(zeta * FOURK + xi);
1366
1367	<	if (errorprop(bT, bA, bI,
1368	<	errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1369	<	lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1367	>	if (errorprop(bT, bA, bI,
1368	>	errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1369	>	lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1370		&btSigma2, &bpSigma2, &brSigma2)) {
1371		bx_err[ind_map] = DRMS_MISSING_FLOAT;
1372		by_err[ind_map] = DRMS_MISSING_FLOAT;
#	Line 1285 | Line 1405 | int readVectorB(DRMS_Record_t *inRec, fl
1405
1406		DRMS_Segment_t *inSeg;
1407		DRMS_Array_t *inArray_ambig;
1408	<	DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1408	>	DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1409
1410		char *ambig;
1411		float *bTotal, *bAzim, *bIncl;
#	Line 1348 | Line 1468 | int readVectorB(DRMS_Record_t *inRec, fl
1468		kOff = ky * bmx + kx;
1469		if (ambig[kOff] % 2) {          // 180
1470		bx_img[iData] *= -1.; by_img[iData] *= -1.;
1471	<	}
1471	>	}
1472		}
1473		}
1474		}
#	Line 1370 | Line 1490 | int readVectorB(DRMS_Record_t *inRec, fl
1490		*
1491		*/
1492
1493	<	int readVectorBErr(DRMS_Record_t *inRec,
1493	>	int readVectorBErr(DRMS_Record_t *inRec,
1494		float *bT, float *bI, float *bA,
1495	<	float *errbT, float *errbI, float *errbA,
1495	>	float *errbT, float *errbI, float *errbA,
1496		float *errbTbI, float *errbTbA, float *errbIbA)
1497		{
1498
#	Line 1380 | Line 1500 | int readVectorBErr(DRMS_Record_t *inRec,
1500
1501		float *data_ptr[9];
1502		char *segName[9] = {"field", "inclination", "azimuth",
1503	<	"field_err", "inclination_err", "azimuth_err",
1504	<	"field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1503	>	"field_err", "inclination_err", "azimuth_err",
1504	>	"field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1505		DRMS_Segment_t *inSegs[9];
1506		DRMS_Array_t *inArrays[9];
1507
#	Line 1417 | Line 1537 | int readVectorBErr(DRMS_Record_t *inRec,
1537		errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG;
1538		errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG;
1539		errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1540	<
1540	>
1541		}
1542
1543		//
1544
1545		for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
1546	<
1546	>
1547		return 0;
1548
1549		}
#	Line 1433 | Line 1553 | int readVectorBErr(DRMS_Record_t *inRec,
1553		* Create Cutout record: top level subroutine
1554		* Do the loops on segments and set the keywords here
1555		* Work is done in writeCutout routine below
1556	<	*
1556	>	*
1557		*/
1558
1559	<	int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1560	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1559	>	int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1560	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1561		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
1562		{
1563
#	Line 1454 | Line 1574 | int createCutRecord(DRMS_Record_t *mharp
1574		break;
1575		}
1576		}
1577	<	if (iHarpSeg != nMharpSegs) return 1;           // if failed
1577	>	if (iHarpSeg != nMharpSegs) {
1578	>	SHOW("Cutout: segment number unmatch\n");
1579	>	return 1;               // if failed
1580	>	}
1581		printf("Magnetogram cutout done.\n");
1582
1583		// Cutout Doppler
#	Line 1495 | Line 1618 | int createCutRecord(DRMS_Record_t *mharp
1618		if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
1619
1620		setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
1621	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
1622	<
1621	>	setKeys(sharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout
1622	>
1623		// Stats
1624	<
1624	>
1625		int nCutSegs = ARRLENGTH(CutSegs);
1626		for (int iSeg = 0; iSeg < nCutSegs; iSeg++) {
1627		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
#	Line 1506 | Line 1629 | int createCutRecord(DRMS_Record_t *mharp
1629		set_statistics(outSeg, outArray, 1);
1630		drms_free_array(outArray);
1631		}
1632	<
1632	>
1633		return 0;
1634
1635	<	}
1635	>	}
1636
1637
1638	<	/*
1638	>	/*
1639		* Get cutout and write segment
1640		* Change DISAMB_AZI to apply disambiguation to azimuth
1641		*
#	Line 1551 | Line 1674 | int writeCutout(DRMS_Record_t *outRec, D
1674		} else {
1675		return 1;
1676		}
1677	<
1677	>
1678		/* Adding disambiguation resolution to cutout azimuth? */
1679	<
1679	>
1680		#if DISAMB_AZI
1681		if (!strcmp(SegName, "azimuth")) {
1682	<	DRMS_Segment_t *disambSeg = drms_segment_lookup(inRec, "disambig");
1682	>	DRMS_Segment_t *disambSeg = NULL;
1683	>	disambSeg = drms_segment_lookup(inRec, "disambig");
1684		if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
1685		DRMS_Array_t *disambArray;
1686		if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1687		disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1688		if (status) {drms_free_array(cutoutArray); return 1;}
1689		} else {
1690	+	drms_free_array(cutoutArray);
1691		return 1;
1692		}
1693		double *azimuth = (double *) cutoutArray->data;
#	Line 1573 | Line 1698 | int writeCutout(DRMS_Record_t *outRec, D
1698		drms_free_array(disambArray);
1699		}
1700		#endif
1701	<
1701	>
1702		/* Write out */
1703
1704		outSeg = drms_segment_lookup(outRec, SegName);
1705		if (!outSeg) return 1;
1706	<	drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);
1706	>	//  drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);      // Jan 02 2013
1707		outSeg->axis[0] = cutoutArray->axis[0];
1708		outSeg->axis[1] = cutoutArray->axis[1];
1709	<	cutoutArray->parent_segment = outSeg;
1709	>	//  cutoutArray->parent_segment = outSeg;
1710		cutoutArray->israw = 0;         // always compressed
1711		cutoutArray->bzero = outSeg->bzero;
1712		cutoutArray->bscale = outSeg->bscale;               // Same as inArray's
#	Line 1594 | Line 1719 | int writeCutout(DRMS_Record_t *outRec, D
1719		}
1720
1721
1722	<	/*
1722	>	/*
1723		* Compute space weather indices, no error checking for now
1724		* Based on M. Bobra's swharp_vectorB.c
1725		* No error checking for now
#	Line 1608 | Line 1733 | void computeSWIndex(struct swIndex *swKe
1733		int nx = mInfo->ncol, ny = mInfo->nrow;
1734		int nxny = nx * ny;
1735		int dims[2] = {nx, ny};
1736	+
1737		// Get bx, by, bz, mask
1738
1739	<	// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1740	<	//DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "bitmap");
1741	<	//DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1742	<	//int *mask = (int *) maskArray->data;          // get the previously made mask array
1743	<
1744	<	//Use conf_disambig map as a threshold on spaceweather quantities
1745	<	DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1739	>	// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1740	>	DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
1741	>	DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
1742	>	int *bitmask = (int *) bitmaskArray->data;              // get the previously made mask array
1743	>
1744	>	//Use conf_disambig map as a threshold on spaceweather quantities
1745	>	DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1746		DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1747		int *mask = (int *) maskArray->data;            // get the previously made mask array
1748	<
1748	>
1749		DRMS_Segment_t *bxSeg = drms_segment_lookup(inRec, BP_SEG_CEA);
1750		DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
1751		float *bx = (float *) bxArray->data;            // bx
#	Line 1632 | Line 1758 | void computeSWIndex(struct swIndex *swKe
1758		DRMS_Segment_t *bzSeg = drms_segment_lookup(inRec, BR_SEG_CEA);
1759		DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
1760		float *bz = (float *) bzArray->data;            // bz
1761	+
1762	+	DRMS_Segment_t *bz_errSeg = drms_segment_lookup(inRec, BR_ERR_SEG_CEA);
1763	+	DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
1764	+	float *bz_err = (float *) bz_errArray->data;            // bz_err
1765	+
1766	+	DRMS_Segment_t *by_errSeg = drms_segment_lookup(inRec, BT_ERR_SEG_CEA);
1767	+	DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
1768	+	float *by_err = (float *) by_errArray->data;            // by_err
1769	+	//for (int i = 0; i < nxny; i++) by_err[i] *= -1;
1770	+
1771	+	DRMS_Segment_t *bx_errSeg = drms_segment_lookup(inRec, BP_ERR_SEG_CEA);
1772	+	DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
1773	+	float *bx_err = (float *) bx_errArray->data;            // bx_err
1774
1775		// Get emphemeris
1776	<
1777	<	//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1778	<	float cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);
1779	<	float dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1780	<	double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);
1781	<	double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);
1782	<	float imcrpix1    = drms_getkey_float(inRec, "IMCRPIX1", &status);
1783	<	float imcrpix2    = drms_getkey_float(inRec, "IMCRPIX2", &status);
1784	<	float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);
1785	<	float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);
1786	<
1787	<	//float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)
1788	<
1789	<	printf("cdelt1=%f\n",cdelt1);
1790	<	printf("rsun_ref=%f\n",rsun_ref);
1791	<	printf("rsun_obs=%f\n",rsun_obs);
1792	<	printf("dsun_obs=%f\n",dsun_obs);
1793	<
1794	<	// Temp arrays
1795	<
1796	<	float *bh = (float *) (malloc(nxny * sizeof(float)));
1797	<	float *bt = (float *) (malloc(nxny * sizeof(float)));
1659	<	float *jz = (float *) (malloc(nxny * sizeof(float)));
1660	<	float *bpx = (float *) (malloc(nxny * sizeof(float)));
1661	<	float *bpy = (float *) (malloc(nxny * sizeof(float)));
1662	<	float *bpz = (float *) (malloc(nxny * sizeof(float)));
1663	<	float *derx = (float *) (malloc(nxny * sizeof(float)));
1664	<	float *dery = (float *) (malloc(nxny * sizeof(float)));
1776	>	float  cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1777	>	float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1778	>	double rsun_ref    = drms_getkey_double(inRec, "RSUN_REF", &status);
1779	>	double rsun_obs    = drms_getkey_double(inRec, "RSUN_OBS", &status);
1780	>	float imcrpix1     = drms_getkey_float(inRec, "IMCRPIX1", &status);
1781	>	float imcrpix2     = drms_getkey_float(inRec, "IMCRPIX2", &status);
1782	>	float crpix1       = drms_getkey_float(inRec, "CRPIX1", &status);
1783	>	float crpix2       = drms_getkey_float(inRec, "CRPIX2", &status);
1784	>
1785	>	// convert cdelt1_orig from degrees to arcsec
1786	>	float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
1787	>
1788	>	// Temp arrays
1789	>	float *bh      = (float *) (malloc(nxny * sizeof(float)));
1790	>	float *bt      = (float *) (malloc(nxny * sizeof(float)));
1791	>	float *jz      = (float *) (malloc(nxny * sizeof(float)));
1792	>	float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
1793	>	float *bpx     = (float *) (malloc(nxny * sizeof(float)));
1794	>	float *bpy     = (float *) (malloc(nxny * sizeof(float)));
1795	>	float *bpz     = (float *) (malloc(nxny * sizeof(float)));
1796	>	float *derx    = (float *) (malloc(nxny * sizeof(float)));
1797	>	float *dery    = (float *) (malloc(nxny * sizeof(float)));
1798		float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
1799		float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
1800		float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
1801		float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
1802		float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1803		float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1804	<
1805	<	// Compute
1806	<
1807	<	if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1808	<	mask, cdelt1, rsun_ref, rsun_obs)){
1804	>	float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
1805	>	float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
1806	>	float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
1807	>	float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
1808	>	float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
1809	>	float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
1810	>	//spaceweather quantities computed
1811	>
1812	>
1813	>	if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),  &(swKeys_ptr->mean_vf_err),
1814	>	&(swKeys_ptr->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1815	>	{
1816		swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
1817		swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1818	+	swKeys_ptr->mean_vf_err = DRMS_MISSING_FLOAT;
1819	+	swKeys_ptr->count_mask  = DRMS_MISSING_INT;
1820		}
1821	<
1821	>
1822		for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
1823	<	greenpot(bpx, bpy, bpz, nx, ny);
1823	>	greenpot(bpx, bpy, bpz, nx, ny);
1824
1825	<	computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask);
1826	<
1827	<	if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask))
1828	<	swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
1825	>	computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
1826	>
1827	>	if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), &(swKeys_ptr->mean_gamma_err),mask, bitmask))
1828	>	{
1829	>	swKeys_ptr->mean_gamma     =  DRMS_MISSING_FLOAT;
1830	>	swKeys_ptr->mean_gamma_err =  DRMS_MISSING_FLOAT;
1831	>	}
1832
1833	<	computeB_total(bx, by, bz, bt, dims, mask);
1833	>	computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);
1834
1835	<	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt))
1835	>	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err)))
1836	>	{
1837		swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
1838	+	swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
1839	+	}
1840
1841	<	if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh))
1841	>	if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh), &(swKeys_ptr->mean_derivative_bh_err), mask, bitmask, derx_bh, dery_bh))
1842	>	{
1843		swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
1844	<
1845	<	if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, derx_bz, dery_bz))
1844	>	swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT;
1845	>	}
1846	>
1847	>	if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err), mask, bitmask, derx_bz, dery_bz))
1848	>	{
1849		swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1850	+	swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT;
1851	+	}
1852
1853	<
1854	<
1855	<	if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask,
1856	<	cdelt1, rsun_ref, rsun_obs, derx, dery)) {
1857	<	swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
1858	<	swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
1853	>	computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
1854	>	derx, dery);
1855	>
1856	>
1857	>	if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &(swKeys_ptr->mean_jz),
1858	>	&(swKeys_ptr->mean_jz_err), &(swKeys_ptr->us_i), &(swKeys_ptr->us_i_err), mask, bitmask, cdelt1,
1859	>	rsun_ref, rsun_obs, derx, dery))
1860	>	{
1861	>	swKeys_ptr->mean_jz            = DRMS_MISSING_FLOAT;
1862	>	swKeys_ptr->us_i               = DRMS_MISSING_FLOAT;
1863	>	swKeys_ptr->mean_jz_err        = DRMS_MISSING_FLOAT;
1864	>	swKeys_ptr->us_i_err           = DRMS_MISSING_FLOAT;
1865		}
1866	<
1867	<	printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
1868	<
1869	<	if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, cdelt1, rsun_ref, rsun_obs))
1870	<	swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
1871	<
1872	<	if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih),
1873	<	&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
1874	<	mask, cdelt1, rsun_ref, rsun_obs)) {
1875	<	swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
1876	<	swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
1877	<	swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
1866	>
1867	>	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))
1868	>	{
1869	>	swKeys_ptr->mean_alpha         = DRMS_MISSING_FLOAT;
1870	>	swKeys_ptr->mean_alpha_err     = DRMS_MISSING_FLOAT;
1871	>	}
1872	>
1873	>	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),
1874	>	&(swKeys_ptr->total_us_ih_err), &(swKeys_ptr->total_abs_ih_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1875	>	{
1876	>	swKeys_ptr->mean_ih            = DRMS_MISSING_FLOAT;
1877	>	swKeys_ptr->total_us_ih        = DRMS_MISSING_FLOAT;
1878	>	swKeys_ptr->total_abs_ih       = DRMS_MISSING_FLOAT;
1879	>	swKeys_ptr->mean_ih_err        = DRMS_MISSING_FLOAT;
1880	>	swKeys_ptr->total_us_ih_err    = DRMS_MISSING_FLOAT;
1881	>	swKeys_ptr->total_abs_ih_err   = DRMS_MISSING_FLOAT;
1882		}
1883	<
1884	<	if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),
1885	<	mask, cdelt1, rsun_ref, rsun_obs))
1886	<	swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
1887	<
1888	<
1725	<	if (computeFreeEnergy(bx, by, bpx, bpy, dims,
1726	<	&(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
1727	<	mask, cdelt1, rsun_ref, rsun_obs)) {
1728	<	swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1729	<	swKeys_ptr->totpot = DRMS_MISSING_FLOAT;
1883	>
1884	>	if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &(swKeys_ptr->totaljz), &(swKeys_ptr->totaljz_err),
1885	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1886	>	{
1887	>	swKeys_ptr->totaljz            = DRMS_MISSING_FLOAT;
1888	>	swKeys_ptr->totaljz_err        = DRMS_MISSING_FLOAT;
1889		}
1890
1891	<	if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,
1892	<	&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),
1893	<	&(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
1894	<	mask)) {
1895	<	swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1891	>	if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims,
1892	>	&(swKeys_ptr->meanpot), &(swKeys_ptr->meanpot_err), &(swKeys_ptr->totpot), &(swKeys_ptr->totpot_err),
1893	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1894	>	{
1895	>	swKeys_ptr->meanpot            = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1896	>	swKeys_ptr->totpot             = DRMS_MISSING_FLOAT;
1897	>	swKeys_ptr->meanpot_err        = DRMS_MISSING_FLOAT;
1898	>	swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;
1899	>	}
1900	>
1901	>	if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,
1902	>	&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),
1903	>	mask, bitmask)) {
1904	>	swKeys_ptr->meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1905		swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
1906	<	swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1739	<	swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;
1906	>	swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;
1907		}
1741	–
1742	–	// Clean up
1908
1909	+	// Clean up the arrays
1910	+
1911	+	drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong
1912		drms_free_array(maskArray);
1913		drms_free_array(bxArray);
1914		drms_free_array(byArray);
1915		drms_free_array(bzArray);
1916
1917	<	free(bh); free(bt); free(jz);
1917	>	free(bh); free(bt); free(jz); free(jz_smooth);
1918		free(bpx); free(bpy); free(bpz);
1919	<	free(derx); free(dery);
1920	<	free(derx_bt); free(dery_bt);
1921	<	free(derx_bz); free(dery_bz);
1919	>	free(derx); free(dery);
1920	>	free(derx_bt); free(dery_bt);
1921	>	free(derx_bz); free(dery_bz);
1922		free(derx_bh); free(dery_bh);
1923	<
1923	>	free(bt_err); free(bh_err);  free(jz_err);
1924	>	free(jz_err_squared); free(jz_rms_err);
1925	>	free(jz_err_squared_smooth);
1926		}
1927
1928	<
1759	<	/*
1928	>	/*
1929		* Set space weather indices, no error checking for now
1930		*
1931		*/
1932
1933		void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
1934		{
1935	<	drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
1935	>	drms_setkey_float(outRec, "USFLUX",  swKeys_ptr->mean_vf);
1936		drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
1937		drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
1938		drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
#	Line 1776 | Line 1945 | void setSWIndex(DRMS_Record_t *outRec, s
1945		drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
1946		drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
1947		drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
1948	<	drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);
1948	>	drms_setkey_float(outRec, "TOTPOT",  swKeys_ptr->totpot);
1949		drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
1950		drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
1951	+	drms_setkey_float(outRec, "CMASK",   swKeys_ptr->count_mask);
1952	+	drms_setkey_float(outRec, "ERRBT",   swKeys_ptr->mean_derivative_btotal_err);
1953	+	drms_setkey_float(outRec, "ERRVF",   swKeys_ptr->mean_vf_err);
1954	+	drms_setkey_float(outRec, "ERRGAM",  swKeys_ptr->mean_gamma_err);
1955	+	drms_setkey_float(outRec, "ERRBH",   swKeys_ptr->mean_derivative_bh_err);
1956	+	drms_setkey_float(outRec, "ERRBZ",   swKeys_ptr->mean_derivative_bz_err);
1957	+	drms_setkey_float(outRec, "ERRJZ",   swKeys_ptr->mean_jz_err);
1958	+	drms_setkey_float(outRec, "ERRUSI",  swKeys_ptr->us_i_err);
1959	+	drms_setkey_float(outRec, "ERRALP",  swKeys_ptr->mean_alpha_err);
1960	+	drms_setkey_float(outRec, "ERRMIH",  swKeys_ptr->mean_ih_err);
1961	+	drms_setkey_float(outRec, "ERRTUI",  swKeys_ptr->total_us_ih_err);
1962	+	drms_setkey_float(outRec, "ERRTAI",  swKeys_ptr->total_abs_ih_err);
1963	+	drms_setkey_float(outRec, "ERRJHT",  swKeys_ptr->totaljz_err);
1964	+	drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);
1965	+	drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);
1966	+	drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);
1967		};
1968
1969	<
1785	<	/*
1969	>	/*
1970		* Set all keywords, no error checking for now
1971		*
1972		*/
1973
1974	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)
1974	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo)
1975		{
1976	<	copy_me_keys(inRec, outRec);
1977	<	copy_patch_keys(inRec, outRec);
1978	<	copy_geo_keys(inRec, outRec);
1979	<	copy_ambig_keys(inRec, outRec);
1976	>	copy_me_keys(inRec, outRec);
1977	>	copy_patch_keys(inRec, outRec);
1978	>	copy_geo_keys(inRec, outRec);
1979	>	copy_ambig_keys(inRec, outRec);
1980	>
1981	>	int status = 0;
1982	>
1983	>	// Change a few geometry keywords for CEA records
1984	>	if (mInfo != NULL) {
1985	>
1986	>	drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
1987	>	drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
1988	>
1989	>	drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
1990	>	drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
1991	>	drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
1992	>	drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
1993	>	drms_setkey_string(outRec, "CUNIT1", "degree");
1994	>	drms_setkey_string(outRec, "CUNIT2", "degree");
1995	>
1996	>	char key[64];
1997	>	snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
1998	>	drms_setkey_string(outRec, "CTYPE1", key);
1999	>	snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
2000	>	drms_setkey_string(outRec, "CTYPE2", key);
2001	>	drms_setkey_float(outRec, "CROTA2", 0.0);
2002	>
2003	>	} else {
2004	>
2005	>	float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status);
2006	>	float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status);
2007	>	float x_ll = drms_getkey_float(inRec, "CRPIX1", &status);
2008	>	float y_ll = drms_getkey_float(inRec, "CRPIX2", &status);
2009	>	// Defined as disk center's pixel address wrt lower-left of cutout
2010	>	drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
2011	>	drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
2012	>	// Always 0.
2013	>	drms_setkey_float(outRec, "CRVAL1", 0);
2014	>	drms_setkey_float(outRec, "CRVAL2", 0);
2015	>
2016	>	}
2017	>
2018	>	char timebuf[1024];
2019	>	float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
2020	>	double val;
2021	>
2022	>	val = drms_getkey_double(inRec, "DATE",&status);
2023	>	drms_setkey_double(outRec, "DATE_B", val);
2024	>	sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
2025	>	drms_setkey_string(outRec, "DATE", timebuf);
2026	>
2027	>	// set cvs commit version into keyword HEADER
2028	>	char *cvsinfo = strdup("$Id$");
2029	>	char *cvsinfo2 = sw_functions_version();
2030	>	char cvsinfoall[2048];
2031	>	strcat(cvsinfoall,cvsinfo);
2032	>	strcat(cvsinfoall,"\n");
2033	>	strcat(cvsinfoall,cvsinfo2);
2034	>	status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
2035	>
2036		};
2037
2038		//

Diff Legend

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