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.8 by xudong, Mon Jan 7 22:03:47 2013 UTC vs.
Revision 1.28 by mbobra, Fri May 16 21:56:11 2014 UTC

#	Line 1 | Line 1
1		/*
2		*  sharp.c
3		*
4	<	*      This module creates the pipeline space weather harps
5	<	*      It is a hard-coded strip-down version of bmap.c
6	<	*      It takes the Mharp and Bharp series and crete the following quantities
7	<	*  Series 1: Sharp_CEA
8	<	*        CEA remapped magnetogram, bitmap, continuum, doppler (same size in map coordinate, need manual spec?)
9	<	*        CEA remapped vector field (Br, Bt, Bp) (same as above)
10	<	*    Space weather indices based on vector cutouts (step 2)
11	<	*  Series 2: Sharp_cutout:
12	<	*        cutouts of magnetogram, bitmap, continuum, doppler (HARP defined, various sizes in CCD pixels)
13	<	*        cutouts of all vector data segments (same as above)
14	<	*      Series 3: Other remaps
15	<	*
4	>	*      This module creates the pipeline Space Weather Active Region Patches (SHARPs).
5	>	*      It is a hard-coded strip-down version of bmap.c.
6	>	*      It takes the Mharp and Bharp series and create the following quantities:
7	>	*
8	>	*      Series 1: Sharp_CEA
9	>	*                CEA remapped magnetogram, bitmap, continuum, doppler (same size in map coordinate, need manual spec?)
10	>	*                CEA remapped vector field (Br, Bt, Bp) (same as above)
11	>	*                Space weather indices based on vector cutouts (step 2)
12	>	*
13	>	*      Series 2: Sharp_cutout:
14	>	*                cutouts of magnetogram, bitmap, continuum, doppler (HARP defined, various sizes in CCD pixels)
15	>	*                cutouts of all vector data segments (same as above)
16	>	*
17		*      Author:
18		*              Xudong Sun; Monica Bobra
19		*
20		*      Version:
21	<	*              v0.0    Jul 02 2012
22	<	*              v0.1    Jul 23 2012
23	<	*              v0.2    Sep 04 2012
24	<	*              v0.3    Dec 18 2012
25	<	*              v0.4    Jan 02 2013
21	>	*              v0.0 Jul 02 2012
22	>	*              v0.1 Jul 23 2012
23	>	*              v0.2 Sep 04 2012
24	>	*              v0.3 Dec 18 2012
25	>	*              v0.4 Jan 02 2013
26	>	*              v0.5 Jan 23 2013
27	>	*              v0.6 Aug 12 2013
28	>	*              v0.7 Jan 02 2014
29	>	*              v0.8 Feb 12 2014
30	>	*                              v0.9 Mar 04 2014
31		*
32		*      Notes:
33		*              v0.0
#	Line 39 | Line 45
45		*              Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"
46		*              v0.4
47		*              Took out convert_inplace(). Was causing all the images to be int
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"
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"
48	>	*              v0.5
49	>	*              Corrected ephemeris keywords, added argument mInfo for setKeys()
50	>	*              v0.6
51	>	*              Changes in remapping of bitmap and conf_disambig, now near neighbor without anti-aliasing
52	>	*              v0.7
53	>	*              Added full disk as "b"
54	>	*              Global flag fullDisk is set if "b" is set
55	>	*              Utilize BharpRS and BharpRec all around
56	>	*              Pass mharpRec to set_keys() too in case of full disk
57	>	*              Fixed Bunit (removed from copy_me_keys(), added loops for Bunits in set_keys() here)
58	>	*              Error for CEA still does account for disambiguation yet
59	>	*              v0.8
60	>	*              Added disambig to azimuth during error propagation
61	>	*              Changed usage for disambig: bit 2 (radial acute) for full disk, bit 0 for patch
62	>	*              Fixed disambig cutout for patch: 0 for even, 7 for odd
63	>	*              v0.9
64	>	*              Fixed unit
65	>	*              Check whether in PATCH of FD mode, so the error propagation uses disambiguated azimuth or not
66	>	*
67	>	*
68	>	*      Example Calls:
69	>	*      [I]   B (full disk disambiguation)
70	>	*      > 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"
71	>	*      [II]  BHARP (patch disambiguation)
72	>	*      > 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"
73		*
74		*
75		*/
#	Line 71 | Line 91
91		#include "errorprop.c"
92		#include "sw_functions.c"
93
94	+	//#include <mkl.h> // Comment out mkl.h, which can only run on solar3
95		#include <mkl_blas.h>
96		#include <mkl_service.h>
97		#include <mkl_lapack.h>
#	Line 89 | Line 110
110		// Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
111		#define NYQVIST         (0.015)
112
113	+	// Some other things
114		#ifndef MIN
115		#define MIN(a,b) (((a)<(b)) ? (a) : (b))
116		#endif
#	Line 116 | Line 138
138		#define INTERP                  0
139		#define dpath    "/home/jsoc/cvs/Development/JSOC"
140
141	+
142		/* ========================================================================================================== */
143
144		// Space weather keywords
145		struct swIndex {
146	<	float mean_vf;
146	>	float mean_vf;
147	>	float count_mask;
148		float absFlux;
149		float mean_hf;
150		float mean_gamma;
151		float mean_derivative_btotal;
152	<	float mean_derivative_bh;
152	>	float mean_derivative_bh;
153		float mean_derivative_bz;
154		float mean_jz;
155		float us_i;
156		float mean_alpha;
157	<	float mean_ih;
158	<	float total_us_ih;
159	<	float total_abs_ih;
160	<	float totaljz;
161	<	float totpot;
162	<	float meanpot;
163	<	float area_w_shear_gt_45;
164	<	float meanshear_angle;
165	<	float area_w_shear_gt_45h;
166	<	float meanshear_angleh;
157	>	float mean_ih;
158	>	float total_us_ih;
159	>	float total_abs_ih;
160	>	float totaljz;
161	>	float totpot;
162	>	float meanpot;
163	>	float area_w_shear_gt_45;
164	>	float meanshear_angle;
165	>	float area_w_shear_gt_45h;
166	>	float meanshear_angleh;
167	>	float mean_derivative_btotal_err;
168	>	float mean_vf_err;
169	>	float mean_gamma_err;
170	>	float mean_derivative_bh_err;
171	>	float mean_derivative_bz_err;
172	>	float mean_jz_err;
173	>	float us_i_err;
174	>	float mean_alpha_err;
175	>	float mean_ih_err;
176	>	float total_us_ih_err;
177	>	float total_abs_ih_err;
178	>	float totaljz_err;
179	>	float meanpot_err;
180	>	float totpot_err;
181	>	float meanshear_angle_err;
182	>	float Rparam;
183		};
184
185		// Mapping method
#	Line 156 | Line 196 | enum projection {
196		lambert
197		};
198
199	<	// Ephemeris
199	>	// WSC code
200	>	char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
201	>	"SIN", "ZEA"};
202	>
203	>	// Ephemeris information
204		struct ephemeris {
205		double disk_lonc, disk_latc;
206		double disk_xc, disk_yc;
#	Line 174 | Line 218 | struct mapInfo {
218		float *xi_out, *zeta_out;       // coordinate on full disk image to sample at
219		};
220
177	–
221		/* ========================================================================================================== */
222
223		/* Get all input data series */
#	Line 190 | Line 233 | int getInputRS_aux(DRMS_RecordSet_t **in
233		/* Find record from record set with given T_rec */
234		int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
235
193	–	// ===================
194	–
236		/* Create CEA record */
237		int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
238	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
238	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
239		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
240
241		/* Mapping single segment, wrapper */
#	Line 217 | Line 258 | int getEphemeris(DRMS_Record_t *inRec, s
258		void findCoord(struct mapInfo *mInfo);
259
260		/* Mapping function */
261	<	int performSampling(float *outData, float *inData, struct mapInfo *mInfo);
261	>	int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
262
263		/* Performing local vector transformation */
264		void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
#	Line 230 | Line 271 | int getBErr(float *bx_err, float *by_err
271		int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
272
273		/* Read variances and covariances of vector magnetograms */
274	<	int readVectorBErr(DRMS_Record_t *bharpRec,
274	>	int readVectorBErr(DRMS_Record_t *bharpRec,
275		float *bT, float *bI, float *bA,
276	<	float *errbT, float *errbI, float *errbA,
276	>	float *errbT, float *errbI, float *errbA,
277		float *errbTbI, float *errbTbA, float *errbIbA);
278
279		// ===================
#	Line 254 | Line 295 | void computeSWIndex(struct swIndex *swKe
295		void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
296
297		/* Set all keywords, no error checking for now */
298	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);
298	>	// Changed Dec 30 XS
299	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
300
301		// ===================
302
#	Line 285 | Line 327 | char *BharpSegs[] = {"inclination", "azi
327		"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
328		"disambig", "conf_disambig"};
329		// For stats
330	<	char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
330	>	char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
331		"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
332		"damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
333		"conv_flag", // fixed
#	Line 294 | Line 336 | char *CutSegs[] = {"magnetogram", "bitma
336		"field_inclination_err", "field_az_err", "inclin_azimuth_err",
337		"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
338		"disambig", "conf_disambig"};
339	<	char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",
340	<	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
341	<	/* ========================================================================================================== */
342	<
339	>	char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
340	>	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA, "conf_disambig"};
341	>	// For Bunits, added Dec 30 XS
342	>	char *CutBunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
343	>	"degree", "degree", "Mx/cm^2", "cm/s", "mA", " ",
344	>	"length units", "DN/s", "DN/s", " ", " ",
345	>	" ",
346	>	" ", " ",
347	>	"degree", "degree", "Mx/cm^2", "cm/s", " ",
348	>	" ", " ", " ",
349	>	" ", " ", " ",
350	>	" ", " "};
351	>	char *CEABunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
352	>	"Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", " "};      // Mar 4 2014 XS
353
354	+	/* ========================================================================================================== */
355
356		char *module_name = "sharp";
357	<	char *version_id = "2012 Dec 18";  /* Version number */
357	>	int seed;
358	>
359	>	int fullDisk;       // full disk mode
360	>	int amb4err;      // Use azimuth disambiguation for error propagation, default is 0 for patch and 1 for FD
361
362		ModuleArgs_t module_args[] =
363		{
364		{ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
365		{ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},
366	+	{ARG_STRING, "b", kNotSpecified, "Input B series, if set, overrides bharp."},
367		{ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},
368		{ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
369		{ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
370		{ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
371	+	{ARG_INT,    "seed", "987654", "Seed for the random number generator."},
372	+	{ARG_INT,   "f_amb4err", "0", "Force using disambiguation in error propagation"},     // Mar 4 2014 XS
373		{ARG_END}
374		};
375
376	<	int DoIt(void)
376	>	int DoIt(void)
377		{
378	<
378	>	int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
379	>	int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
380	>
381		int status = DRMS_SUCCESS;
382		int nrecs, irec;
383
384	<	char *mharpQuery, *bharpQuery;
384	>	char *mharpQuery, *bharpQuery, *bQuery;
385		char *dopQuery, *contQuery;
386		char *sharpCeaQuery, *sharpCutQuery;
387
#	Line 331 | Line 392 | int DoIt(void)
392
393		mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
394		bharpQuery = (char *) params_get_str(&cmdparams, "bharp");
395	+	bQuery = (char *) params_get_str(&cmdparams, "b");
396		dopQuery = (char *) params_get_str(&cmdparams, "dop");
397		contQuery = (char *) params_get_str(&cmdparams, "cont");
398		sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
399		sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
400	+
401	+	seed = params_get_int(&cmdparams, "seed");
402	+	int f_amb4err = params_get_int(&cmdparams, "f_amb4err");
403
404		/* Get input data, check everything */
405
406	<	if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
407	<	DIE("Input harp data error.");
406	>	// Full disk mode if "b" is set
407	>	if (strcmp(bQuery, kNotSpecified)) {
408	>	fullDisk = 1;
409	>	bharpQuery = bQuery;
410	>	//        SHOW(bharpQuery); SHOW("\n");
411	>	SHOW("Full disk mode\n");
412	>	} else {
413	>	fullDisk = 0;
414	>	SHOW("Harp mode\n");
415	>	}
416	>
417	>	// Mar 4 2014
418	>	if (f_amb4err == 0) {         // no forcing, 0 for patch and 1 for FD
419	>	amb4err = fullDisk ? 1 : 0;
420	>	} else {
421	>	amb4err = 1;
422	>	}
423	>	printf("amb4err=%d\n", amb4err);
424	>
425	>	// Bharp point to B if full disk
426	>	if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
427	>	DIE("Input harp data error.");
428		nrecs = mharpRS->n;
429
430	<	if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
430	>	if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
431		DIE("Input doppler data error.");
432
433	<	if (getInputRS_aux(&contRS, contQuery, mharpRS))
433	>	if (getInputRS_aux(&contRS, contQuery, mharpRS))
434		DIE("Input continuum data error.");
435
436		/* Start */
437
438		printf("==============\nStart. %d image(s) in total.\n", nrecs);
439	<
439	>
440		for (irec = 0; irec < nrecs; irec++) {
441
442		/* Records in work */
443
444		DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;
445	+
446		mharpRec = mharpRS->records[irec];
447	<	bharpRec = bharpRS->records[irec];
448	<
449	<	TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
450	<
447	>
448	>	TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
449	>
450	>	if (!fullDisk) {
451	>	bharpRec = bharpRS->records[irec];
452	>	} else {
453	>	if (getInputRec_aux(&bharpRec, bharpRS, trec)) {     // Bharp point to full disk B
454	>	printf("Fetching B failed, image #%d skipped.\n", irec);
455	>	continue;
456	>	}
457	>	}
458	>
459		struct swIndex swKeys;
460
461		DRMS_Record_t *dopRec = NULL, *contRec = NULL;
462	+
463		if (getInputRec_aux(&dopRec, dopRS, trec)) {
464		printf("Fetching Doppler failed, image #%d skipped.\n", irec);
465		continue;
#	Line 373 | Line 468 | int DoIt(void)
468		printf("Fetching continuum failed, image #%d skipped.\n", irec);
469		continue;
470		}
471	<
471	>
472		/* Create CEA record */
473
474		DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
#	Line 411 | Line 506 | int DoIt(void)
506		printf("Image #%d done.\n", irec);
507
508		} // irec
509	<
509	>
510
511		drms_close_records(mharpRS, DRMS_FREE_RECORD);
512		drms_close_records(bharpRS, DRMS_FREE_RECORD);
#	Line 443 | Line 538 | int getInputRS(DRMS_RecordSet_t **mharpR
538		*mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
539		if (status || (*mharpRS_ptr)->n == 0) return 1;
540
541	<	*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
542	<	if (status || (*bharpRS_ptr)->n == 0) return 1;
543	<
544	<	if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
541	>	if (fullDisk) {
542	>	if (getInputRS_aux(bharpRS_ptr, bharpQuery, *mharpRS_ptr)) return 1;
543	>	} else {
544	>	*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
545	>	if (status || (*bharpRS_ptr)->n == 0) return 1;
546	>	if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
547	>	}
548
549		return 0;
550
#	Line 472 | Line 570 | int compareHarp(DRMS_RecordSet_t *mharpR
570		for (int i = 0; i < nrecs; i++) {
571		mharpRec_t = mharpRS->records[i];
572		bharpRec_t = bharpRS->records[i];
573	<	if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
573	>	if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
574		drms_getkey_int(bharpRec_t, "HARPNUM", &status)) ||
575	<	(drms_getkey_time(mharpRec_t, "T_REC", &status) !=
575	>	(drms_getkey_time(mharpRec_t, "T_REC", &status) !=
576		drms_getkey_time(bharpRec_t, "T_REC", &status)))
577		{
578		return 1;
#	Line 485 | Line 583 | int compareHarp(DRMS_RecordSet_t *mharpR
583
584		}
585
586	<	/*
586	>	/*
587		* Get other data series, check all T_REC are available
588		*
589		*/
#	Line 529 | Line 627 | int getInputRS_aux(DRMS_RecordSet_t **in
627
628		}
629
630	<	/*
630	>	/*
631		* Find record from record set with given T_rec
632		*
633		*/
#	Line 559 | Line 657 | int getInputRec_aux(DRMS_Record_t **inRe
657		*
658		*/
659
660	<	int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
661	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
660	>	int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
661	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
662		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
663		{
664
#	Line 572 | Line 670 | int createCeaRecord(DRMS_Record_t *mharp
670		mInfo.proj = (enum projection) cyleqa;          // projection method
671		mInfo.xscale = XSCALE;
672		mInfo.yscale = YSCALE;
673	<	mInfo.nbin = NBIN;
673	>
674	>	int ncol0, nrow0;           // oversampled map size
675
676		// Get ephemeris
677
#	Line 588 | Line 687 | int createCeaRecord(DRMS_Record_t *mharp
687		return 1;
688		}
689
690	+	// ========================================
691	+	// Do this for all bitmaps, Aug 12 2013 XS
692	+	// ========================================
693	+
694	+	mInfo.nbin = 1;                     // for bitmaps. suppress anti-aliasing
695	+	ncol0 = mInfo.ncol;
696	+	nrow0 = mInfo.nrow;
697	+
698	+	mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
699	+	mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
700	+
701	+	findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
702	+
703	+	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
704	+	SHOW("CEA: mapping bitmap error\n");
705	+	return 1;
706	+	}
707	+	printf("Bitmap mapping done.\n");
708	+
709	+	if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
710	+	SHOW("CEA: mapping conf_disambig error\n");
711	+	return 1;
712	+	}
713	+	printf("Conf disambig mapping done.\n");
714	+
715	+	free(mInfo.xi_out);
716	+	free(mInfo.zeta_out);
717	+
718	+	// ========================================
719	+	// Do this again for floats, Aug 12 2013 XS
720	+	// ========================================
721		// Create xi_out, zeta_out array in mInfo:
722		// Coordinates to sample in original full disk image
723
724	<	int ncol0, nrow0;               // oversampled map size
724	>	mInfo.nbin = NBIN;
725		ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
726		nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
727
#	Line 601 | Line 731 | int createCeaRecord(DRMS_Record_t *mharp
731		findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
732
733		// Mapping single segment: Mharp, etc.
734	<
734	>
735		if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
736		SHOW("CEA: mapping magnetogram error\n");
737		return 1;
738		}
609	–	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
610	–	SHOW("CEA: mapping bitmap error\n");
611	–	return 1;
612	–	}
739		printf("Magnetogram mapping done.\n");
740	<
740	>
741		if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
742		SHOW("CEA: mapping dopplergram error\n");
743		return 1;
#	Line 624 | Line 750 | int createCeaRecord(DRMS_Record_t *mharp
750		}
751		printf("Intensitygram mapping done.\n");
752
627	–	if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
628	–	SHOW("CEA: mapping conf_disambig error\n");
629	–	return 1;
630	–	}
631	–	printf("Conf disambig mapping done.\n");
632	–
753		// Mapping vector B
754
755		if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
#	Line 651 | Line 771 | int createCeaRecord(DRMS_Record_t *mharp
771		drms_copykey(sharpRec, mharpRec, "T_REC");
772		drms_copykey(sharpRec, mharpRec, "HARPNUM");
773
774	+	if (fullDisk) {
775	+	DRMS_Link_t *bLink = hcon_lookup_lower(&sharpRec->links, "B");
776	+	if (bLink) drms_link_set("B", sharpRec, bharpRec);
777	+	} else {
778	+	DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
779	+	if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
780	+	}
781		DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
782		if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
656	–	DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
657	–	if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
783
784	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
784	>	setKeys(sharpRec, mharpRec, bharpRec, &mInfo);            // Set all other keywords
785		drms_copykey(sharpRec, mharpRec, "QUALITY");            // copied from los records
786
787		// Space weather
#	Line 684 | Line 809 | int createCeaRecord(DRMS_Record_t *mharp
809		}
810
811
812	<	/*
812	>	/*
813		* Mapping a single segment
814		* Read in full disk image, utilize mapImage for mapping
815		* then write the segment out, segName same in in/out Rec
#	Line 698 | Line 823 | int mapScaler(DRMS_Record_t *sharpRec, D
823		int status = 0;
824		int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
825		int dims[2] = {nx, ny};
826	+	int interpOpt = INTERP;         // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
827
828		// Input full disk array
829
#	Line 709 | Line 835 | int mapScaler(DRMS_Record_t *sharpRec, D
835		inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
836		if (!inArray) return 1;
837
838	+	if (!strcmp(segName, "conf_disambig") || !strcmp(segName, "bitmap")) {
839	+	// Moved out so it works for FD conf_disambig as well
840	+	// Jan 2 2014 XS
841	+	interpOpt = 3;          // Aug 12 XS, near neighbor
842	+	}
843	+
844		float *inData;
845		int xsz = inArray->axis[0], ysz = inArray->axis[1];
846		if ((xsz != FOURK) || (ysz != FOURK)) {         // for bitmap, make tmp full disk
#	Line 731 | Line 863 | int mapScaler(DRMS_Record_t *sharpRec, D
863		// Mapping
864
865		float *map = (float *) (malloc(nxny * sizeof(float)));
866	<	if (performSampling(map, inData, mInfo))
866	>	if (performSampling(map, inData, mInfo, interpOpt))             // Add interpOpt for different types, Aug 12 XS
867		{if (inArray) drms_free_array(inArray); free(map); return 1;}
868
869		// Write out
870
871	<	DRMS_Segment_t *outSeg = NULL;
871	>	DRMS_Segment_t *outSeg = NULL;
872		outSeg = drms_segment_lookup(sharpRec, segName);
873		if (!outSeg) return 1;
874
875	<	//      DRMS_Type_t arrayType = outSeg->info->type;
875	>	//  DRMS_Type_t arrayType = outSeg->info->type;
876		DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
877		if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
878
879		// convert to needed data type
880
881	<	//      drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);         // Jan 02 2013
881	>	//  drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);         // Jan 02 2013
882
883		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
884	<	//      outArray->parent_segment = outSeg;
884	>	//  outArray->parent_segment = outSeg;
885		outArray->israw = 0;            // always compressed
886		outArray->bzero = outSeg->bzero;
887		outArray->bscale = outSeg->bscale;
#	Line 765 | Line 897 | int mapScaler(DRMS_Record_t *sharpRec, D
897		}
898
899
900	<	/*
900	>	/*
901		* Mapping vector magnetogram
902		*
903		*/
#	Line 794 | Line 926 | int mapVectorB(DRMS_Record_t *sharpRec,
926		float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation
927
928		bx_map = (float *) (malloc(nxny * sizeof(float)));
929	<	if (performSampling(bx_map, bx_img, mInfo))
929	>	if (performSampling(bx_map, bx_img, mInfo, INTERP))
930		{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
931
932		by_map = (float *) (malloc(nxny * sizeof(float)));
933	<	if (performSampling(by_map, by_img, mInfo))
933	>	if (performSampling(by_map, by_img, mInfo, INTERP))
934		{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
935
936		bz_map = (float *) (malloc(nxny * sizeof(float)));
937	<	if (performSampling(bz_map, bz_img, mInfo))
937	>	if (performSampling(bz_map, bz_img, mInfo, INTERP))
938		{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
939
940		free(bx_img); free(by_img); free(bz_img);
#	Line 825 | Line 957 | int mapVectorB(DRMS_Record_t *sharpRec,
957		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
958		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
959		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
960	<	//              outArray->parent_segment = outSeg;
960	>	//              outArray->parent_segment = outSeg;
961		outArray->israw = 0;
962		outArray->bzero = outSeg->bzero;
963		outArray->bscale = outSeg->bscale;
#	Line 841 | Line 973 | int mapVectorB(DRMS_Record_t *sharpRec,
973		}
974
975
976	<	/*
976	>	/*
977		* Mapping vector magnetogram errors
978		*
979		*/
#	Line 877 | Line 1009 | int mapVectorBErr(DRMS_Record_t *sharpRe
1009		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
1010		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
1011		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
1012	<	//              outArray->parent_segment = outSeg;
1012	>	//              outArray->parent_segment = outSeg;
1013		outArray->israw = 0;
1014		outArray->bzero = outSeg->bzero;
1015		outArray->bscale = outSeg->bscale;
#	Line 894 | Line 1026 | int mapVectorBErr(DRMS_Record_t *sharpRe
1026
1027
1028
1029	<	/*
1029	>	/*
1030		* Determine reference point coordinate and patch size according to keywords
1031		* xc, yc are the coordinate of patch center, in degrees
1032		* ncol and nrow are the final size
#	Line 922 | Line 1054 | int findPosition(DRMS_Record_t *inRec, s
1054		// We compute minlon & minlat then by
1055		// LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
1056
1057	<	float psize = drms_getkey_float(inRec, "SIZE", &status);
1058	<	if (psize != psize) {
1059	<	TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1;
1057	>	//  float psize = drms_getkey_float(inRec, "SIZE", &status);
1058	>	//  if (psize != psize) {
1059	>
1060	>	if (minlon != minlon || maxlon != maxlon) {         // check lons instead of SIZE
1061	>	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
1062		double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
1063		char firstRecQuery[100], t0_str[100];
1064		sprint_time(t0_str, t0, "TAI", 0);
#	Line 964 | Line 1098 | int getEphemeris(DRMS_Record_t *inRec, s
1098		int status = 0;
1099
1100		float crota2 = drms_getkey_float(inRec, "CROTA2", &status);     // rotation
1101	<	double sina = sin(crota2 * RADSINDEG);
1101	>	double sina = sin(crota2 * RADSINDEG);
1102		double cosa = cos(crota2 * RADSINDEG);
1103
1104		ephem->pa = - crota2 * RADSINDEG;
#	Line 976 | Line 1110 | int getEphemeris(DRMS_Record_t *inRec, s
1110		float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
1111		float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
1112		float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
979	–	printf("cdelt=%f\n",cdelt);
1113		ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0
1114		ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
1115
#	Line 1041 | Line 1174 | void findCoord(struct mapInfo *mInfo)
1174		x = (col0 + 0.5 - ncol0/2.) * xscale0;          // in rad
1175		y = (row0 + 0.5 - nrow0/2.) * yscale0;
1176
1177	<	/* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1177	>	/* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1178		* the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
1179		* projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
1180	<	* is the heliographic longitude and latitude of the map center. Both are in degree.
1180	>	* is the heliographic longitude and latitude of the map center. Both are in degree.
1181		*/
1182
1183		if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
#	Line 1057 | Line 1190 | void findCoord(struct mapInfo *mInfo)
1190		* image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
1191		*/
1192
1193	<	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1193	>	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1194		disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1195		xi_out[ind_map] = -1;
1196		zeta_out[ind_map] = -1;
#	Line 1073 | Line 1206 | void findCoord(struct mapInfo *mInfo)
1206		}
1207
1208
1209	<	/*
1209	>	/*
1210		* Sampling function
1211		* oversampling by nbin, then binning using a Gaussian
1212		* save results in outData, always of float type
1213		*
1214		*/
1215
1216	<	int performSampling(float *outData, float *inData, struct mapInfo *mInfo)
1216	>	int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
1217		{
1218
1219		int status = 0;
1220	+	int ind_map;
1221
1222		int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN
1223		int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
1224
1225	<	float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1225	>	// Changed Aug 12 2013, XS, for bitmaps
1226	>	float *outData0;
1227	>	if (interpOpt == 3 && mInfo->nbin == 1) {
1228	>	outData0 = outData;
1229	>	} else {
1230	>	outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1231	>	}
1232
1233		float *xi_out = mInfo->xi_out;
1234		float *zeta_out = mInfo->zeta_out;
#	Line 1096 | Line 1236 | int performSampling(float *outData, floa
1236		// Interpolation
1237
1238		struct fint_struct pars;
1239	<	int interpOpt = INTERP;         // Use Wiener by default, 6 order, 1 constraint
1239	>	// Aug 12 2013, passed in as argument now
1240
1241		switch (interpOpt) {
1242		case 0:                 // Wiener, 6 order, 1 constraint
#	Line 1108 | Line 1248 | int performSampling(float *outData, floa
1248		case 2:                 // Bilinear
1249		init_finterpolate_linear(&pars, 1.);
1250		break;
1251	+	case 3:                 // Near neighbor
1252	+	break;
1253		default:
1254		return 1;
1255		}
1256
1257	<	finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1258	<	FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1257	>	printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
1258	>	if (interpOpt == 3) {                   // Aug 6 2013, Xudong
1259	>	for (int row0 = 0; row0 < nrow0; row0++) {
1260	>	for (int col0 = 0; col0 < ncol0; col0++) {
1261	>	ind_map = row0 * ncol0 + col0;
1262	>	outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
1263	>	}
1264	>	}
1265	>	} else {
1266	>	finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1267	>	FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1268	>	}
1269
1270		// Rebinning, smoothing
1271
1272	<	frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
1273	<	free(outData0);         // Dec 18 2012
1272	>	if (interpOpt == 3 && mInfo->nbin == 1) {
1273	>	return 0;
1274	>	} else {
1275	>	frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
1276	>	free(outData0);         // Dec 18 2012
1277	>	}
1278
1279		//
1280
#	Line 1127 | Line 1283 | int performSampling(float *outData, floa
1283		}
1284
1285
1286	<	/*
1286	>	/*
1287		* Performing local vector transformation
1288		*  xyz: z refers to vertical (radial) component, x EW (phi), y NS
1289		*
#	Line 1193 | Line 1349 | void vectorTransform(float *bx_map, floa
1349
1350
1351
1352	<	/*
1352	>	/*
1353		* Map and propogate vector field errors
1354		*
1355		*/
#	Line 1218 | Line 1374 | int getBErr(float *bx_err, float *by_err
1374		float *errbTbA = (float *) (malloc(FOURK2 * sizeof(float)));
1375		float *errbIbA = (float *) (malloc(FOURK2 * sizeof(float)));
1376
1377	<	if (readVectorBErr(inRec,
1377	>	if (readVectorBErr(inRec,
1378		bT, bI, bA,
1379	<	errbT, errbI, errbA,
1379	>	errbT, errbI, errbA,
1380		errbTbI, errbTbA, errbIbA)) {
1381		printf("Read full disk variances & covariances error\n");
1382		free(bT); free(bI); free(bA);
#	Line 1273 | Line 1429 | int getBErr(float *bx_err, float *by_err
1429		continue;
1430		}
1431
1432	<	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1432	>	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1433		disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1434		bx_err[ind_map] = DRMS_MISSING_FLOAT;
1435	<	bx_err[ind_map] = DRMS_MISSING_FLOAT;
1436	<	bx_err[ind_map] = DRMS_MISSING_FLOAT;
1435	>	by_err[ind_map] = DRMS_MISSING_FLOAT;
1436	>	bz_err[ind_map] = DRMS_MISSING_FLOAT;       // Mar 7
1437		continue;
1438		}
1439
#	Line 1286 | Line 1442 | int getBErr(float *bx_err, float *by_err
1442
1443		ind_img = round(zeta * FOURK + xi);
1444
1445	<	if (errorprop(bT, bA, bI,
1446	<	errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1447	<	lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1445	>	if (errorprop(bT, bA, bI,
1446	>	errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1447	>	lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1448		&btSigma2, &bpSigma2, &brSigma2)) {
1449		bx_err[ind_map] = DRMS_MISSING_FLOAT;
1450		by_err[ind_map] = DRMS_MISSING_FLOAT;
#	Line 1356 | Line 1512 | int readVectorB(DRMS_Record_t *inRec, fl
1512
1513		int llx, lly;           // lower-left corner
1514		int bmx, bmy;           // bitmap size
1515	<
1516	<	llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1517	<	lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1518	<
1519	<	bmx = inArray_ambig->axis[0];
1520	<	bmy = inArray_ambig->axis[1];
1515	>
1516	>	if (fullDisk) {
1517	>	llx = lly = 0;
1518	>	bmx = bmy = FOURK;
1519	>	} else {
1520	>	llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1521	>	lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1522	>	bmx = inArray_ambig->axis[0];
1523	>	bmy = inArray_ambig->axis[1];
1524	>	}
1525
1526		int kx, ky, kOff;
1527		int ix = 0, jy = 0, yOff = 0, iData = 0;
1528		int xDim = FOURK, yDim = FOURK;
1529	+	int amb = 0;
1530
1531		for (jy = 0; jy < yDim; jy++)
1532		{
#	Line 1388 | Line 1549 | int readVectorB(DRMS_Record_t *inRec, fl
1549		continue;
1550		} else {
1551		kOff = ky * bmx + kx;
1552	<	if (ambig[kOff] % 2) {          // 180
1552	>	//                              if (ambig[kOff] % 2) {          // 180
1553	>	// Feb 12 2014, use bit #2 for full disk, lowest bit for patch
1554	>	if (fullDisk) { amb = (ambig[kOff] / 4) % 2; } else { amb = ambig[kOff] % 2; }
1555	>	if (amb) {                              // Feb 12 2014, use bit #2
1556		bx_img[iData] *= -1.; by_img[iData] *= -1.;
1557	<	}
1557	>	}
1558		}
1559		}
1560		}
#	Line 1412 | Line 1576 | int readVectorB(DRMS_Record_t *inRec, fl
1576		*
1577		*/
1578
1579	<	int readVectorBErr(DRMS_Record_t *inRec,
1579	>	int readVectorBErr(DRMS_Record_t *inRec,
1580		float *bT, float *bI, float *bA,
1581	<	float *errbT, float *errbI, float *errbA,
1581	>	float *errbT, float *errbI, float *errbA,
1582		float *errbTbI, float *errbTbA, float *errbIbA)
1583		{
1584
#	Line 1428 | Line 1592 | int readVectorBErr(DRMS_Record_t *inRec,
1592		DRMS_Array_t *inArrays[9];
1593
1594		// Read full disk images
1595	+	// Do we need disambig? Dec 30 XS
1596
1597		for (int iSeg = 0; iSeg < 9; iSeg++) {
1598
#	Line 1441 | Line 1606 | int readVectorBErr(DRMS_Record_t *inRec,
1606		float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];
1607		float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];
1608
1609	+	// Add disambig, Feb 12 2014
1610	+
1611	+	DRMS_Segment_t *inSeg;
1612	+	DRMS_Array_t *inArray_ambig;
1613	+
1614	+	if (amb4err) {              // Mar 4 2014
1615	+
1616	+	inSeg = drms_segment_lookup(inRec, "disambig");
1617	+	inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
1618	+	if (status) return 1;
1619	+	char *ambig = (char *)inArray_ambig->data;
1620	+
1621	+	int llx, lly;           // lower-left corner
1622	+	int bmx, bmy;           // bitmap size
1623	+
1624	+	if (fullDisk) {
1625	+	llx = lly = 0;
1626	+	bmx = bmy = FOURK;
1627	+	} else {
1628	+	llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1629	+	lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1630	+	bmx = inArray_ambig->axis[0];
1631	+	bmy = inArray_ambig->axis[1];
1632	+	}
1633	+
1634	+	int idx, idx_a;
1635	+	int amb;
1636	+
1637	+	for (int j = 0; j < bmy; j++) {
1638	+	for (int i = 0; i < bmx; i++) {
1639	+	idx_a = j * bmx + i;
1640	+	idx = (j + lly) * FOURK + (i + llx);
1641	+	// Feb 12 2014, use bit #2 for full disk, lowest bit for patch
1642	+	if (fullDisk) { amb = (ambig[idx_a] / 4) % 2; } else { amb = ambig[idx_a] % 2; }
1643	+	if (amb) { bA0[idx] += 180.; }
1644	+	}
1645	+	}
1646	+
1647	+	}
1648	+
1649		// Convert errors to variances, correlation coefficients to covariances
1650
1651		for (int i = 0; i < FOURK2; i++) {
#	Line 1449 | Line 1654 | int readVectorBErr(DRMS_Record_t *inRec,
1654		if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;
1655
1656		bT[i] = bT0[i];
1657	<	bI[i] = bI0[i];
1657	>	bI[i] = bI0[i];         // in deg, coverted in errorprop
1658		bA[i] = bA0[i];
1659
1660		errbT[i] = errbT0[i] * errbT0[i];
#	Line 1465 | Line 1670 | int readVectorBErr(DRMS_Record_t *inRec,
1670		//
1671
1672		for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
1673	+	if (amb4err) drms_free_array(inArray_ambig);            // Feb 12; Mar 04 2014
1674
1675		return 0;
1676
#	Line 1475 | Line 1681 | int readVectorBErr(DRMS_Record_t *inRec,
1681		* Create Cutout record: top level subroutine
1682		* Do the loops on segments and set the keywords here
1683		* Work is done in writeCutout routine below
1684	<	*
1684	>	*
1685		*/
1686
1687	<	int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1688	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1687	>	int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1688	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1689		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
1690		{
1691
#	Line 1540 | Line 1746 | int createCutRecord(DRMS_Record_t *mharp
1746		if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
1747
1748		setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
1749	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
1749	>	setKeys(sharpRec, mharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout
1750
1751		// Stats
1752
#	Line 1554 | Line 1760 | int createCutRecord(DRMS_Record_t *mharp
1760
1761		return 0;
1762
1763	<	}
1763	>	}
1764
1765
1766	<	/*
1766	>	/*
1767		* Get cutout and write segment
1768		* Change DISAMB_AZI to apply disambiguation to azimuth
1769		*
#	Line 1597 | Line 1803 | int writeCutout(DRMS_Record_t *outRec, D
1803		return 1;
1804		}
1805
1806	+	// Feb 12 2014, fool-proof, for patch, change everything to 0 or 7!!!
1807	+	// This is a fix for disambiguation before Aug 2013
1808	+
1809	+	if (!strcmp(SegName, "disambig") && !fullDisk) {
1810	+	double *disamb = (double *) (cutoutArray->data);
1811	+	for (int i = 0; i < nxny; i++) {
1812	+	if (((int)disamb[i]) % 2) { disamb[i] = 7; } else { disamb[i] = 0; }
1813	+	}
1814	+	}
1815	+
1816		/* Adding disambiguation resolution to cutout azimuth? */
1817
1818		#if DISAMB_AZI
1819	+	int amb;
1820		if (!strcmp(SegName, "azimuth")) {
1821		DRMS_Segment_t *disambSeg = NULL;
1822		disambSeg = drms_segment_lookup(inRec, "disambig");
1823		if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
1824		DRMS_Array_t *disambArray;
1825	<	if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1826	<	disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1827	<	if (status) {drms_free_array(cutoutArray); return 1;}
1828	<	} else {
1829	<	drms_free_array(cutoutArray);
1830	<	return 1;
1831	<	}
1825	>	if (fullDisk) { // Jan 2 2014 XS
1826	>	disambArray = drms_segment_readslice(disambSeg, DRMS_TYPE_CHAR, ll, ur, &status);
1827	>	if (status) return 1;
1828	>	} else {
1829	>	if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1830	>	disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1831	>	if (status) {drms_free_array(cutoutArray); return 1;}
1832	>	} else {
1833	>	drms_free_array(cutoutArray);
1834	>	return 1;
1835	>	}
1836	>	}
1837		double *azimuth = (double *) cutoutArray->data;
1838		char *disamb = (char *) disambArray->data;
1839		for (int n = 0; n < nxny; n++) {
1840	<	if (disamb[n]) azimuth[n] += 180.;
1840	>	//                  if (disamb[n] % 2) azimuth[n] += 180.;      // Nov 12 2013 Fixed!!!
1841	>	// Feb 12 2014, use bit #2 for full disk, lowest bit for patch
1842	>	if (fullDisk) { amb = (disamb[n] / 4) % 2; } else { amb = disamb[n] % 2; }
1843	>	if (amb) azimuth[n] += 180.;
1844		}
1845		drms_free_array(disambArray);
1846		}
#	Line 1625 | Line 1850 | int writeCutout(DRMS_Record_t *outRec, D
1850
1851		outSeg = drms_segment_lookup(outRec, SegName);
1852		if (!outSeg) return 1;
1853	<	//      drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);      // Jan 02 2013
1853	>	//  drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);      // Jan 02 2013
1854		outSeg->axis[0] = cutoutArray->axis[0];
1855		outSeg->axis[1] = cutoutArray->axis[1];
1856	<	//      cutoutArray->parent_segment = outSeg;
1856	>	//  cutoutArray->parent_segment = outSeg;
1857		cutoutArray->israw = 0;         // always compressed
1858		cutoutArray->bzero = outSeg->bzero;
1859		cutoutArray->bscale = outSeg->bscale;               // Same as inArray's
#	Line 1641 | Line 1866 | int writeCutout(DRMS_Record_t *outRec, D
1866		}
1867
1868
1869	<	/*
1869	>	/*
1870		* Compute space weather indices, no error checking for now
1871		* Based on M. Bobra's swharp_vectorB.c
1872		* No error checking for now
#	Line 1655 | Line 1880 | void computeSWIndex(struct swIndex *swKe
1880		int nx = mInfo->ncol, ny = mInfo->nrow;
1881		int nxny = nx * ny;
1882		int dims[2] = {nx, ny};
1883	+
1884		// Get bx, by, bz, mask
1885
1886		// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
#	Line 1662 | Line 1888 | void computeSWIndex(struct swIndex *swKe
1888		DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
1889		int *bitmask = (int *) bitmaskArray->data;              // get the previously made mask array
1890
1891	<	//Use conf_disambig map as a threshold on spaceweather quantities
1892	<	DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1891	>	//Use conf_disambig map as a threshold on spaceweather quantities
1892	>	DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1893		DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1894		int *mask = (int *) maskArray->data;            // get the previously made mask array
1895
#	Line 1679 | Line 1905 | void computeSWIndex(struct swIndex *swKe
1905		DRMS_Segment_t *bzSeg = drms_segment_lookup(inRec, BR_SEG_CEA);
1906		DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
1907		float *bz = (float *) bzArray->data;            // bz
1908	+
1909	+	//Use magnetogram map to compute R
1910	+	DRMS_Segment_t *losSeg = drms_segment_lookup(inRec, "magnetogram");
1911	+	DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
1912	+	float *los = (float *) losArray->data;          // los
1913	+
1914	+	DRMS_Segment_t *bz_errSeg = drms_segment_lookup(inRec, BR_ERR_SEG_CEA);
1915	+	DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
1916	+	float *bz_err = (float *) bz_errArray->data;            // bz_err
1917	+
1918	+	DRMS_Segment_t *by_errSeg = drms_segment_lookup(inRec, BT_ERR_SEG_CEA);
1919	+	DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
1920	+	float *by_err = (float *) by_errArray->data;            // by_err
1921	+	//for (int i = 0; i < nxny; i++) by_err[i] *= -1;
1922	+
1923	+	DRMS_Segment_t *bx_errSeg = drms_segment_lookup(inRec, BP_ERR_SEG_CEA);
1924	+	DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
1925	+	float *bx_err = (float *) bx_errArray->data;            // bx_err
1926
1927		// Get emphemeris
1928	<
1929	<	//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1930	<	float cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);
1931	<	float dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1932	<	double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);
1933	<	double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);
1934	<	float imcrpix1    = drms_getkey_float(inRec, "IMCRPIX1", &status);
1935	<	float imcrpix2    = drms_getkey_float(inRec, "IMCRPIX2", &status);
1936	<	float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);
1937	<	float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);
1938	<
1939	<	//float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)
1940	<
1941	<	printf("cdelt1=%f\n",cdelt1);
1942	<	printf("rsun_ref=%f\n",rsun_ref);
1943	<	printf("rsun_obs=%f\n",rsun_obs);
1944	<	printf("dsun_obs=%f\n",dsun_obs);
1945	<
1946	<	// Temp arrays
1947	<
1948	<	float *bh = (float *) (malloc(nxny * sizeof(float)));
1949	<	float *bt = (float *) (malloc(nxny * sizeof(float)));
1950	<	float *jz = (float *) (malloc(nxny * sizeof(float)));
1951	<	float *bpx = (float *) (malloc(nxny * sizeof(float)));
1952	<	float *bpy = (float *) (malloc(nxny * sizeof(float)));
1953	<	float *bpz = (float *) (malloc(nxny * sizeof(float)));
1954	<	float *derx = (float *) (malloc(nxny * sizeof(float)));
1955	<	float *dery = (float *) (malloc(nxny * sizeof(float)));
1928	>	float  cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1929	>	float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1930	>	double rsun_ref    = drms_getkey_double(inRec, "RSUN_REF", &status);
1931	>	double rsun_obs    = drms_getkey_double(inRec, "RSUN_OBS", &status);
1932	>	float imcrpix1     = drms_getkey_float(inRec, "IMCRPIX1", &status);
1933	>	float imcrpix2     = drms_getkey_float(inRec, "IMCRPIX2", &status);
1934	>	float crpix1       = drms_getkey_float(inRec, "CRPIX1", &status);
1935	>	float crpix2       = drms_getkey_float(inRec, "CRPIX2", &status);
1936	>
1937	>	// convert cdelt1_orig from degrees to arcsec
1938	>	float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
1939	>
1940	>	// define some values for the R calculation
1941	>	int scale = round(2.0/cdelt1);
1942	>	int nx1 = nx/scale;
1943	>	int ny1 = ny/scale;
1944	>
1945	>	if (nx1 > floor((nx-1)/scale + 1) )
1946	>	DIE("X-dimension of output array in fsample() is too large.");
1947	>	if (ny1 > floor((ny-1)/scale + 1) )
1948	>	DIE("Y-dimension of output array in fsample() is too large.");
1949	>
1950	>	// Temp arrays
1951	>	float *bh      = (float *) (malloc(nxny * sizeof(float)));
1952	>	float *bt      = (float *) (malloc(nxny * sizeof(float)));
1953	>	float *jz      = (float *) (malloc(nxny * sizeof(float)));
1954	>	float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
1955	>	float *bpx     = (float *) (malloc(nxny * sizeof(float)));
1956	>	float *bpy     = (float *) (malloc(nxny * sizeof(float)));
1957	>	float *bpz     = (float *) (malloc(nxny * sizeof(float)));
1958	>	float *derx    = (float *) (malloc(nxny * sizeof(float)));
1959	>	float *dery    = (float *) (malloc(nxny * sizeof(float)));
1960		float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
1961		float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
1962		float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
1963		float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
1964		float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1965		float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1966	<
1967	<	// Compute
1968	<
1969	<	if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1970	<	mask, bitmask, cdelt1, rsun_ref, rsun_obs)){
1966	>	float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
1967	>	float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
1968	>	float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
1969	>	float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
1970	>	float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
1971	>	float *jz_rms_err = (float *) (malloc(nxny * sizeof(float)));
1972	>
1973	>	// malloc some arrays for the R parameter calculation
1974	>	float *rim = (float *)malloc(nx1*ny1*sizeof(float));
1975	>	float *p1p0 = (float *)malloc(nx1*ny1*sizeof(float));
1976	>	float *p1n0 = (float *)malloc(nx1*ny1*sizeof(float));
1977	>	float *p1p = (float *)malloc(nx1*ny1*sizeof(float));
1978	>	float *p1n = (float *)malloc(nx1*ny1*sizeof(float));
1979	>	float *p1 = (float *)malloc(nx1*ny1*sizeof(float));
1980	>	float *pmap = (float *)malloc(nx1*ny1*sizeof(float));
1981	>
1982	>	//spaceweather quantities computed
1983	>	if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),  &(swKeys_ptr->mean_vf_err),
1984	>	&(swKeys_ptr->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1985	>	{
1986		swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
1987		swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1988	+	swKeys_ptr->mean_vf_err = DRMS_MISSING_FLOAT;
1989	+	swKeys_ptr->count_mask  = DRMS_MISSING_INT;
1990		}
1991	<
1991	>
1992		for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
1993	<	greenpot(bpx, bpy, bpz, nx, ny);
1729	<
1730	<	computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
1993	>	greenpot(bpx, bpy, bpz, nx, ny);
1994
1995	<	if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask, bitmask))
1996	<	swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
1995	>	computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
1996	>
1997	>	if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), &(swKeys_ptr->mean_gamma_err),mask, bitmask))
1998	>	{
1999	>	swKeys_ptr->mean_gamma     =  DRMS_MISSING_FLOAT;
2000	>	swKeys_ptr->mean_gamma_err =  DRMS_MISSING_FLOAT;
2001	>	}
2002
2003	<	computeB_total(bx, by, bz, bt, dims, mask, bitmask);
2003	>	computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);
2004
2005	<	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt))
2005	>	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err)))
2006	>	{
2007		swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
2008	+	swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
2009	+	}
2010
2011	<	if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, bitmask, derx_bh, dery_bh))
2011	>	if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh), &(swKeys_ptr->mean_derivative_bh_err), mask, bitmask, derx_bh, dery_bh))
2012	>	{
2013		swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
2014	<
2015	<	if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, bitmask, derx_bz, dery_bz))
2014	>	swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT;
2015	>	}
2016	>
2017	>	if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err), mask, bitmask, derx_bz, dery_bz))
2018	>	{
2019		swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2020	+	swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT;
2021	+	}
2022
2023	<
2024	<
2025	<	if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask,
2026	<	cdelt1, rsun_ref, rsun_obs, derx, dery)) {
2027	<	swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
2028	<	swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
2023	>	computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
2024	>	derx, dery);
2025	>
2026	>
2027	>	if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &(swKeys_ptr->mean_jz),
2028	>	&(swKeys_ptr->mean_jz_err), &(swKeys_ptr->us_i), &(swKeys_ptr->us_i_err), mask, bitmask, cdelt1,
2029	>	rsun_ref, rsun_obs, derx, dery))
2030	>	{
2031	>	swKeys_ptr->mean_jz            = DRMS_MISSING_FLOAT;
2032	>	swKeys_ptr->us_i               = DRMS_MISSING_FLOAT;
2033	>	swKeys_ptr->mean_jz_err        = DRMS_MISSING_FLOAT;
2034	>	swKeys_ptr->us_i_err           = DRMS_MISSING_FLOAT;
2035		}
2036	<
2037	<	printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
2038	<
2039	<	if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2040	<	swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
2041	<
2042	<	if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih),
2043	<	&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
2044	<	mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
2045	<	swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
2046	<	swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
2047	<	swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
2036	>
2037	>	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))
2038	>	{
2039	>	swKeys_ptr->mean_alpha         = DRMS_MISSING_FLOAT;
2040	>	swKeys_ptr->mean_alpha_err     = DRMS_MISSING_FLOAT;
2041	>	}
2042	>
2043	>	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),
2044	>	&(swKeys_ptr->total_us_ih_err), &(swKeys_ptr->total_abs_ih_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2045	>	{
2046	>	swKeys_ptr->mean_ih            = DRMS_MISSING_FLOAT;
2047	>	swKeys_ptr->total_us_ih        = DRMS_MISSING_FLOAT;
2048	>	swKeys_ptr->total_abs_ih       = DRMS_MISSING_FLOAT;
2049	>	swKeys_ptr->mean_ih_err        = DRMS_MISSING_FLOAT;
2050	>	swKeys_ptr->total_us_ih_err    = DRMS_MISSING_FLOAT;
2051	>	swKeys_ptr->total_abs_ih_err   = DRMS_MISSING_FLOAT;
2052		}
2053	<
2054	<	if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),
2053	>
2054	>	if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &(swKeys_ptr->totaljz), &(swKeys_ptr->totaljz_err),
2055		mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2056	<	swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
2057	<
2058	<
1772	<	if (computeFreeEnergy(bx, by, bpx, bpy, dims,
1773	<	&(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
1774	<	mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
1775	<	swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1776	<	swKeys_ptr->totpot = DRMS_MISSING_FLOAT;
2056	>	{
2057	>	swKeys_ptr->totaljz            = DRMS_MISSING_FLOAT;
2058	>	swKeys_ptr->totaljz_err        = DRMS_MISSING_FLOAT;
2059		}
2060
2061	<	if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,
2062	<	&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),
2063	<	&(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
2061	>	if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims,
2062	>	&(swKeys_ptr->meanpot), &(swKeys_ptr->meanpot_err), &(swKeys_ptr->totpot), &(swKeys_ptr->totpot_err),
2063	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2064	>	{
2065	>	swKeys_ptr->meanpot            = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2066	>	swKeys_ptr->totpot             = DRMS_MISSING_FLOAT;
2067	>	swKeys_ptr->meanpot_err        = DRMS_MISSING_FLOAT;
2068	>	swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;
2069	>	}
2070	>
2071	>
2072	>	if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,
2073	>	&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),
2074		mask, bitmask)) {
2075	<	swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2075	>	swKeys_ptr->meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2076		swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
2077	<	swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2078	<	swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;
2077	>	swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;
2078	>	}
2079	>
2080	>	if (computeR(bz_err, los , dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,
2081	>	p1p, p1n, p1, pmap, nx1, ny1))
2082	>	{
2083	>	swKeys_ptr->Rparam = DRMS_MISSING_FLOAT;                // If fail, fill in NaN
2084		}
2085	+
2086
2087	<	// Clean up
2087	>	// Clean up the arrays
2088
2089		drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong
2090		drms_free_array(maskArray);
2091		drms_free_array(bxArray);
2092		drms_free_array(byArray);
2093		drms_free_array(bzArray);
2094	+	drms_free_array(losArray);          // Mar 7
2095	+	drms_free_array(bx_errArray);
2096	+	drms_free_array(by_errArray);
2097	+	drms_free_array(bz_errArray);
2098
2099	<	free(bh); free(bt); free(jz);
2099	>	free(bh); free(bt); free(jz); free(jz_smooth);
2100		free(bpx); free(bpy); free(bpz);
2101	<	free(derx); free(dery);
2102	<	free(derx_bt); free(dery_bt);
2103	<	free(derx_bz); free(dery_bz);
2101	>	free(derx); free(dery);
2102	>	free(derx_bt); free(dery_bt);
2103	>	free(derx_bz); free(dery_bz);
2104		free(derx_bh); free(dery_bh);
2105	<
2105	>	free(bt_err); free(bh_err);  free(jz_err);
2106	>	free(jz_err_squared); free(jz_rms_err);
2107	>	free(jz_err_squared_smooth);
2108	>
2109	>	free(rim);
2110	>	free(p1p0);
2111	>	free(p1n0);
2112	>	free(p1p);
2113	>	free(p1n);
2114	>	free(p1);
2115	>	free(pmap);
2116	>
2117		}
2118
2119	<
1807	<	/*
2119	>	/*
2120		* Set space weather indices, no error checking for now
2121		*
2122		*/
2123
2124		void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
2125		{
2126	<	drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
2127	<	drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
2128	<	drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
2129	<	drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
2130	<	drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
2131	<	drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
2132	<	drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
2133	<	drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
2134	<	drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
2135	<	drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
2136	<	drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
2137	<	drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
2138	<	drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
2139	<	drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);
2140	<	drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
2141	<	drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
2126	>	drms_setkey_float(outRec, "USFLUX",  swKeys_ptr->mean_vf);
2127	>	drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
2128	>	drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
2129	>	drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
2130	>	drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
2131	>	drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
2132	>	drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
2133	>	drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
2134	>	drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
2135	>	drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
2136	>	drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
2137	>	drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
2138	>	drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
2139	>	drms_setkey_float(outRec, "TOTPOT",  swKeys_ptr->totpot);
2140	>	drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
2141	>	drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
2142	>	drms_setkey_float(outRec, "CMASK",   swKeys_ptr->count_mask);
2143	>	drms_setkey_float(outRec, "ERRBT",   swKeys_ptr->mean_derivative_btotal_err);
2144	>	drms_setkey_float(outRec, "ERRVF",   swKeys_ptr->mean_vf_err);
2145	>	drms_setkey_float(outRec, "ERRGAM",  swKeys_ptr->mean_gamma_err);
2146	>	drms_setkey_float(outRec, "ERRBH",   swKeys_ptr->mean_derivative_bh_err);
2147	>	drms_setkey_float(outRec, "ERRBZ",   swKeys_ptr->mean_derivative_bz_err);
2148	>	drms_setkey_float(outRec, "ERRJZ",   swKeys_ptr->mean_jz_err);
2149	>	drms_setkey_float(outRec, "ERRUSI",  swKeys_ptr->us_i_err);
2150	>	drms_setkey_float(outRec, "ERRALP",  swKeys_ptr->mean_alpha_err);
2151	>	drms_setkey_float(outRec, "ERRMIH",  swKeys_ptr->mean_ih_err);
2152	>	drms_setkey_float(outRec, "ERRTUI",  swKeys_ptr->total_us_ih_err);
2153	>	drms_setkey_float(outRec, "ERRTAI",  swKeys_ptr->total_abs_ih_err);
2154	>	drms_setkey_float(outRec, "ERRJHT",  swKeys_ptr->totaljz_err);
2155	>	drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);
2156	>	drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);
2157	>	drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);
2158	>	drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);
2159		};
2160
2161	<
1833	<	/*
2161	>	/*
2162		* Set all keywords, no error checking for now
2163		*
2164		*/
2165
2166	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)
2166	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
2167		{
2168	<	copy_me_keys(inRec, outRec);
2169	<	copy_patch_keys(inRec, outRec);
2170	<	copy_geo_keys(inRec, outRec);
2171	<	copy_ambig_keys(inRec, outRec);
2172	<
2173	<	char timebuf[1024];
2174	<	float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1847	<	double val;
1848	<	int status = DRMS_SUCCESS;
2168	>
2169	>	copy_me_keys(bharpRec, outRec);
2170	>	copy_patch_keys(mharpRec, outRec);      // Dec 30
2171	>	copy_geo_keys(mharpRec, outRec);        // Dec 30
2172	>	copy_ambig_keys(bharpRec, outRec);
2173	>
2174	>	int status = 0;
2175
2176	<	val = drms_getkey_double(inRec, "DATE",&status);
2177	<	drms_setkey_double(outRec, "DATE_B", val);
2178	<	sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
2179	<	drms_setkey_string(outRec, "DATE", timebuf);
2180	<
2181	<	// set cvs commit version into keyword HEADER
2182	<	char *cvsinfo = strdup("$Header$");
2183	<	//   status = drms_setkey_string(outRec, "HEADER", cvsinfo);
2184	<	status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);
2176	>	// Change a few geometry keywords for CEA & cutout records
2177	>	if (mInfo != NULL) {        // CEA
2178	>
2179	>	drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
2180	>	drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
2181	>
2182	>	drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
2183	>	drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
2184	>	drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
2185	>	drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
2186	>	drms_setkey_string(outRec, "CUNIT1", "degree");
2187	>	drms_setkey_string(outRec, "CUNIT2", "degree");
2188	>
2189	>	char key[64];
2190	>	snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
2191	>	drms_setkey_string(outRec, "CTYPE1", key);
2192	>	snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
2193	>	drms_setkey_string(outRec, "CTYPE2", key);
2194	>	drms_setkey_float(outRec, "CROTA2", 0.0);
2195	>
2196	>	// Jan 2 2014 XS
2197	>	int nSeg = ARRLENGTH(CEASegs);
2198	>	for (int iSeg = 0; iSeg < nSeg; iSeg++) {
2199	>	DRMS_Segment_t *outSeg = NULL;
2200	>	outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);
2201	>	if (!outSeg) continue;
2202	>	// Set Bunit
2203	>	char bunit_xxx[20];
2204	>	sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
2205	>	//printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);
2206	>	drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);
2207	>	}
2208	>
2209	>	} else {        // Cutout
2210	>
2211	>	float disk_xc, disk_yc;
2212	>	if (fullDisk) {
2213	>	disk_xc = drms_getkey_float(bharpRec, "CRPIX1", &status);
2214	>	disk_yc = drms_getkey_float(bharpRec, "CRPIX2", &status);
2215	>	} else {
2216	>	disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);
2217	>	disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);
2218	>	}
2219	>	float x_ll = drms_getkey_float(mharpRec, "CRPIX1", &status);
2220	>	float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &status);
2221	>	// Defined as disk center's pixel address wrt lower-left of cutout
2222	>	drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
2223	>	drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
2224	>	// Always 0.
2225	>	drms_setkey_float(outRec, "CRVAL1", 0);
2226	>	drms_setkey_float(outRec, "CRVAL2", 0);
2227	>
2228	>	// Jan 2 2014 XS
2229	>	int nSeg = ARRLENGTH(CutSegs);
2230	>	for (int iSeg = 0; iSeg < nSeg; iSeg++) {
2231	>	DRMS_Segment_t *outSeg = NULL;
2232	>	outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);
2233	>	if (!outSeg) continue;
2234	>	// Set Bunit
2235	>	char bunit_xxx[20];
2236	>	sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
2237	>	//printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);
2238	>	drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);
2239	>	}
2240	>
2241	>
2242	>	}
2243	>
2244	>	// Mar 19 XS
2245	>	if (fullDisk) {
2246	>	drms_setkey_int(outRec, "AMBPATCH", 0);
2247	>	drms_setkey_int(outRec, "AMBWEAK", 2);
2248	>	} else {
2249	>	drms_setkey_int(outRec, "AMBPATCH", 1);
2250	>	}
2251	>
2252	>	TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
2253	>	tnow = (double)time(NULL);
2254	>	tnow += UNIX_epoch;
2255	>
2256	>	val = drms_getkey_time(bharpRec, "DATE", &status);
2257	>	drms_setkey_time(outRec, "DATE_B", val);
2258	>	drms_setkey_time(outRec, "DATE", tnow);
2259	>
2260	>	// set cvs commit version into keyword HEADER
2261	>	char *cvsinfo  = strdup("$Id$");
2262	>	char *cvsinfo2 = sw_functions_version();
2263	>	char cvsinfoall[2048];
2264	>	strcat(cvsinfoall,cvsinfo);
2265	>	strcat(cvsinfoall,"\n");
2266	>	strcat(cvsinfoall,cvsinfo2);
2267	>	status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
2268
2269		};
2270
1862	–	//
1863	–	//
2271
2272		/* ############# Nearest neighbour interpolation ############### */
2273

Diff Legend

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