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.3 by xudong, Fri Sep 7 02:14:02 2012 UTC vs.
Revision 1.35 by mbobra, Mon Mar 2 21:41:15 2015 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
2	>	*  sharp.c
3		*
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
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 33 | Line 41
41		*              SW indices fixed
42		*              Added doppler and continuum
43		*              Added other keywords: HEADER (populated by cvs build version), DATE_B
44	<	*
45	<	*      Example:
46	<	*      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
47	<	"bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
48	<	"dop=hmi.V_720s[2012.02.20_10:00]" \
49	<	"cont=hmi.Ic_720s[2012.02.20_10:00]" \
50	<	"sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
51	<	*      For comparison:
52	<	*      bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
53	<	"out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
44	>	*              v0.3
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	>	*              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		*/
76	<
76	>
77		#include <stdio.h>
78		#include <stdlib.h>
79		#include <time.h>
#	Line 64 | 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 82 | 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 93 | Line 122
122		#define SHOW(msg) {printf("%s", msg); fflush(stdout);}
123
124		#define kNotSpecified "Not Specified"
125	<
125	>
126		// Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
127		// and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
128		// 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 109 | 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;
147	<	float absFlux;
148	<	float mean_hf;
149	<	float mean_gamma;
150	<	float mean_derivative_btotal;
151	<	float mean_derivative_bh;
152	<	float mean_derivative_bz;
153	<	float mean_jz;
154	<	float us_i;
155	<	float mean_alpha;
156	<	float mean_ih;
157	<	float total_us_ih;
158	<	float total_abs_ih;
159	<	float totaljz;
160	<	float totpot;
161	<	float meanpot;
162	<	float area_w_shear_gt_45;
163	<	float meanshear_angle;
164	<	float area_w_shear_gt_45h;
165	<	float meanshear_angleh;
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;
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;
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	>	float totfx;
184	>	float totfy;
185	>	float totfz;
186	>	float totbsq;
187	>	float epsx;
188	>	float epsy;
189	>	float epsz;
190		};
191
192		// Mapping method
#	Line 149 | Line 203 | enum projection {
203		lambert
204		};
205
206	<	// Ephemeris
206	>	// WSC code
207	>	char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
208	>	"SIN", "ZEA"};
209	>
210	>	// Ephemeris information
211		struct ephemeris {
212		double disk_lonc, disk_latc;
213		double disk_xc, disk_yc;
#	Line 167 | Line 225 | struct mapInfo {
225		float *xi_out, *zeta_out;       // coordinate on full disk image to sample at
226		};
227
170	–
228		/* ========================================================================================================== */
229
230		/* Get all input data series */
#	Line 183 | Line 240 | int getInputRS_aux(DRMS_RecordSet_t **in
240		/* Find record from record set with given T_rec */
241		int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
242
186	–	// ===================
187	–
243		/* Create CEA record */
244		int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
245	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
245	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
246		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
247
248		/* Mapping single segment, wrapper */
#	Line 210 | Line 265 | int getEphemeris(DRMS_Record_t *inRec, s
265		void findCoord(struct mapInfo *mInfo);
266
267		/* Mapping function */
268	<	int performSampling(float *outData, float *inData, struct mapInfo *mInfo);
268	>	int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
269
270		/* Performing local vector transformation */
271		void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
272
273		/* Map and propogate errors */
274		int getBErr(float *bx_err, float *by_err, float *bz_err,
275	<	DRMS_Record_t *inRec, struct mapInfo *mInfo);
275	>	DRMS_Record_t *inRec, struct mapInfo *mInfo);
276
277		/* Read full disk vector magnetogram */
278		int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
279
280		/* Read variances and covariances of vector magnetograms */
281	<	int readVectorBErr(DRMS_Record_t *bharpRec,
281	>	int readVectorBErr(DRMS_Record_t *bharpRec,
282		float *bT, float *bI, float *bA,
283	<	float *errbT, float *errbI, float *errbA,
283	>	float *errbT, float *errbI, float *errbA,
284		float *errbTbI, float *errbTbA, float *errbIbA);
285
286		// ===================
#	Line 247 | Line 302 | void computeSWIndex(struct swIndex *swKe
302		void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
303
304		/* Set all keywords, no error checking for now */
305	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);
305	>	// Changed Dec 30 XS
306	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
307
308		// ===================
309
#	Line 278 | Line 334 | char *BharpSegs[] = {"inclination", "azi
334		"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
335		"disambig", "conf_disambig"};
336		// For stats
337	<	char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
337	>	char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
338		"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
339		"damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
340		"conv_flag", // fixed
#	Line 287 | Line 343 | char *CutSegs[] = {"magnetogram", "bitma
343		"field_inclination_err", "field_az_err", "inclin_azimuth_err",
344		"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
345		"disambig", "conf_disambig"};
346	<	char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",
347	<	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
348	<	/* ========================================================================================================== */
349	<
346	>	char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
347	>	BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA, "conf_disambig"};
348	>	// For Bunits, added Dec 30 XS
349	>	char *CutBunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
350	>	"degree", "degree", "Mx/cm^2", "cm/s", "mA", " ",
351	>	"length units", "DN/s", "DN/s", " ", " ",
352	>	" ",
353	>	" ", " ",
354	>	"degree", "degree", "Mx/cm^2", "cm/s", " ",
355	>	" ", " ", " ",
356	>	" ", " ", " ",
357	>	" ", " "};
358	>	char *CEABunits[] = {"Mx/cm^2", " ", "cm/s", "DN/s",
359	>	"Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", "Mx/cm^2", " "};      // Mar 4 2014 XS
360
361	+	/* ========================================================================================================== */
362
363		char *module_name = "sharp";
364	<	char *version_id = "2012 Jul 02";  /* Version number */
364	>	int seed;
365	>
366	>	int fullDisk;       // full disk mode
367	>	int amb4err;      // Use azimuth disambiguation for error propagation, default is 0 for patch and 1 for FD
368
369		ModuleArgs_t module_args[] =
370		{
371		{ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
372		{ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},
373	+	{ARG_STRING, "b", kNotSpecified, "Input B series, if set, overrides bharp."},
374		{ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},
375		{ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
376		{ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
377		{ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
378	+	{ARG_INT,    "seed", "987654", "Seed for the random number generator."},
379	+	{ARG_INT,   "f_amb4err", "0", "Force using disambiguation in error propagation"},     // Mar 4 2014 XS
380		{ARG_END}
381		};
382
383	<	int DoIt(void)
383	>	int DoIt(void)
384		{
385	<
385	>	int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
386	>	int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
387	>
388		int status = DRMS_SUCCESS;
389		int nrecs, irec;
390
391	<	char *mharpQuery, *bharpQuery;
391	>	char *mharpQuery, *bharpQuery, *bQuery;
392		char *dopQuery, *contQuery;
393		char *sharpCeaQuery, *sharpCutQuery;
394
#	Line 324 | Line 399 | int DoIt(void)
399
400		mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
401		bharpQuery = (char *) params_get_str(&cmdparams, "bharp");
402	+	bQuery = (char *) params_get_str(&cmdparams, "b");
403		dopQuery = (char *) params_get_str(&cmdparams, "dop");
404		contQuery = (char *) params_get_str(&cmdparams, "cont");
405		sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
406		sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
407	+
408	+	seed = params_get_int(&cmdparams, "seed");
409	+	int f_amb4err = params_get_int(&cmdparams, "f_amb4err");
410
411		/* Get input data, check everything */
412
413	<	if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
414	<	DIE("Input harp data error.");
413	>	// Full disk mode if "b" is set
414	>	if (strcmp(bQuery, kNotSpecified)) {
415	>	fullDisk = 1;
416	>	bharpQuery = bQuery;
417	>	//        SHOW(bharpQuery); SHOW("\n");
418	>	SHOW("Full disk mode\n");
419	>	} else {
420	>	fullDisk = 0;
421	>	SHOW("Harp mode\n");
422	>	}
423	>
424	>	// Mar 4 2014
425	>	if (f_amb4err == 0) {         // no forcing, 0 for patch and 1 for FD
426	>	amb4err = fullDisk ? 1 : 0;
427	>	} else {
428	>	amb4err = 1;
429	>	}
430	>	printf("amb4err=%d\n", amb4err);
431	>
432	>	// Bharp point to B if full disk
433	>	if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
434	>	DIE("Input harp data error.");
435		nrecs = mharpRS->n;
436
437	<	if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
437	>	if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
438		DIE("Input doppler data error.");
439
440	<	if (getInputRS_aux(&contRS, contQuery, mharpRS))
440	>	if (getInputRS_aux(&contRS, contQuery, mharpRS))
441		DIE("Input continuum data error.");
442
443		/* Start */
444
445		printf("==============\nStart. %d image(s) in total.\n", nrecs);
446	<
446	>
447		for (irec = 0; irec < nrecs; irec++) {
448
449		/* Records in work */
450
451		DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;
452	+
453		mharpRec = mharpRS->records[irec];
454	<	bharpRec = bharpRS->records[irec];
455	<
456	<	TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
457	<
454	>
455	>	TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
456	>
457	>	if (!fullDisk) {
458	>	bharpRec = bharpRS->records[irec];
459	>	} else {
460	>	if (getInputRec_aux(&bharpRec, bharpRS, trec)) {     // Bharp point to full disk B
461	>	printf("Fetching B failed, image #%d skipped.\n", irec);
462	>	continue;
463	>	}
464	>	}
465	>
466		struct swIndex swKeys;
467
468		DRMS_Record_t *dopRec = NULL, *contRec = NULL;
469	+
470		if (getInputRec_aux(&dopRec, dopRS, trec)) {
471		printf("Fetching Doppler failed, image #%d skipped.\n", irec);
472		continue;
#	Line 366 | Line 475 | int DoIt(void)
475		printf("Fetching continuum failed, image #%d skipped.\n", irec);
476		continue;
477		}
478	<
478	>
479		/* Create CEA record */
480
481		DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
#	Line 404 | Line 513 | int DoIt(void)
513		printf("Image #%d done.\n", irec);
514
515		} // irec
516	+
517
518		drms_close_records(mharpRS, DRMS_FREE_RECORD);
519		drms_close_records(bharpRS, DRMS_FREE_RECORD);
520	+	drms_close_records(dopRS, DRMS_FREE_RECORD);                            // Dec 18 2012
521	+	drms_close_records(contRS, DRMS_FREE_RECORD);                           // Dec 18 2012
522
523		return 0;
524	<
524	>
525		}       // DoIt
526
527
#	Line 433 | Line 545 | int getInputRS(DRMS_RecordSet_t **mharpR
545		*mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
546		if (status || (*mharpRS_ptr)->n == 0) return 1;
547
548	<	*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
549	<	if (status || (*bharpRS_ptr)->n == 0) return 1;
550	<
551	<	if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
548	>	if (fullDisk) {
549	>	if (getInputRS_aux(bharpRS_ptr, bharpQuery, *mharpRS_ptr)) return 1;
550	>	} else {
551	>	*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
552	>	if (status || (*bharpRS_ptr)->n == 0) return 1;
553	>	if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
554	>	}
555
556		return 0;
557
#	Line 462 | Line 577 | int compareHarp(DRMS_RecordSet_t *mharpR
577		for (int i = 0; i < nrecs; i++) {
578		mharpRec_t = mharpRS->records[i];
579		bharpRec_t = bharpRS->records[i];
580	<	if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
580	>	if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
581		drms_getkey_int(bharpRec_t, "HARPNUM", &status)) ||
582	<	(drms_getkey_time(mharpRec_t, "T_REC", &status) !=
582	>	(drms_getkey_time(mharpRec_t, "T_REC", &status) !=
583		drms_getkey_time(bharpRec_t, "T_REC", &status)))
584		{
585		return 1;
#	Line 475 | Line 590 | int compareHarp(DRMS_RecordSet_t *mharpR
590
591		}
592
593	<	/*
593	>	/*
594		* Get other data series, check all T_REC are available
595		*
596		*/
#	Line 490 | Line 605 | int getInputRS_aux(DRMS_RecordSet_t **in
605
606		// Check if all T_rec are available, need to match both ways
607		int n = harpRS->n, n0 = (*inRS_ptr)->n;
608	<
608	>
609		for (int i0 = 0; i0 < n0; i0++) {
610		DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
611		TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
#	Line 519 | Line 634 | int getInputRS_aux(DRMS_RecordSet_t **in
634
635		}
636
637	<	/*
637	>	/*
638		* Find record from record set with given T_rec
639		*
640		*/
#	Line 549 | Line 664 | int getInputRec_aux(DRMS_Record_t **inRe
664		*
665		*/
666
667	<	int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
668	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
667	>	int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
668	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
669		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
670		{
671
#	Line 562 | Line 677 | int createCeaRecord(DRMS_Record_t *mharp
677		mInfo.proj = (enum projection) cyleqa;          // projection method
678		mInfo.xscale = XSCALE;
679		mInfo.yscale = YSCALE;
680	<	mInfo.nbin = NBIN;
680	>
681	>	int ncol0, nrow0;           // oversampled map size
682
683		// Get ephemeris
684	<
684	>
685		if (getEphemeris(mharpRec, &(mInfo.ephem))) {
686		SHOW("CEA: get ephemeris error\n");
687		return 1;
#	Line 578 | Line 694 | int createCeaRecord(DRMS_Record_t *mharp
694		return 1;
695		}
696
697	+	// ========================================
698	+	// Do this for all bitmaps, Aug 12 2013 XS
699	+	// ========================================
700	+
701	+	mInfo.nbin = 1;                     // for bitmaps. suppress anti-aliasing
702	+	ncol0 = mInfo.ncol;
703	+	nrow0 = mInfo.nrow;
704	+
705	+	mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
706	+	mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
707	+
708	+	findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
709	+
710	+	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
711	+	SHOW("CEA: mapping bitmap error\n");
712	+	return 1;
713	+	}
714	+	printf("Bitmap mapping done.\n");
715	+
716	+	if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
717	+	SHOW("CEA: mapping conf_disambig error\n");
718	+	return 1;
719	+	}
720	+	printf("Conf disambig mapping done.\n");
721	+
722	+	free(mInfo.xi_out);
723	+	free(mInfo.zeta_out);
724	+
725	+	// ========================================
726	+	// Do this again for floats, Aug 12 2013 XS
727	+	// ========================================
728		// Create xi_out, zeta_out array in mInfo:
729		// Coordinates to sample in original full disk image
730
731	<	int ncol0, nrow0;               // oversampled map size
731	>	mInfo.nbin = NBIN;
732		ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
733		nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
734
#	Line 591 | Line 738 | int createCeaRecord(DRMS_Record_t *mharp
738		findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
739
740		// Mapping single segment: Mharp, etc.
741	<
741	>
742		if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
743		SHOW("CEA: mapping magnetogram error\n");
744		return 1;
745		}
599	–	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
600	–	SHOW("CEA: mapping bitmap error\n");
601	–	return 1;
602	–	}
746		printf("Magnetogram mapping done.\n");
747	<
747	>
748		if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
749		SHOW("CEA: mapping dopplergram error\n");
750		return 1;
#	Line 613 | Line 756 | int createCeaRecord(DRMS_Record_t *mharp
756		return 1;
757		}
758		printf("Intensitygram mapping done.\n");
759	<
617	<	if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
618	<	SHOW("CEA: mapping conf_disambig error\n");
619	<	return 1;
620	<	}
621	<	printf("Conf disambig mapping done.\n");
622	<
759	>
760		// Mapping vector B
761
762		if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
#	Line 641 | Line 778 | int createCeaRecord(DRMS_Record_t *mharp
778		drms_copykey(sharpRec, mharpRec, "T_REC");
779		drms_copykey(sharpRec, mharpRec, "HARPNUM");
780
781	+	if (fullDisk) {
782	+	DRMS_Link_t *bLink = hcon_lookup_lower(&sharpRec->links, "B");
783	+	if (bLink) drms_link_set("B", sharpRec, bharpRec);
784	+	} else {
785	+	DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
786	+	if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
787	+	}
788		DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
789		if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
646	–	DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
647	–	if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
790
791	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
791	>	setKeys(sharpRec, mharpRec, bharpRec, &mInfo);            // Set all other keywords
792		drms_copykey(sharpRec, mharpRec, "QUALITY");            // copied from los records
793	<
793	>
794		// Space weather
795
796		computeSWIndex(swKeys_ptr, sharpRec, &mInfo);           // compute it!
797		printf("Space weather indices done.\n");
798
799		setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
800	<
800	>
801		// Stats
802
803		int nCEASegs = ARRLENGTH(CEASegs);
#	Line 663 | Line 805 | int createCeaRecord(DRMS_Record_t *mharp
805		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
806		DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
807		int stat = set_statistics(outSeg, outArray, 1);
808	<	//              printf("%d => %d\n", iSeg, stat);
808	>	//              printf("%d => %d\n", iSeg, stat);
809		drms_free_array(outArray);
810		}
811
#	Line 674 | Line 816 | int createCeaRecord(DRMS_Record_t *mharp
816		}
817
818
819	<	/*
819	>	/*
820		* Mapping a single segment
821		* Read in full disk image, utilize mapImage for mapping
822		* then write the segment out, segName same in in/out Rec
#	Line 688 | Line 830 | int mapScaler(DRMS_Record_t *sharpRec, D
830		int status = 0;
831		int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
832		int dims[2] = {nx, ny};
833	+	int interpOpt = INTERP;         // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
834
835		// Input full disk array
836
#	Line 699 | Line 842 | int mapScaler(DRMS_Record_t *sharpRec, D
842		inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
843		if (!inArray) return 1;
844
845	+	if (!strcmp(segName, "conf_disambig") || !strcmp(segName, "bitmap")) {
846	+	// Moved out so it works for FD conf_disambig as well
847	+	// Jan 2 2014 XS
848	+	interpOpt = 3;          // Aug 12 XS, near neighbor
849	+	}
850	+
851		float *inData;
852		int xsz = inArray->axis[0], ysz = inArray->axis[1];
853		if ((xsz != FOURK) || (ysz != FOURK)) {         // for bitmap, make tmp full disk
#	Line 721 | Line 870 | int mapScaler(DRMS_Record_t *sharpRec, D
870		// Mapping
871
872		float *map = (float *) (malloc(nxny * sizeof(float)));
873	<	if (performSampling(map, inData, mInfo))
874	<	{if (inArray) drms_free_array(inArray); free(map); return 1;}
873	>	if (performSampling(map, inData, mInfo, interpOpt))             // Add interpOpt for different types, Aug 12 XS
874	>	{if (inArray) drms_free_array(inArray); free(map); return 1;}
875
876		// Write out
877
878	<	DRMS_Segment_t *outSeg = NULL;
878	>	DRMS_Segment_t *outSeg = NULL;
879		outSeg = drms_segment_lookup(sharpRec, segName);
880		if (!outSeg) return 1;
881
882	<	DRMS_Type_t arrayType = outSeg->info->type;
882	>	//  DRMS_Type_t arrayType = outSeg->info->type;
883		DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
884		if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
885
886		// convert to needed data type
887
888	<	drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);
888	>	//  drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);         // Jan 02 2013
889
890		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
891	<	outArray->parent_segment = outSeg;
891	>	//  outArray->parent_segment = outSeg;
892		outArray->israw = 0;            // always compressed
893		outArray->bzero = outSeg->bzero;
894		outArray->bscale = outSeg->bscale;
895	<
895	>
896		status = drms_segment_write(outSeg, outArray, 0);
897		if (status) return 0;
898
899		if (inArray) drms_free_array(inArray);
900	+	if ((xsz != FOURK) || (ysz != FOURK)) free(inData);                     // Dec 18 2012
901		if (outArray) drms_free_array(outArray);
902		return 0;
903
904		}
905
906
907	<	/*
907	>	/*
908		* Mapping vector magnetogram
909		*
910		*/
#	Line 781 | Line 931 | int mapVectorB(DRMS_Record_t *sharpRec,
931		// Mapping
932
933		float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation
934	<
934	>
935		bx_map = (float *) (malloc(nxny * sizeof(float)));
936	<	if (performSampling(bx_map, bx_img, mInfo))
937	<	{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
938	<
936	>	if (performSampling(bx_map, bx_img, mInfo, INTERP))
937	>	{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
938	>
939		by_map = (float *) (malloc(nxny * sizeof(float)));
940	<	if (performSampling(by_map, by_img, mInfo))
941	<	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
942	<
940	>	if (performSampling(by_map, by_img, mInfo, INTERP))
941	>	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
942	>
943		bz_map = (float *) (malloc(nxny * sizeof(float)));
944	<	if (performSampling(bz_map, bz_img, mInfo))
945	<	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
944	>	if (performSampling(bz_map, bz_img, mInfo, INTERP))
945	>	{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
946
947		free(bx_img); free(by_img); free(bz_img);
948
#	Line 814 | Line 964 | int mapVectorB(DRMS_Record_t *sharpRec,
964		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
965		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
966		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
967	<	outArray->parent_segment = outSeg;
967	>	//              outArray->parent_segment = outSeg;
968		outArray->israw = 0;
969		outArray->bzero = outSeg->bzero;
970		outArray->bscale = outSeg->bscale;
#	Line 830 | Line 980 | int mapVectorB(DRMS_Record_t *sharpRec,
980		}
981
982
983	<	/*
983	>	/*
984		* Mapping vector magnetogram errors
985		*
986		*/
#	Line 853 | Line 1003 | int mapVectorBErr(DRMS_Record_t *sharpRe
1003		free(bx_err); free(by_err); free(bz_err);
1004		return 1;
1005		}
1006	<
1006	>
1007		// Write out
1008
1009		DRMS_Segment_t *outSeg;
#	Line 866 | Line 1016 | int mapVectorBErr(DRMS_Record_t *sharpRe
1016		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
1017		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
1018		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
1019	<	outArray->parent_segment = outSeg;
1019	>	//              outArray->parent_segment = outSeg;
1020		outArray->israw = 0;
1021		outArray->bzero = outSeg->bzero;
1022		outArray->bscale = outSeg->bscale;
#	Line 883 | Line 1033 | int mapVectorBErr(DRMS_Record_t *sharpRe
1033
1034
1035
1036	<	/*
1036	>	/*
1037		* Determine reference point coordinate and patch size according to keywords
1038		* xc, yc are the coordinate of patch center, in degrees
1039		* ncol and nrow are the final size
#	Line 899 | Line 1049 | int findPosition(DRMS_Record_t *inRec, s
1049		float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
1050
1051		/* Center coord */
1052	+	// Changed into double Jun 16 2014 XS
1053
1054	<	float minlon = drms_getkey_float(inRec, "LONDTMIN", &status); if (status) return 1;             // Stonyhurst lon
1055	<	float maxlon = drms_getkey_float(inRec, "LONDTMAX", &status); if (status) return 1;
1056	<	float minlat = drms_getkey_float(inRec, "LATDTMIN", &status); if (status) return 1;
1057	<	float maxlat = drms_getkey_float(inRec, "LATDTMAX", &status); if (status) return 1;
1054	>	double minlon = drms_getkey_double(inRec, "LONDTMIN", &status); if (status) return 1;           // Stonyhurst lon
1055	>	double maxlon = drms_getkey_double(inRec, "LONDTMAX", &status); if (status) return 1;
1056	>	double minlat = drms_getkey_double(inRec, "LATDTMIN", &status); if (status) return 1;
1057	>	double maxlat = drms_getkey_double(inRec, "LATDTMAX", &status); if (status) return 1;
1058
1059		// A bug fixer for HARP (per M. Turmon)
1060		// When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
#	Line 911 | Line 1062 | int findPosition(DRMS_Record_t *inRec, s
1062		// We compute minlon & minlat then by
1063		// LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
1064
1065	<	float psize = drms_getkey_float(inRec, "SIZE", &status);
1066	<	if (psize != psize) {
1067	<	TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1;
1065	>	//  float psize = drms_getkey_float(inRec, "SIZE", &status);
1066	>	//  if (psize != psize) {
1067	>
1068	>	if (minlon != minlon || maxlon != maxlon) {         // check lons instead of SIZE
1069	>	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
1070		double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
1071		char firstRecQuery[100], t0_str[100];
1072		sprint_time(t0_str, t0, "TAI", 0);
#	Line 932 | Line 1085 | int findPosition(DRMS_Record_t *inRec, s
1085		mInfo->yc = (maxlat + minlat) / 2.;
1086
1087		/* Size */
1088	+	// Rounded to 1.d3 precision first. Jun 16 2014 XS
1089
1090	<	mInfo->ncol = round((maxlon - minlon) / mInfo->xscale);
1091	<	mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
1090	>	mInfo->ncol = round(round((maxlon - minlon) * 1.e3) / 1.e3 / mInfo->xscale);
1091	>	mInfo->nrow = round(round((maxlat - minlat) * 1.e3) / 1.e3 / mInfo->yscale);
1092
1093		return 0;
1094
#	Line 953 | Line 1107 | int getEphemeris(DRMS_Record_t *inRec, s
1107		int status = 0;
1108
1109		float crota2 = drms_getkey_float(inRec, "CROTA2", &status);     // rotation
1110	<	double sina = sin(crota2 * RADSINDEG);
1110	>	double sina = sin(crota2 * RADSINDEG);
1111		double cosa = cos(crota2 * RADSINDEG);
1112
1113		ephem->pa = - crota2 * RADSINDEG;
#	Line 965 | Line 1119 | int getEphemeris(DRMS_Record_t *inRec, s
1119		float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
1120		float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
1121		float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
968	–	printf("cdelt=%f\n",cdelt);
1122		ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0
1123		ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
1124
#	Line 1030 | Line 1183 | void findCoord(struct mapInfo *mInfo)
1183		x = (col0 + 0.5 - ncol0/2.) * xscale0;          // in rad
1184		y = (row0 + 0.5 - nrow0/2.) * yscale0;
1185
1186	<	/* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1186	>	/* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1187		* the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
1188		* projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
1189	<	* is the heliographic longitude and latitude of the map center. Both are in degree.
1189	>	* is the heliographic longitude and latitude of the map center. Both are in degree.
1190		*/
1191
1192		if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
#	Line 1046 | Line 1199 | void findCoord(struct mapInfo *mInfo)
1199		* image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
1200		*/
1201
1202	<	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1202	>	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1203		disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1204		xi_out[ind_map] = -1;
1205		zeta_out[ind_map] = -1;
#	Line 1062 | Line 1215 | void findCoord(struct mapInfo *mInfo)
1215		}
1216
1217
1218	<	/*
1218	>	/*
1219		* Sampling function
1220		* oversampling by nbin, then binning using a Gaussian
1221		* save results in outData, always of float type
1222		*
1223		*/
1224
1225	<	int performSampling(float *outData, float *inData, struct mapInfo *mInfo)
1225	>	int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
1226		{
1227
1228		int status = 0;
1229	+	int ind_map;
1230
1231		int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN
1232		int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
1233
1234	<	float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1234	>	// Changed Aug 12 2013, XS, for bitmaps
1235	>	float *outData0;
1236	>	if (interpOpt == 3 && mInfo->nbin == 1) {
1237	>	outData0 = outData;
1238	>	} else {
1239	>	outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1240	>	}
1241
1242		float *xi_out = mInfo->xi_out;
1243		float *zeta_out = mInfo->zeta_out;
1244	<
1244	>
1245		// Interpolation
1246
1247		struct fint_struct pars;
1248	<	int interpOpt = INTERP;         // Use Wiener by default, 6 order, 1 constraint
1248	>	// Aug 12 2013, passed in as argument now
1249
1250		switch (interpOpt) {
1251		case 0:                 // Wiener, 6 order, 1 constraint
#	Line 1097 | Line 1257 | int performSampling(float *outData, floa
1257		case 2:                 // Bilinear
1258		init_finterpolate_linear(&pars, 1.);
1259		break;
1260	+	case 3:                 // Near neighbor
1261	+	break;
1262		default:
1263		return 1;
1264		}
1265
1266	<	finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1267	<	FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1266	>	printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
1267	>	if (interpOpt == 3) {                   // Aug 6 2013, Xudong
1268	>	for (int row0 = 0; row0 < nrow0; row0++) {
1269	>	for (int col0 = 0; col0 < ncol0; col0++) {
1270	>	ind_map = row0 * ncol0 + col0;
1271	>	outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
1272	>	}
1273	>	}
1274	>	} else {
1275	>	finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1276	>	FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1277	>	}
1278
1279		// Rebinning, smoothing
1280
1281	<	frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
1281	>	if (interpOpt == 3 && mInfo->nbin == 1) {
1282	>	return 0;
1283	>	} else {
1284	>	frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
1285	>	free(outData0);         // Dec 18 2012
1286	>	}
1287
1288		//
1289
#	Line 1115 | Line 1292 | int performSampling(float *outData, floa
1292		}
1293
1294
1295	<	/*
1295	>	/*
1296		* Performing local vector transformation
1297		*  xyz: z refers to vertical (radial) component, x EW (phi), y NS
1298		*
#	Line 1176 | Line 1353 | void vectorTransform(float *bx_map, floa
1353
1354		}
1355		}
1356	<
1356	>
1357		}
1358
1359
1360
1361	<	/*
1361	>	/*
1362		* Map and propogate vector field errors
1363		*
1364		*/
1365
1366		int getBErr(float *bx_err, float *by_err, float *bz_err,
1367	<	DRMS_Record_t *inRec, struct mapInfo *mInfo)
1367	>	DRMS_Record_t *inRec, struct mapInfo *mInfo)
1368		{
1369
1370		int status = 0;
#	Line 1206 | Line 1383 | int getBErr(float *bx_err, float *by_err
1383		float *errbTbA = (float *) (malloc(FOURK2 * sizeof(float)));
1384		float *errbIbA = (float *) (malloc(FOURK2 * sizeof(float)));
1385
1386	<	if (readVectorBErr(inRec,
1386	>	if (readVectorBErr(inRec,
1387		bT, bI, bA,
1388	<	errbT, errbI, errbA,
1388	>	errbT, errbI, errbA,
1389		errbTbI, errbTbA, errbIbA)) {
1390		printf("Read full disk variances & covariances error\n");
1391		free(bT); free(bI); free(bA);
#	Line 1245 | Line 1422 | int getBErr(float *bx_err, float *by_err
1422		int ind_map, ind_img;
1423
1424		double bpSigma2, btSigma2, brSigma2;            // variances after prop
1425	<
1425	>
1426		for (int row = 0; row < mInfo->nrow; row++) {
1427		for (int col = 0; col < mInfo->ncol; col++) {
1428
#	Line 1261 | Line 1438 | int getBErr(float *bx_err, float *by_err
1438		continue;
1439		}
1440
1441	<	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1441	>	if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1442		disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1443		bx_err[ind_map] = DRMS_MISSING_FLOAT;
1444	<	bx_err[ind_map] = DRMS_MISSING_FLOAT;
1445	<	bx_err[ind_map] = DRMS_MISSING_FLOAT;
1444	>	by_err[ind_map] = DRMS_MISSING_FLOAT;
1445	>	bz_err[ind_map] = DRMS_MISSING_FLOAT;       // Mar 7
1446		continue;
1447		}
1448
#	Line 1274 | Line 1451 | int getBErr(float *bx_err, float *by_err
1451
1452		ind_img = round(zeta * FOURK + xi);
1453
1454	<	if (errorprop(bT, bA, bI,
1455	<	errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1456	<	lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1454	>	if (errorprop(bT, bA, bI,
1455	>	errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1456	>	lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1457		&btSigma2, &bpSigma2, &brSigma2)) {
1458		bx_err[ind_map] = DRMS_MISSING_FLOAT;
1459		by_err[ind_map] = DRMS_MISSING_FLOAT;
#	Line 1315 | Line 1492 | int readVectorB(DRMS_Record_t *inRec, fl
1492
1493		DRMS_Segment_t *inSeg;
1494		DRMS_Array_t *inArray_ambig;
1495	<	DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1495	>	DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1496
1497		char *ambig;
1498		float *bTotal, *bAzim, *bIncl;
#	Line 1344 | Line 1521 | int readVectorB(DRMS_Record_t *inRec, fl
1521
1522		int llx, lly;           // lower-left corner
1523		int bmx, bmy;           // bitmap size
1524	<
1525	<	llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1526	<	lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1527	<
1528	<	bmx = inArray_ambig->axis[0];
1529	<	bmy = inArray_ambig->axis[1];
1524	>
1525	>	if (fullDisk) {
1526	>	llx = lly = 0;
1527	>	bmx = bmy = FOURK;
1528	>	} else {
1529	>	llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1530	>	lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1531	>	bmx = inArray_ambig->axis[0];
1532	>	bmy = inArray_ambig->axis[1];
1533	>	}
1534
1535		int kx, ky, kOff;
1536		int ix = 0, jy = 0, yOff = 0, iData = 0;
1537		int xDim = FOURK, yDim = FOURK;
1538	+	int amb = 0;
1539
1540		for (jy = 0; jy < yDim; jy++)
1541		{
#	Line 1376 | Line 1558 | int readVectorB(DRMS_Record_t *inRec, fl
1558		continue;
1559		} else {
1560		kOff = ky * bmx + kx;
1561	<	if (ambig[kOff] % 2) {          // 180
1561	>	//                              if (ambig[kOff] % 2) {          // 180
1562	>	// Feb 12 2014, use bit #2 for full disk, lowest bit for patch
1563	>	if (fullDisk) { amb = (ambig[kOff] / 4) % 2; } else { amb = ambig[kOff] % 2; }
1564	>	if (amb) {                              // Feb 12 2014, use bit #2
1565		bx_img[iData] *= -1.; by_img[iData] *= -1.;
1566	<	}
1566	>	}
1567		}
1568		}
1569		}
#	Line 1400 | Line 1585 | int readVectorB(DRMS_Record_t *inRec, fl
1585		*
1586		*/
1587
1588	<	int readVectorBErr(DRMS_Record_t *inRec,
1588	>	int readVectorBErr(DRMS_Record_t *inRec,
1589		float *bT, float *bI, float *bA,
1590	<	float *errbT, float *errbI, float *errbA,
1590	>	float *errbT, float *errbI, float *errbA,
1591		float *errbTbI, float *errbTbA, float *errbIbA)
1592		{
1593
#	Line 1410 | Line 1595 | int readVectorBErr(DRMS_Record_t *inRec,
1595
1596		float *data_ptr[9];
1597		char *segName[9] = {"field", "inclination", "azimuth",
1598	<	"field_err", "inclination_err", "azimuth_err",
1599	<	"field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1598	>	"field_err", "inclination_err", "azimuth_err",
1599	>	"field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1600		DRMS_Segment_t *inSegs[9];
1601		DRMS_Array_t *inArrays[9];
1602
1603		// Read full disk images
1604	+	// Do we need disambig? Dec 30 XS
1605
1606		for (int iSeg = 0; iSeg < 9; iSeg++) {
1607
#	Line 1429 | Line 1615 | int readVectorBErr(DRMS_Record_t *inRec,
1615		float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];
1616		float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];
1617
1618	+	// Add disambig, Feb 12 2014
1619	+
1620	+	DRMS_Segment_t *inSeg;
1621	+	DRMS_Array_t *inArray_ambig;
1622	+
1623	+	if (amb4err) {              // Mar 4 2014
1624	+
1625	+	inSeg = drms_segment_lookup(inRec, "disambig");
1626	+	inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
1627	+	if (status) return 1;
1628	+	char *ambig = (char *)inArray_ambig->data;
1629	+
1630	+	int llx, lly;           // lower-left corner
1631	+	int bmx, bmy;           // bitmap size
1632	+
1633	+	if (fullDisk) {
1634	+	llx = lly = 0;
1635	+	bmx = bmy = FOURK;
1636	+	} else {
1637	+	llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1638	+	lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1639	+	bmx = inArray_ambig->axis[0];
1640	+	bmy = inArray_ambig->axis[1];
1641	+	}
1642	+
1643	+	int idx, idx_a;
1644	+	int amb;
1645	+
1646	+	for (int j = 0; j < bmy; j++) {
1647	+	for (int i = 0; i < bmx; i++) {
1648	+	idx_a = j * bmx + i;
1649	+	idx = (j + lly) * FOURK + (i + llx);
1650	+	// Feb 12 2014, use bit #2 for full disk, lowest bit for patch
1651	+	if (fullDisk) { amb = (ambig[idx_a] / 4) % 2; } else { amb = ambig[idx_a] % 2; }
1652	+	if (amb) { bA0[idx] += 180.; }
1653	+	}
1654	+	}
1655	+
1656	+	}
1657	+
1658		// Convert errors to variances, correlation coefficients to covariances
1659
1660		for (int i = 0; i < FOURK2; i++) {
#	Line 1437 | Line 1663 | int readVectorBErr(DRMS_Record_t *inRec,
1663		if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;
1664
1665		bT[i] = bT0[i];
1666	<	bI[i] = bI0[i];
1666	>	bI[i] = bI0[i];         // in deg, coverted in errorprop
1667		bA[i] = bA0[i];
1668
1669		errbT[i] = errbT0[i] * errbT0[i];
#	Line 1447 | Line 1673 | int readVectorBErr(DRMS_Record_t *inRec,
1673		errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG;
1674		errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG;
1675		errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1676	<
1676	>
1677		}
1678
1679		//
1680
1681		for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
1682	<
1682	>	if (amb4err) drms_free_array(inArray_ambig);            // Feb 12; Mar 04 2014
1683	>
1684		return 0;
1685
1686		}
#	Line 1463 | Line 1690 | int readVectorBErr(DRMS_Record_t *inRec,
1690		* Create Cutout record: top level subroutine
1691		* Do the loops on segments and set the keywords here
1692		* Work is done in writeCutout routine below
1693	<	*
1693	>	*
1694		*/
1695
1696	<	int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1697	<	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1696	>	int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1697	>	DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1698		DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
1699		{
1700
#	Line 1528 | Line 1755 | int createCutRecord(DRMS_Record_t *mharp
1755		if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
1756
1757		setSWIndex(sharpRec, swKeys_ptr);       // Set space weather indices
1758	<	setKeys(sharpRec, bharpRec);            // Set all other keywords
1759	<
1758	>	setKeys(sharpRec, mharpRec, bharpRec, NULL);              // Set all other keywords, NULL specifies cutout
1759	>
1760		// Stats
1761	<
1761	>
1762		int nCutSegs = ARRLENGTH(CutSegs);
1763		for (int iSeg = 0; iSeg < nCutSegs; iSeg++) {
1764		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
#	Line 1539 | Line 1766 | int createCutRecord(DRMS_Record_t *mharp
1766		set_statistics(outSeg, outArray, 1);
1767		drms_free_array(outArray);
1768		}
1769	<
1769	>
1770		return 0;
1771
1772	<	}
1772	>	}
1773
1774
1775	<	/*
1775	>	/*
1776		* Get cutout and write segment
1777		* Change DISAMB_AZI to apply disambiguation to azimuth
1778		*
#	Line 1584 | Line 1811 | int writeCutout(DRMS_Record_t *outRec, D
1811		} else {
1812		return 1;
1813		}
1814	<
1814	>
1815	>	// Feb 12 2014, fool-proof, for patch, change everything to 0 or 7!!!
1816	>	// This is a fix for disambiguation before Aug 2013
1817	>
1818	>	if (!strcmp(SegName, "disambig") && !fullDisk) {
1819	>	double *disamb = (double *) (cutoutArray->data);
1820	>	for (int i = 0; i < nxny; i++) {
1821	>	if (((int)disamb[i]) % 2) { disamb[i] = 7; } else { disamb[i] = 0; }
1822	>	}
1823	>	}
1824	>
1825		/* Adding disambiguation resolution to cutout azimuth? */
1826	<
1826	>
1827		#if DISAMB_AZI
1828	+	int amb;
1829		if (!strcmp(SegName, "azimuth")) {
1830	<	DRMS_Segment_t *disambSeg = drms_segment_lookup(inRec, "disambig");
1830	>	DRMS_Segment_t *disambSeg = NULL;
1831	>	disambSeg = drms_segment_lookup(inRec, "disambig");
1832		if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
1833		DRMS_Array_t *disambArray;
1834	<	if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1835	<	disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1836	<	if (status) {drms_free_array(cutoutArray); return 1;}
1837	<	} else {
1838	<	return 1;
1839	<	}
1834	>	if (fullDisk) { // Jan 2 2014 XS
1835	>	disambArray = drms_segment_readslice(disambSeg, DRMS_TYPE_CHAR, ll, ur, &status);
1836	>	if (status) return 1;
1837	>	} else {
1838	>	if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1839	>	disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1840	>	if (status) {drms_free_array(cutoutArray); return 1;}
1841	>	} else {
1842	>	drms_free_array(cutoutArray);
1843	>	return 1;
1844	>	}
1845	>	}
1846		double *azimuth = (double *) cutoutArray->data;
1847		char *disamb = (char *) disambArray->data;
1848		for (int n = 0; n < nxny; n++) {
1849	<	if (disamb[n]) azimuth[n] += 180.;
1849	>	//                  if (disamb[n] % 2) azimuth[n] += 180.;      // Nov 12 2013 Fixed!!!
1850	>	// Feb 12 2014, use bit #2 for full disk, lowest bit for patch
1851	>	if (fullDisk) { amb = (disamb[n] / 4) % 2; } else { amb = disamb[n] % 2; }
1852	>	if (amb) azimuth[n] += 180.;
1853		}
1854		drms_free_array(disambArray);
1855		}
1856		#endif
1857	<
1857	>
1858		/* Write out */
1859
1860		outSeg = drms_segment_lookup(outRec, SegName);
1861		if (!outSeg) return 1;
1862	<	drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);
1862	>	//  drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);      // Jan 02 2013
1863		outSeg->axis[0] = cutoutArray->axis[0];
1864		outSeg->axis[1] = cutoutArray->axis[1];
1865	<	cutoutArray->parent_segment = outSeg;
1865	>	//  cutoutArray->parent_segment = outSeg;
1866		cutoutArray->israw = 0;         // always compressed
1867		cutoutArray->bzero = outSeg->bzero;
1868		cutoutArray->bscale = outSeg->bscale;               // Same as inArray's
#	Line 1627 | Line 1875 | int writeCutout(DRMS_Record_t *outRec, D
1875		}
1876
1877
1878	<	/*
1878	>	/*
1879		* Compute space weather indices, no error checking for now
1880		* Based on M. Bobra's swharp_vectorB.c
1881		* No error checking for now
#	Line 1641 | Line 1889 | void computeSWIndex(struct swIndex *swKe
1889		int nx = mInfo->ncol, ny = mInfo->nrow;
1890		int nxny = nx * ny;
1891		int dims[2] = {nx, ny};
1892	+
1893		// Get bx, by, bz, mask
1894
1895	<	// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1896	<	//DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "bitmap");
1897	<	//DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1898	<	//int *mask = (int *) maskArray->data;          // get the previously made mask array
1899	<
1900	<	//Use conf_disambig map as a threshold on spaceweather quantities
1901	<	DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1895	>	// Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1896	>	DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
1897	>	DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
1898	>	int *bitmask = (int *) bitmaskArray->data;              // get the previously made mask array
1899	>
1900	>	//Use conf_disambig map as a threshold on spaceweather quantities
1901	>	DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1902		DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1903		int *mask = (int *) maskArray->data;            // get the previously made mask array
1904	<
1904	>
1905		DRMS_Segment_t *bxSeg = drms_segment_lookup(inRec, BP_SEG_CEA);
1906		DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
1907		float *bx = (float *) bxArray->data;            // bx
#	Line 1665 | Line 1914 | void computeSWIndex(struct swIndex *swKe
1914		DRMS_Segment_t *bzSeg = drms_segment_lookup(inRec, BR_SEG_CEA);
1915		DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
1916		float *bz = (float *) bzArray->data;            // bz
1917	+
1918	+	//Use magnetogram map to compute R
1919	+	DRMS_Segment_t *losSeg = drms_segment_lookup(inRec, "magnetogram");
1920	+	DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
1921	+	float *los = (float *) losArray->data;          // los
1922	+
1923	+	DRMS_Segment_t *bz_errSeg = drms_segment_lookup(inRec, BR_ERR_SEG_CEA);
1924	+	DRMS_Array_t *bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status);
1925	+	float *bz_err = (float *) bz_errArray->data;            // bz_err
1926	+
1927	+	DRMS_Segment_t *by_errSeg = drms_segment_lookup(inRec, BT_ERR_SEG_CEA);
1928	+	DRMS_Array_t *by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status);
1929	+	float *by_err = (float *) by_errArray->data;            // by_err
1930	+	//for (int i = 0; i < nxny; i++) by_err[i] *= -1;
1931	+
1932	+	DRMS_Segment_t *bx_errSeg = drms_segment_lookup(inRec, BP_ERR_SEG_CEA);
1933	+	DRMS_Array_t *bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status);
1934	+	float *bx_err = (float *) bx_errArray->data;            // bx_err
1935
1936		// Get emphemeris
1937	<
1938	<	//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1939	<	float cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);
1940	<	float dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1941	<	double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);
1942	<	double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);
1943	<	float imcrpix1    = drms_getkey_float(inRec, "IMCRPIX1", &status);
1944	<	float imcrpix2    = drms_getkey_float(inRec, "IMCRPIX2", &status);
1945	<	float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);
1946	<	float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);
1947	<
1948	<	//float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)
1949	<
1950	<	printf("cdelt1=%f\n",cdelt1);
1951	<	printf("rsun_ref=%f\n",rsun_ref);
1952	<	printf("rsun_obs=%f\n",rsun_obs);
1953	<	printf("dsun_obs=%f\n",dsun_obs);
1954	<
1955	<	// Temp arrays
1956	<
1957	<	float *bh = (float *) (malloc(nxny * sizeof(float)));
1958	<	float *bt = (float *) (malloc(nxny * sizeof(float)));
1959	<	float *jz = (float *) (malloc(nxny * sizeof(float)));
1960	<	float *bpx = (float *) (malloc(nxny * sizeof(float)));
1961	<	float *bpy = (float *) (malloc(nxny * sizeof(float)));
1962	<	float *bpz = (float *) (malloc(nxny * sizeof(float)));
1963	<	float *derx = (float *) (malloc(nxny * sizeof(float)));
1697	<	float *dery = (float *) (malloc(nxny * sizeof(float)));
1937	>	float  cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1938	>	float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1939	>	double rsun_ref    = drms_getkey_double(inRec, "RSUN_REF", &status);
1940	>	double rsun_obs    = drms_getkey_double(inRec, "RSUN_OBS", &status);
1941	>	float imcrpix1     = drms_getkey_float(inRec, "IMCRPIX1", &status);
1942	>	float imcrpix2     = drms_getkey_float(inRec, "IMCRPIX2", &status);
1943	>	float crpix1       = drms_getkey_float(inRec, "CRPIX1", &status);
1944	>	float crpix2       = drms_getkey_float(inRec, "CRPIX2", &status);
1945	>
1946	>	// convert cdelt1_orig from degrees to arcsec
1947	>	float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
1948	>
1949	>	//if (nx1 > floor((nx-1)/scale + 1) )
1950	>	//      DIE("X-dimension of output array in fsample() is too large.");
1951	>	//if (ny1 > floor((ny-1)/scale + 1) )
1952	>	//      DIE("Y-dimension of output array in fsample() is too large.");
1953	>
1954	>	// Temp arrays
1955	>	float *bh      = (float *) (malloc(nxny * sizeof(float)));
1956	>	float *bt      = (float *) (malloc(nxny * sizeof(float)));
1957	>	float *jz      = (float *) (malloc(nxny * sizeof(float)));
1958	>	float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
1959	>	float *bpx     = (float *) (malloc(nxny * sizeof(float)));
1960	>	float *bpy     = (float *) (malloc(nxny * sizeof(float)));
1961	>	float *bpz     = (float *) (malloc(nxny * sizeof(float)));
1962	>	float *derx    = (float *) (malloc(nxny * sizeof(float)));
1963	>	float *dery    = (float *) (malloc(nxny * sizeof(float)));
1964		float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
1965		float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
1966		float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
1967		float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
1968		float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1969		float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1970	<
1971	<	// Compute
1972	<
1973	<	if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1974	<	mask, cdelt1, rsun_ref, rsun_obs)){
1970	>	float *bt_err  = (float *) (malloc(nxny * sizeof(float)));
1971	>	float *bh_err  = (float *) (malloc(nxny * sizeof(float)));
1972	>	float *jz_err  = (float *) (malloc(nxny * sizeof(float)));
1973	>	float *jz_err_squared = (float *) (malloc(nxny * sizeof(float)));
1974	>	float *jz_err_squared_smooth = (float *) (malloc(nxny * sizeof(float)));
1975	>	float *jz_rms_err  = (float *) (malloc(nxny * sizeof(float)));
1976	>	float *err_term1   = (float *) (calloc(nxny, sizeof(float)));
1977	>	float *err_term2   = (float *) (calloc(nxny, sizeof(float)));
1978	>	float *err_termA   = (float *) (calloc(nxny, sizeof(float)));
1979	>	float *err_termB   = (float *) (calloc(nxny, sizeof(float)));
1980	>	float *err_termAt  = (float *) (calloc(nxny, sizeof(float)));
1981	>	float *err_termBt  = (float *) (calloc(nxny, sizeof(float)));
1982	>	float *err_termAh  = (float *) (calloc(nxny, sizeof(float)));
1983	>	float *err_termBh  = (float *) (calloc(nxny, sizeof(float)));
1984	>
1985	>	// define some values for the R calculation
1986	>	int scale = round(2.0/cdelt1);
1987	>	int nx1 = nx/scale;
1988	>	int ny1 = ny/scale;
1989	>	int nxp = nx1+40;
1990	>	int nyp = ny1+40;
1991	>	float *rim     = (float *)malloc(nx1*ny1*sizeof(float));
1992	>	float *p1p0    = (float *)malloc(nx1*ny1*sizeof(float));
1993	>	float *p1n0    = (float *)malloc(nx1*ny1*sizeof(float));
1994	>	float *p1p     = (float *)malloc(nx1*ny1*sizeof(float));
1995	>	float *p1n     = (float *)malloc(nx1*ny1*sizeof(float));
1996	>	float *p1      = (float *)malloc(nx1*ny1*sizeof(float));
1997	>	float *pmap    = (float *)malloc(nxp*nyp*sizeof(float));
1998	>	float *p1pad   = (float *)malloc(nxp*nyp*sizeof(float));
1999	>	float *pmapn   = (float *)malloc(nx1*ny1*sizeof(float));
2000	>
2001	>	// define some arrays for the lorentz force calculation
2002	>	float *fx = (float *) (malloc(nxny * sizeof(float)));
2003	>	float *fy = (float *) (malloc(nxny * sizeof(float)));
2004	>	float *fz = (float *) (malloc(nxny * sizeof(float)));
2005	>
2006	>
2007	>	//spaceweather quantities computed
2008	>	if (computeAbsFlux(bz_err, bz , dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),  &(swKeys_ptr->mean_vf_err),
2009	>	&(swKeys_ptr->count_mask), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2010	>	{
2011		swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;               // If fail, fill in NaN
2012		swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
2013	+	swKeys_ptr->mean_vf_err = DRMS_MISSING_FLOAT;
2014	+	swKeys_ptr->count_mask  = DRMS_MISSING_INT;
2015		}
2016	<
2016	>
2017		for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
2018	<	greenpot(bpx, bpy, bpz, nx, ny);
2018	>	greenpot(bpx, bpy, bpz, nx, ny);
2019
2020	<	computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask);
2021	<
2022	<	if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask))
2023	<	swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
2024	<
2025	<	computeB_total(bx, by, bz, bt, dims, mask);
2026	<
2027	<	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt))
2020	>	computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
2021	>
2022	>	if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), &(swKeys_ptr->mean_gamma_err),mask, bitmask))
2023	>	{
2024	>	swKeys_ptr->mean_gamma     =  DRMS_MISSING_FLOAT;
2025	>	swKeys_ptr->mean_gamma_err =  DRMS_MISSING_FLOAT;
2026	>	}
2027	>
2028	>	computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask);
2029	>
2030	>	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt,
2031	>	dery_bt, bt_err, &(swKeys_ptr->mean_derivative_btotal_err), err_termAt, err_termBt))
2032	>	{
2033		swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
2034	+	swKeys_ptr->mean_derivative_btotal_err = DRMS_MISSING_FLOAT;
2035	+	}
2036
2037	<	if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh))
2037	>	if (computeBhderivative(bh, bh_err, dims, &(swKeys_ptr->mean_derivative_bh),
2038	>	&(swKeys_ptr->mean_derivative_bh_err), mask, bitmask, derx_bh, dery_bh, err_termAh, err_termBh))
2039	>	{
2040		swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
2041	<
2042	<	if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, derx_bz, dery_bz))
2041	>	swKeys_ptr->mean_derivative_bh_err = DRMS_MISSING_FLOAT;
2042	>	}
2043	>
2044	>	if (computeBzderivative(bz, bz_err, dims, &(swKeys_ptr->mean_derivative_bz), &(swKeys_ptr->mean_derivative_bz_err),
2045	>	mask, bitmask, derx_bz, dery_bz, err_termA, err_termB))
2046	>	{
2047		swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2048	+	swKeys_ptr->mean_derivative_bz_err = DRMS_MISSING_FLOAT;
2049	+	}
2050
2051	<
2052	<
2053	<	if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask,
2054	<	cdelt1, rsun_ref, rsun_obs, derx, dery)) {
2055	<	swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
2056	<	swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
2051	>	computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs,
2052	>	derx, dery, err_term1, err_term2);
2053	>
2054	>
2055	>	if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &(swKeys_ptr->mean_jz),
2056	>	&(swKeys_ptr->mean_jz_err), &(swKeys_ptr->us_i), &(swKeys_ptr->us_i_err), mask, bitmask, cdelt1,
2057	>	rsun_ref, rsun_obs, derx, dery))
2058	>	{
2059	>	swKeys_ptr->mean_jz            = DRMS_MISSING_FLOAT;
2060	>	swKeys_ptr->us_i               = DRMS_MISSING_FLOAT;
2061	>	swKeys_ptr->mean_jz_err        = DRMS_MISSING_FLOAT;
2062	>	swKeys_ptr->us_i_err           = DRMS_MISSING_FLOAT;
2063		}
2064	<
2065	<	printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
2066	<
2067	<	if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, cdelt1, rsun_ref, rsun_obs))
2068	<	swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
2069	<
2070	<	if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih),
2071	<	&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
2072	<	mask, cdelt1, rsun_ref, rsun_obs)) {
2073	<	swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
2074	<	swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
2075	<	swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
2064	>
2065	>	if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &(swKeys_ptr->mean_alpha), &(swKeys_ptr->mean_alpha_err),
2066	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2067	>	{
2068	>	swKeys_ptr->mean_alpha         = DRMS_MISSING_FLOAT;
2069	>	swKeys_ptr->mean_alpha_err     = DRMS_MISSING_FLOAT;
2070	>	}
2071	>
2072	>	if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &(swKeys_ptr->mean_ih), &(swKeys_ptr->mean_ih_err),
2073	>	&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
2074	>	&(swKeys_ptr->total_us_ih_err), &(swKeys_ptr->total_abs_ih_err), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2075	>	{
2076	>	swKeys_ptr->mean_ih            = DRMS_MISSING_FLOAT;
2077	>	swKeys_ptr->total_us_ih        = DRMS_MISSING_FLOAT;
2078	>	swKeys_ptr->total_abs_ih       = DRMS_MISSING_FLOAT;
2079	>	swKeys_ptr->mean_ih_err        = DRMS_MISSING_FLOAT;
2080	>	swKeys_ptr->total_us_ih_err    = DRMS_MISSING_FLOAT;
2081	>	swKeys_ptr->total_abs_ih_err   = DRMS_MISSING_FLOAT;
2082		}
2083	<
2084	<	if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),
2085	<	mask, cdelt1, rsun_ref, rsun_obs))
2086	<	swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
2087	<
2088	<
1758	<	if (computeFreeEnergy(bx, by, bpx, bpy, dims,
1759	<	&(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
1760	<	mask, cdelt1, rsun_ref, rsun_obs)) {
1761	<	swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1762	<	swKeys_ptr->totpot = DRMS_MISSING_FLOAT;
2083	>
2084	>	if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &(swKeys_ptr->totaljz), &(swKeys_ptr->totaljz_err),
2085	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2086	>	{
2087	>	swKeys_ptr->totaljz            = DRMS_MISSING_FLOAT;
2088	>	swKeys_ptr->totaljz_err        = DRMS_MISSING_FLOAT;
2089		}
2090
2091	<	if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,
2092	<	&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),
2093	<	&(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
2094	<	mask)) {
2095	<	swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2091	>	if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims,
2092	>	&(swKeys_ptr->meanpot), &(swKeys_ptr->meanpot_err), &(swKeys_ptr->totpot), &(swKeys_ptr->totpot_err),
2093	>	mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2094	>	{
2095	>	swKeys_ptr->meanpot            = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2096	>	swKeys_ptr->totpot             = DRMS_MISSING_FLOAT;
2097	>	swKeys_ptr->meanpot_err        = DRMS_MISSING_FLOAT;
2098	>	swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;
2099	>	}
2100	>
2101	>
2102	>	if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,
2103	>	&(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),
2104	>	mask, bitmask))
2105	>	{
2106	>	swKeys_ptr->meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2107		swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
2108	<	swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
2109	<	swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;
2108	>	swKeys_ptr->meanshear_angle_err= DRMS_MISSING_FLOAT;
2109	>	}
2110	>
2111	>	if (computeR(bz_err, los , dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,
2112	>	p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
2113	>	{
2114	>	swKeys_ptr->Rparam = DRMS_MISSING_FLOAT;                // If fail, fill in NaN
2115	>	}
2116	>
2117	>
2118	>	if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &(swKeys_ptr->totfx), &(swKeys_ptr->totfy), &(swKeys_ptr->totfz), &(swKeys_ptr->totbsq),
2119	>	&(swKeys_ptr->epsx), &(swKeys_ptr->epsy), &(swKeys_ptr->epsz), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
2120	>	{
2121	>	swKeys_ptr->totfx             = DRMS_MISSING_FLOAT;
2122	>	swKeys_ptr->totfy             = DRMS_MISSING_FLOAT;
2123	>	swKeys_ptr->totfz             = DRMS_MISSING_FLOAT;
2124	>	swKeys_ptr->totbsq            = DRMS_MISSING_FLOAT;
2125	>	swKeys_ptr->epsx              = DRMS_MISSING_FLOAT;
2126	>	swKeys_ptr->epsy              = DRMS_MISSING_FLOAT;
2127	>	swKeys_ptr->epsz              = DRMS_MISSING_FLOAT;
2128	>
2129		}
1774	–
1775	–	// Clean up
2130
2131	+	// Clean up the arrays
2132	+
2133	+	drms_free_array(bitmaskArray);          // Dec 18 2012 Xudong
2134		drms_free_array(maskArray);
2135		drms_free_array(bxArray);
2136		drms_free_array(byArray);
2137		drms_free_array(bzArray);
2138	+	drms_free_array(losArray);              // Mar 7
2139	+	drms_free_array(bx_errArray);
2140	+	drms_free_array(by_errArray);
2141	+	drms_free_array(bz_errArray);
2142
2143	<	free(bh); free(bt); free(jz);
2143	>	free(bh); free(bt); free(jz); free(jz_smooth);
2144		free(bpx); free(bpy); free(bpz);
2145	<	free(derx); free(dery);
2146	<	free(derx_bt); free(dery_bt);
2147	<	free(derx_bz); free(dery_bz);
2145	>	free(derx); free(dery);
2146	>	free(derx_bt); free(dery_bt);
2147	>	free(derx_bz); free(dery_bz);
2148		free(derx_bh); free(dery_bh);
2149	<
2150	<	}
2149	>	free(bt_err); free(bh_err);  free(jz_err);
2150	>	free(jz_err_squared); free(jz_rms_err);
2151	>	free(jz_err_squared_smooth);
2152	>
2153	>	// free the arrays that are related to the numerical derivatives
2154	>	free(err_term2);
2155	>	free(err_term1);
2156	>	free(err_termB);
2157	>	free(err_termA);
2158	>	free(err_termBt);
2159	>	free(err_termAt);
2160	>	free(err_termBh);
2161	>	free(err_termAh);
2162	>
2163	>	// free the arrays that are related to the r calculation
2164	>	free(rim);
2165	>	free(p1p0);
2166	>	free(p1n0);
2167	>	free(p1p);
2168	>	free(p1n);
2169	>	free(p1);
2170	>	free(pmap);
2171	>	free(p1pad);
2172	>	free(pmapn);
2173
2174	+	// free the arrays that are related to the lorentz calculation
2175	+	free(fx); free(fy); free(fz);
2176	+	}
2177
2178	<	/*
2178	>	/*
2179		* Set space weather indices, no error checking for now
2180		*
2181		*/
2182
2183		void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
2184		{
2185	<	drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
2186	<	drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
2187	<	drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
2188	<	drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
2189	<	drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
2190	<	drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
2191	<	drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
2192	<	drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
2193	<	drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
2194	<	drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
2195	<	drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
2196	<	drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
2197	<	drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
2198	<	drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);
2199	<	drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
2200	<	drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
2185	>	drms_setkey_float(outRec, "USFLUX",  swKeys_ptr->mean_vf);
2186	>	drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
2187	>	drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
2188	>	drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
2189	>	drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
2190	>	drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
2191	>	drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
2192	>	drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
2193	>	drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
2194	>	drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
2195	>	drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
2196	>	drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
2197	>	drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
2198	>	drms_setkey_float(outRec, "TOTPOT",  swKeys_ptr->totpot);
2199	>	drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
2200	>	drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
2201	>	drms_setkey_float(outRec, "CMASK",   swKeys_ptr->count_mask);
2202	>	drms_setkey_float(outRec, "ERRBT",   swKeys_ptr->mean_derivative_btotal_err);
2203	>	drms_setkey_float(outRec, "ERRVF",   swKeys_ptr->mean_vf_err);
2204	>	drms_setkey_float(outRec, "ERRGAM",  swKeys_ptr->mean_gamma_err);
2205	>	drms_setkey_float(outRec, "ERRBH",   swKeys_ptr->mean_derivative_bh_err);
2206	>	drms_setkey_float(outRec, "ERRBZ",   swKeys_ptr->mean_derivative_bz_err);
2207	>	drms_setkey_float(outRec, "ERRJZ",   swKeys_ptr->mean_jz_err);
2208	>	drms_setkey_float(outRec, "ERRUSI",  swKeys_ptr->us_i_err);
2209	>	drms_setkey_float(outRec, "ERRALP",  swKeys_ptr->mean_alpha_err);
2210	>	drms_setkey_float(outRec, "ERRMIH",  swKeys_ptr->mean_ih_err);
2211	>	drms_setkey_float(outRec, "ERRTUI",  swKeys_ptr->total_us_ih_err);
2212	>	drms_setkey_float(outRec, "ERRTAI",  swKeys_ptr->total_abs_ih_err);
2213	>	drms_setkey_float(outRec, "ERRJHT",  swKeys_ptr->totaljz_err);
2214	>	drms_setkey_float(outRec, "ERRMPOT", swKeys_ptr->meanpot_err);
2215	>	drms_setkey_float(outRec, "ERRTPOT", swKeys_ptr->totpot_err);
2216	>	drms_setkey_float(outRec, "ERRMSHA", swKeys_ptr->meanshear_angle_err);
2217	>	drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);
2218	>	drms_setkey_float(outRec, "TOTFX",   swKeys_ptr->totfx);
2219	>	drms_setkey_float(outRec, "TOTFY",   swKeys_ptr->totfy);
2220	>	drms_setkey_float(outRec, "TOTFZ",   swKeys_ptr->totfz);
2221	>	drms_setkey_float(outRec, "TOTBSQ",  swKeys_ptr->totbsq);
2222	>	drms_setkey_float(outRec, "EPSX",    swKeys_ptr->epsx);
2223	>	drms_setkey_float(outRec, "EPSY",    swKeys_ptr->epsy);
2224	>	drms_setkey_float(outRec, "EPSZ",    swKeys_ptr->epsz);
2225		};
2226
2227	<
1818	<	/*
2227	>	/*
2228		* Set all keywords, no error checking for now
2229		*
2230		*/
2231
2232	<	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)
2232	>	void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
2233		{
2234	<	copy_me_keys(inRec, outRec);
2235	<	copy_patch_keys(inRec, outRec);
2236	<	copy_geo_keys(inRec, outRec);
2237	<	copy_ambig_keys(inRec, outRec);
2238	<
2239	<	char timebuf[1024];
2240	<	float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
2241	<	double val;
2242	<	int status = DRMS_SUCCESS;
2243	<
2244	<	val = drms_getkey_double(inRec, "DATE",&status);
2245	<	drms_setkey_double(outRec, "DATE_B", val);
2246	<	sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
2247	<	drms_setkey_string(outRec, "DATE", timebuf);
2248	<
2249	<	// set cvs commit version into keyword HEADER
2250	<	char *cvsinfo = strdup("$Header$");
2251	<	status = drms_setkey_string(outRec, "HEADER", cvsinfo);
2252	<
2234	>
2235	>	copy_me_keys(bharpRec, outRec);
2236	>	copy_patch_keys(mharpRec, outRec);      // Dec 30
2237	>	copy_geo_keys(mharpRec, outRec);        // Dec 30
2238	>	copy_ambig_keys(bharpRec, outRec);
2239	>
2240	>	int status = 0;
2241	>
2242	>	// Change a few geometry keywords for CEA & cutout records
2243	>	if (mInfo != NULL) {        // CEA
2244	>
2245	>	drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
2246	>	drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
2247	>
2248	>	drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
2249	>	drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
2250	>	drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
2251	>	drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
2252	>	drms_setkey_string(outRec, "CUNIT1", "degree");
2253	>	drms_setkey_string(outRec, "CUNIT2", "degree");
2254	>
2255	>	char key[64];
2256	>	snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
2257	>	drms_setkey_string(outRec, "CTYPE1", key);
2258	>	snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
2259	>	drms_setkey_string(outRec, "CTYPE2", key);
2260	>	drms_setkey_float(outRec, "CROTA2", 0.0);
2261	>
2262	>	// Jan 2 2014 XS
2263	>	int nSeg = ARRLENGTH(CEASegs);
2264	>	for (int iSeg = 0; iSeg < nSeg; iSeg++) {
2265	>	DRMS_Segment_t *outSeg = NULL;
2266	>	outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);
2267	>	if (!outSeg) continue;
2268	>	// Set Bunit
2269	>	char bunit_xxx[20];
2270	>	sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
2271	>	//printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);
2272	>	drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);
2273	>	}
2274	>
2275	>	} else {        // Cutout
2276	>
2277	>	float disk_xc, disk_yc;
2278	>	if (fullDisk) {
2279	>	disk_xc = drms_getkey_float(bharpRec, "CRPIX1", &status);
2280	>	disk_yc = drms_getkey_float(bharpRec, "CRPIX2", &status);
2281	>	} else {
2282	>	disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);
2283	>	disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);
2284	>	}
2285	>	float x_ll = drms_getkey_float(mharpRec, "CRPIX1", &status);
2286	>	float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &status);
2287	>	// Defined as disk center's pixel address wrt lower-left of cutout
2288	>	drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
2289	>	drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
2290	>	// Always 0.
2291	>	drms_setkey_float(outRec, "CRVAL1", 0);
2292	>	drms_setkey_float(outRec, "CRVAL2", 0);
2293	>
2294	>	// Jan 2 2014 XS
2295	>	int nSeg = ARRLENGTH(CutSegs);
2296	>	for (int iSeg = 0; iSeg < nSeg; iSeg++) {
2297	>	DRMS_Segment_t *outSeg = NULL;
2298	>	outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);
2299	>	if (!outSeg) continue;
2300	>	// Set Bunit
2301	>	char bunit_xxx[20];
2302	>	sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
2303	>	//printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);
2304	>	drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);
2305	>	}
2306	>
2307	>
2308	>	}
2309	>
2310	>	// Mar 19 XS
2311	>	if (fullDisk) {
2312	>	drms_setkey_int(outRec, "AMBPATCH", 0);
2313	>	drms_setkey_int(outRec, "AMBWEAK", 2);
2314	>	} else {
2315	>	drms_setkey_int(outRec, "AMBPATCH", 1);
2316	>	}
2317	>
2318	>	TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
2319	>	tnow = (double)time(NULL);
2320	>	tnow += UNIX_epoch;
2321	>
2322	>	val = drms_getkey_time(bharpRec, "DATE", &status);
2323	>	drms_setkey_time(outRec, "DATE_B", val);
2324	>	drms_setkey_time(outRec, "DATE", tnow);
2325	>
2326	>	// set cvs commit version into keyword HEADER
2327	>	char *cvsinfo  = strdup("$Id$");
2328	>	char *cvsinfo2 = sw_functions_version();
2329	>	char cvsinfoall[2048];
2330	>	strcat(cvsinfoall,cvsinfo);
2331	>	strcat(cvsinfoall,"\n");
2332	>	strcat(cvsinfoall,cvsinfo2);
2333	>	status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
2334	>
2335		};
2336
1846	–	//
1847	–	//
2337
2338		/* ############# Nearest neighbour interpolation ############### */
2339
#	Line 1861 | Line 2350 | float nnb (float *f, int nx, int ny, dou
2350
2351		}
2352
2353	+
2354		/* ################## Wrapper for Jesper's rebin code ################## */
2355
2356		void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)

Diff Legend

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