(file) Return to sharp.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

   1 xudong 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             *
  16             *	Author:
  17             *		Xudong Sun; Monica Bobra
  18             *
  19             *	Version:
  20             *		v0.0	Jul 02 2012
  21             *		v0.1	Jul 23 2012
  22 xudong 1.3  *              v0.2    Sep 04 2012
  23 xudong 1.1  *
  24             *	Notes:
  25             *		v0.0
  26             *		Mharp & Bharp must be fully specified; other input are series names only
  27             *		All input records need to match, otherwise quit
  28             *		Mapping parameters depend on keywords of each record only, not necessarily consistent for now
  29             *		Cutout doesn't work for char segments yet (drms_segment_readslice bug)
  30             *		SW indices require ephemeris info which is not passed properly as of now
  31             *		v0.1
  32             *		Fixed char I/O thanks to Art
  33             *		SW indices fixed
  34             *		Added doppler and continuum
  35 xudong 1.3  *              Added other keywords: HEADER (populated by cvs build version), DATE_B
  36 xudong 1.1  *
  37             *	Example:
  38             *	sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
  39                      "bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
  40            		  "dop=hmi.V_720s[2012.02.20_10:00]" \
  41            		  "cont=hmi.Ic_720s[2012.02.20_10:00]" \
  42                      "sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
  43             *	For comparison:
  44             *	bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
  45            		 "out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
  46             *
  47             *
  48             */
  49                                                                           
  50            #include <stdio.h>
  51            #include <stdlib.h>
  52            #include <time.h>
  53            #include <sys/time.h>
  54            #include <math.h>
  55            #include <string.h>
  56            #include "jsoc_main.h"
  57 xudong 1.1 #include "astro.h"
  58            #include "fstats.h"
  59            #include "cartography.c"
  60            #include "fresize.h"
  61            #include "finterpolate.h"
  62            #include "img2helioVector.c"
  63            #include "copy_me_keys.c"
  64            #include "errorprop.c"
  65            #include "sw_functions.c"
  66            
  67            #include <mkl_blas.h>
  68            #include <mkl_service.h>
  69            #include <mkl_lapack.h>
  70            #include <mkl_vml_functions.h>
  71            #include <omp.h>
  72            
  73            #define PI              (M_PI)
  74            #define RADSINDEG		(PI/180.)
  75            #define RAD2ARCSEC		(648000./M_PI)
  76            #define SECINDAY		(86400.)
  77            #define FOURK			(4096)
  78 xudong 1.1 #define FOURK2    (16777216)
  79            
  80            #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
  81            
  82            // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
  83            #define NYQVIST		(0.015)
  84            
  85            #ifndef MIN
  86            #define MIN(a,b) (((a)<(b)) ? (a) : (b))
  87            #endif
  88            #ifndef MAX
  89            #define MAX(a,b) (((a)>(b)) ? (a) : (b))
  90            #endif
  91            
  92            #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
  93            #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
  94            
  95            #define kNotSpecified "Not Specified"
  96                                                    
  97            // Macros for WCS transformations.  assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
  98            // and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
  99 xudong 1.1 // PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun from disk center.
 100            #define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
 101            #define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
 102            #define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
 103            #define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
 104            
 105            #define DISAMB_AZI		1
 106            #define XSCALE			0.03
 107            #define YSCALE			0.03
 108            #define NBIN			3
 109            #define INTERP			0
 110            #define dpath    "/home/jsoc/cvs/Development/JSOC"
 111            
 112            /* ========================================================================================================== */
 113            
 114            // Space weather keywords
 115            struct swIndex {
 116            	float mean_vf; 
 117            	float absFlux;
 118            	float mean_hf;
 119            	float mean_gamma;
 120 xudong 1.1 	float mean_derivative_btotal;
 121            	float mean_derivative_bh; 
 122            	float mean_derivative_bz;
 123            	float mean_jz;
 124            	float us_i;
 125            	float mean_alpha;
 126            	float mean_ih;
 127            	float total_us_ih;
 128            	float total_abs_ih;
 129            	float totaljz;
 130            	float totpot;
 131            	float meanpot;
 132            	float area_w_shear_gt_45;
 133            	float meanshear_angle;
 134            	float area_w_shear_gt_45h;
 135            	float meanshear_angleh; 
 136            };
 137            
 138            // Mapping method
 139            enum projection {
 140            	carree,
 141 xudong 1.1 	cassini,
 142            	mercator,
 143            	cyleqa,
 144            	sineqa,
 145            	gnomonic,
 146            	postel,
 147            	stereographic,
 148            	orthographic,
 149            	lambert
 150            };
 151            
 152            // Ephemeris
 153            struct ephemeris {
 154            	double disk_lonc, disk_latc;
 155            	double disk_xc, disk_yc;
 156            	double rSun, asd, pa;
 157            };
 158            
 159            // Mapping information
 160            struct mapInfo {
 161            	float xc, yc;		// reference point: center
 162 xudong 1.1 	int nrow, ncol;		// size
 163            	float xscale, yscale;	// scale
 164            	int nbin;
 165            	enum projection proj;	// projection method
 166            	struct ephemeris ephem;		// ephemeris info
 167            	float *xi_out, *zeta_out;	// coordinate on full disk image to sample at
 168            };
 169            
 170            
 171            /* ========================================================================================================== */
 172            
 173            /* Get all input data series */
 174            int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
 175            			   char *mharpQuery, char *bharpQuery);
 176            
 177            /* Check if Mharp and Bharp match */
 178            int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS);
 179            
 180            /* Get other data series */
 181            int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS);
 182            
 183 xudong 1.1 /* Find record from record set with given T_rec */
 184            int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
 185            
 186            // ===================
 187            
 188            /* Create CEA record */
 189            int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
 190            					DRMS_Record_t *dopRec, DRMS_Record_t *contRec, 
 191            					DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
 192            
 193            /* Mapping single segment, wrapper */
 194            int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
 195            			  struct mapInfo *mInfo, char *segName);
 196            
 197            /* Mapping vector magnetogram, wrapper */
 198            int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
 199            
 200            /* Mapping errors of vector magnetogram, wraper */
 201            int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
 202            
 203            /* Determine reference point coordinate and patch size according to input */
 204 xudong 1.1 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo);
 205            
 206            /* Get ephemeris information */
 207            int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem);
 208            
 209            /* Compute the coordinates at which the full disk image is sampled */
 210            void findCoord(struct mapInfo *mInfo);
 211            
 212            /* Mapping function */
 213            int performSampling(float *outData, float *inData, struct mapInfo *mInfo);
 214            
 215            /* Performing local vector transformation */
 216            void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
 217            
 218            /* Map and propogate errors */
 219            int getBErr(float *bx_err, float *by_err, float *bz_err,
 220            			 DRMS_Record_t *inRec, struct mapInfo *mInfo);
 221            
 222            /* Read full disk vector magnetogram */
 223            int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
 224            
 225 xudong 1.1 /* Read variances and covariances of vector magnetograms */
 226            int readVectorBErr(DRMS_Record_t *bharpRec, 
 227            				   float *bT, float *bI, float *bA,
 228            				   float *errbT, float *errbI, float *errbA, 
 229            				   float *errbTbI, float *errbTbA, float *errbIbA);
 230            
 231            // ===================
 232            
 233            /* Create Cutout record */
 234            int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
 235            					DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
 236            					DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
 237            
 238            /* Get cutout and write segment */
 239            int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName);
 240            
 241            // ===================
 242            
 243            /* Compute space weather indices, no error checking for now */
 244            void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo);
 245            
 246 xudong 1.1 /* Set space weather indices, no error checking for now */
 247            void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
 248            
 249            /* Set all keywords, no error checking for now */
 250            void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);
 251            
 252            // ===================
 253            
 254            /* Nearest neighbor interpolation */
 255            float nnb (float *f, int nx, int ny, double x, double y);
 256            
 257            /* Wrapper for Jesper's rebin code */
 258            void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss);
 259            
 260            /* ========================================================================================================== */
 261            
 262            /* Remap segment names */
 263            #define BR_SEG_CEA		"Br"
 264            #define BT_SEG_CEA		"Bt"
 265            #define BP_SEG_CEA		"Bp"
 266            #define BR_ERR_SEG_CEA		"Br_err"
 267 xudong 1.1 #define BT_ERR_SEG_CEA		"Bt_err"
 268            #define BP_ERR_SEG_CEA		"Bp_err"
 269            
 270            /* Cutout segment names, input identical to output */
 271            char *MharpSegs[] = {"magnetogram", "bitmap"};
 272            char *BharpSegs[] = {"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
 273            	"damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
 274            	"conv_flag", // fixed
 275            	"info_map", "confid_map",
 276            	"inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
 277            	"field_inclination_err", "field_az_err", "inclin_azimuth_err",
 278            	"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
 279            	"disambig", "conf_disambig"};
 280            // For stats
 281            char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", 
 282            	"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
 283            	"damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
 284            	"conv_flag", // fixed
 285            	"info_map", "confid_map",
 286            	"inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
 287            	"field_inclination_err", "field_az_err", "inclin_azimuth_err",
 288 xudong 1.1 	"field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
 289            	"disambig", "conf_disambig"};
 290            char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig", 
 291            					BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
 292            /* ========================================================================================================== */
 293            
 294            
 295            
 296            char *module_name = "sharp";
 297            char *version_id = "2012 Jul 02";  /* Version number */
 298            
 299            ModuleArgs_t module_args[] =
 300            {
 301            	{ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
 302            	{ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},
 303            	{ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},
 304            	{ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
 305            	{ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
 306            	{ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
 307            	{ARG_END}
 308            };
 309 xudong 1.1 
 310            int DoIt(void)                    
 311            {
 312            	
 313            	int status = DRMS_SUCCESS;
 314            	int nrecs, irec;
 315            	
 316            	char *mharpQuery, *bharpQuery;
 317            	char *dopQuery, *contQuery;
 318            	char *sharpCeaQuery, *sharpCutQuery;
 319            	
 320            	DRMS_RecordSet_t *mharpRS = NULL, *bharpRS = NULL;
 321            	DRMS_RecordSet_t *dopRS = NULL, *contRS = NULL;
 322            	
 323            	/* Get parameters */
 324                
 325            	mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
 326            	bharpQuery = (char *) params_get_str(&cmdparams, "bharp");
 327            	dopQuery = (char *) params_get_str(&cmdparams, "dop");
 328            	contQuery = (char *) params_get_str(&cmdparams, "cont");
 329            	sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
 330 xudong 1.1 	sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
 331            	
 332            	/* Get input data, check everything */
 333            	
 334            	if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
 335            		DIE("Input harp data error.");
 336            	nrecs = mharpRS->n;
 337            	
 338            	if (getInputRS_aux(&dopRS, dopQuery, mharpRS)) 
 339            		DIE("Input doppler data error.");
 340            	
 341            	if (getInputRS_aux(&contRS, contQuery, mharpRS)) 
 342            		DIE("Input continuum data error.");
 343            	
 344            	/* Start */
 345            	
 346            	printf("==============\nStart. %d image(s) in total.\n", nrecs);
 347            	
 348            	for (irec = 0; irec < nrecs; irec++) {
 349            		
 350            		/* Records in work */
 351 xudong 1.1 		
 352            		DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;
 353            		mharpRec = mharpRS->records[irec];
 354            		bharpRec = bharpRS->records[irec];
 355            		
 356            		TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
 357            	
 358            		struct swIndex swKeys;
 359            		
 360            		DRMS_Record_t *dopRec = NULL, *contRec = NULL;
 361            		if (getInputRec_aux(&dopRec, dopRS, trec)) {
 362            			printf("Fetching Doppler failed, image #%d skipped.\n", irec);
 363            			continue;
 364            		}
 365            		if (getInputRec_aux(&contRec, contRS, trec)) {
 366            			printf("Fetching continuum failed, image #%d skipped.\n", irec);
 367            			continue;
 368            		}
 369            		
 370            		/* Create CEA record */
 371            		
 372 xudong 1.1 		DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
 373            		if (status) {		// if failed
 374            			printf("Creating CEA failed, image #%d skipped.\n", irec);
 375            			continue;
 376            		}
 377            		
 378            		if (createCeaRecord(mharpRec, bharpRec, dopRec, contRec, sharpCeaRec, &swKeys)) {		// do the work
 379            			printf("Creating CEA failed, image #%d skipped.\n", irec);
 380            			drms_close_record(sharpCeaRec, DRMS_FREE_RECORD);
 381            			continue;
 382            		}		// swKeys updated here
 383            		
 384            		drms_close_record(sharpCeaRec, DRMS_INSERT_RECORD);
 385            		
 386            		/* Create Cutout record */
 387            		
 388            		DRMS_Record_t *sharpCutRec = drms_create_record(drms_env, sharpCutQuery, DRMS_PERMANENT, &status);
 389            		if (status) {		// if failed
 390            			printf("Creating cutout failed, image #%d skipped.\n", irec);
 391            			continue;
 392            		}
 393 xudong 1.1 		
 394            		if (createCutRecord(mharpRec, bharpRec, dopRec, contRec, sharpCutRec, &swKeys)) {		// do the work
 395            			printf("Creating cutout failed, image #%d skipped.\n", irec);
 396            			drms_close_record(sharpCutRec, DRMS_FREE_RECORD);
 397            			continue;
 398            		}		// swKeys used here
 399            		
 400            		drms_close_record(sharpCutRec, DRMS_INSERT_RECORD);
 401            		
 402            		/* Done */
 403            		
 404            		printf("Image #%d done.\n", irec);
 405            		
 406            	} // irec
 407            	
 408            	drms_close_records(mharpRS, DRMS_FREE_RECORD);
 409            	drms_close_records(bharpRS, DRMS_FREE_RECORD);
 410            	
 411            	return 0;
 412            
 413            }	// DoIt
 414 xudong 1.1 
 415            
 416            // ===================================================================
 417            // ===================================================================
 418            // ===================================================================
 419            
 420            
 421            /*
 422             * Get input data series, including mHarp and bharp
 423             * Need all records to match, otherwise quit
 424             *
 425             */
 426            
 427            int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
 428            			   char *mharpQuery, char *bharpQuery)
 429            {
 430            	
 431            	int status = 0;
 432            	
 433            	*mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
 434                if (status || (*mharpRS_ptr)->n == 0) return 1;
 435 xudong 1.1 	
 436            	*bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
 437                if (status || (*bharpRS_ptr)->n == 0) return 1;
 438            
 439            	if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
 440            	
 441            	return 0;
 442            	
 443            }
 444            
 445            /*
 446             * Check if Mharp and Bharp match
 447             *
 448             */
 449            
 450            int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS)
 451            {
 452            	
 453            	int status = 0;
 454            	int nrecs = mharpRS->n;
 455            	
 456 xudong 1.1 	DRMS_Record_t *mharpRec_t = NULL, *bharpRec_t = NULL;		// temporary recs for utility
 457            	
 458                if (bharpRS->n != nrecs) {
 459            		return 1;		// return 1 if different
 460            	}
 461            	
 462            	for (int i = 0; i < nrecs; i++) {
 463            		mharpRec_t = mharpRS->records[i];
 464            		bharpRec_t = bharpRS->records[i];
 465            		if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) != 
 466            			 drms_getkey_int(bharpRec_t, "HARPNUM", &status)) ||
 467            			(drms_getkey_time(mharpRec_t, "T_REC", &status) != 
 468            			 drms_getkey_time(bharpRec_t, "T_REC", &status)))
 469            		{
 470            			return 1;
 471            		}
 472            	}
 473            	
 474            	return 0;
 475            	
 476            }
 477 xudong 1.1 
 478            /* 
 479             * Get other data series, check all T_REC are available
 480             *
 481             */
 482            
 483            int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS)
 484            {
 485            	
 486            	int status = 0;
 487            	
 488            	*inRS_ptr = drms_open_records(drms_env, inQuery, &status);
 489            	if (status || (*inRS_ptr)->n == 0) return status;
 490            	
 491            	// Check if all T_rec are available, need to match both ways
 492            	int n = harpRS->n, n0 = (*inRS_ptr)->n;
 493            
 494            	for (int i0 = 0; i0 < n0; i0++) {
 495            		DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
 496            		TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
 497            		TIME trec = 0;
 498 xudong 1.1 		for (int i = 0; i < n; i++) {
 499            			DRMS_Record_t *harpRec = harpRS->records[i];
 500            			trec = drms_getkey_time(harpRec, "T_REC", &status);
 501            			if (fabs(trec0 - trec) < 10) break;
 502            		}
 503            		if (fabs(trec0 - trec) >= 10) return 1;
 504            	}
 505            	
 506            	for (int i = 0; i < n; i++) {
 507            		DRMS_Record_t *harpRec = harpRS->records[i];
 508            		TIME trec = drms_getkey_time(harpRec, "T_REC", &status);
 509            		TIME trec0 = 0;
 510            		for (int i0 = 0; i0 < n0; i0++) {
 511            			DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
 512            			trec0 = drms_getkey_time(inRec, "T_REC", &status);
 513            			if (fabs(trec0 - trec) < 10) break;
 514            		}
 515            		if (fabs(trec0 - trec) >= 10) return 1;
 516            	}
 517            	
 518            	return 0;
 519 xudong 1.1 	
 520            }
 521            
 522            /* 
 523             * Find record from record set with given T_rec
 524             *
 525             */
 526            
 527            int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec)
 528            {
 529            	
 530            	int status = 0;
 531            	
 532            	int n = inRS->n;
 533            	for (int i = 0; i < n; i++) {
 534            		*inRec_ptr = inRS->records[i];
 535            		TIME trec0 = drms_getkey_time((*inRec_ptr), "T_REC", &status);
 536            		if (fabs(trec0 - trec) < 10) return 0;
 537            	}
 538            	
 539            	return 1;
 540 xudong 1.1 	
 541            }
 542            
 543            
 544            
 545            
 546            /*
 547             * Create CEA record: top level subroutine
 548             * Also compute all the space weather keywords here
 549             *
 550             */
 551            
 552            int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, 
 553            					DRMS_Record_t *dopRec, DRMS_Record_t *contRec, 
 554            					DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
 555            {
 556            	
 557            	int status = 0;
 558            	DRMS_Segment_t *inSeg;
 559            	DRMS_Array_t *inArray;
 560            	
 561 xudong 1.1 	struct mapInfo mInfo;
 562            	mInfo.proj = (enum projection) cyleqa;		// projection method
 563            	mInfo.xscale = XSCALE;
 564            	mInfo.yscale = YSCALE;
 565            	mInfo.nbin = NBIN;
 566            	
 567            	// Get ephemeris
 568            
 569 xudong 1.3 	if (getEphemeris(mharpRec, &(mInfo.ephem))) {
 570            		SHOW("CEA: get ephemeris error\n");
 571            		return 1;
 572            	}
 573 xudong 1.1 	
 574            	// Find position
 575            	
 576 xudong 1.3 	if (findPosition(mharpRec, &mInfo)) {
 577            		SHOW("CEA: find position error\n");
 578            		return 1;
 579            	}
 580 xudong 1.1 	
 581            	// Create xi_out, zeta_out array in mInfo:
 582            	// Coordinates to sample in original full disk image
 583            	
 584            	int ncol0, nrow0;		// oversampled map size
 585            	ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2;	// pad with nbin/2 on edge to avoid NAN
 586            	nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
 587            	
 588            	mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
 589            	mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
 590            	
 591            	findCoord(&mInfo);		// compute it here so it could be shared by the following 4 functions
 592            	
 593            	// Mapping single segment: Mharp, etc.
 594            	
 595 xudong 1.3 	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
 596            		SHOW("CEA: mapping magnetogram error\n");
 597            		return 1;
 598            	}
 599            	if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
 600            		SHOW("CEA: mapping bitmap error\n");
 601            		return 1;
 602            	}
 603 xudong 1.1 	printf("Magnetogram mapping done.\n");
 604            	
 605 xudong 1.3 	if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
 606            		SHOW("CEA: mapping dopplergram error\n");
 607            		return 1;
 608            	}
 609 xudong 1.1 	printf("Dopplergram mapping done.\n");
 610            	
 611 xudong 1.3 	if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
 612            		SHOW("CEA: mapping continuum error\n");
 613            		return 1;
 614            	}
 615 xudong 1.1 	printf("Intensitygram mapping done.\n");
 616            
 617 xudong 1.3 	if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
 618            		SHOW("CEA: mapping conf_disambig error\n");
 619            		return 1;
 620            	}
 621 xudong 1.2 	printf("Conf disambig mapping done.\n");
 622 xudong 1.1 
 623            	// Mapping vector B
 624            	
 625 xudong 1.3 	if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
 626            		SHOW("CEA: mapping vector B error\n");
 627            		return 1;
 628            	}
 629 xudong 1.1 	printf("Vector B mapping done.\n");
 630            	
 631            	// Mapping vector B errors
 632            	
 633 xudong 1.3 	if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) {
 634            		SHOW("CEA: mapping vector B uncertainty error\n");
 635            		return 1;
 636            	}
 637 xudong 1.1 	printf("Vector B error done.\n");
 638            	
 639            	// Keywords & Links
 640            	
 641            	drms_copykey(sharpRec, mharpRec, "T_REC");
 642            	drms_copykey(sharpRec, mharpRec, "HARPNUM");
 643            	
 644            	DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
 645            	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);
 648            	
 649            	setKeys(sharpRec, bharpRec);		// Set all other keywords
 650 xudong 1.3 	drms_copykey(sharpRec, mharpRec, "QUALITY");		// copied from los records
 651 xudong 1.1 
 652            	// Space weather
 653            	
 654            	computeSWIndex(swKeys_ptr, sharpRec, &mInfo);		// compute it!
 655            	printf("Space weather indices done.\n");
 656            	
 657            	setSWIndex(sharpRec, swKeys_ptr);	// Set space weather indices
 658            
 659            	// Stats
 660            	
 661            	int nCEASegs = ARRLENGTH(CEASegs);
 662            	for (int iSeg = 0; iSeg < nCEASegs; iSeg++) {
 663            		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
 664            		DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
 665            		int stat = set_statistics(outSeg, outArray, 1);
 666            //		printf("%d => %d\n", iSeg, stat);
 667            		drms_free_array(outArray);
 668            	}
 669            	
 670            	free(mInfo.xi_out);
 671            	free(mInfo.zeta_out);
 672 xudong 1.1 	return 0;
 673            	
 674            }
 675            
 676            
 677            /* 
 678             * Mapping a single segment
 679             * Read in full disk image, utilize mapImage for mapping
 680             * then write the segment out, segName same in in/out Rec
 681             *
 682             */
 683            
 684            int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
 685            			  struct mapInfo *mInfo, char *segName)
 686            {
 687            	
 688            	int status = 0;
 689            	int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
 690            	int dims[2] = {nx, ny};
 691            	
 692            	// Input full disk array
 693 xudong 1.1 	
 694            	DRMS_Segment_t *inSeg = NULL;
 695            	inSeg = drms_segment_lookup(inRec, segName);
 696            	if (!inSeg) return 1;
 697            	
 698            	DRMS_Array_t *inArray = NULL;
 699            	inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
 700            	if (!inArray) return 1;
 701            	
 702            	float *inData;
 703            	int xsz = inArray->axis[0], ysz = inArray->axis[1];
 704            	if ((xsz != FOURK) || (ysz != FOURK)) {		// for bitmap, make tmp full disk
 705            		float *inData0 = (float *) inArray->data;
 706            		inData = (float *) (calloc(FOURK2, sizeof(float)));
 707            		int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
 708            		int y0 = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1;
 709            		int ind_map;
 710            		for (int row = 0; row < ysz; row++) {
 711            			for (int col = 0; col < xsz; col++) {
 712            				ind_map = (row + y0) * FOURK + (col + x0);
 713            				inData[ind_map] = inData0[row * xsz + col];
 714 xudong 1.1 			}
 715            		}
 716            		drms_free_array(inArray); inArray = NULL;
 717            	} else {
 718            		inData = (float *) inArray->data;
 719            	}
 720            	
 721            	// Mapping
 722            	
 723            	float *map = (float *) (malloc(nxny * sizeof(float)));
 724            	if (performSampling(map, inData, mInfo))
 725            		{if (inArray) drms_free_array(inArray); free(map); return 1;}
 726            	
 727            	// Write out
 728            	
 729            	DRMS_Segment_t *outSeg = NULL;	
 730            	outSeg = drms_segment_lookup(sharpRec, segName);
 731            	if (!outSeg) return 1;
 732            	
 733            	DRMS_Type_t arrayType = outSeg->info->type;
 734            	DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
 735 xudong 1.1 	if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
 736            	
 737            	// convert to needed data type
 738            	
 739            	drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);
 740            	
 741            	outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
 742            	outArray->parent_segment = outSeg;
 743            	outArray->israw = 0;		// always compressed
 744            	outArray->bzero = outSeg->bzero;
 745            	outArray->bscale = outSeg->bscale;
 746            
 747            	status = drms_segment_write(outSeg, outArray, 0);
 748            	if (status) return 0;
 749            	
 750            	if (inArray) drms_free_array(inArray);
 751            	if (outArray) drms_free_array(outArray);
 752            	return 0;
 753            	
 754            }
 755            
 756 xudong 1.1 
 757            /* 
 758             * Mapping vector magnetogram
 759             *
 760             */
 761            
 762            int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
 763            {
 764            	
 765            	int status = 0;
 766            	int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
 767            	int dims[2] = {nx, ny};
 768            	
 769            	// Read in segments, filling factor assume to be 1
 770            	
 771            	float *bx_img = (float *) (malloc(FOURK2 * sizeof(float)));
 772            	float *by_img = (float *) (malloc(FOURK2 * sizeof(float)));
 773            	float *bz_img = (float *) (malloc(FOURK2 * sizeof(float)));
 774            	
 775            	if (readVectorB(bharpRec, bx_img, by_img, bz_img)) {
 776            		printf("Read full disk image error\n");
 777 xudong 1.1 		free(bx_img); free(by_img); free(bz_img);
 778            		return 1;
 779            	}
 780            	
 781            	// Mapping
 782            	
 783            	float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;	// intermediate maps, in CCD bxyz representation
 784            
 785            	bx_map = (float *) (malloc(nxny * sizeof(float)));
 786            	if (performSampling(bx_map, bx_img, mInfo))
 787            		{free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
 788            
 789            	by_map = (float *) (malloc(nxny * sizeof(float)));
 790            	if (performSampling(by_map, by_img, mInfo))
 791            		{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
 792            
 793            	bz_map = (float *) (malloc(nxny * sizeof(float)));
 794            	if (performSampling(bz_map, bz_img, mInfo))
 795            		{free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
 796            	
 797            	free(bx_img); free(by_img); free(bz_img);
 798 xudong 1.1 	
 799            	// Vector transform
 800            	
 801            	vectorTransform(bx_map, by_map, bz_map, mInfo);
 802            	
 803            	for (int i = 0; i < nxny; i++) by_map[i] *= -1;		// positive theta pointing south
 804            	
 805            	// Write out
 806            	
 807            	DRMS_Segment_t *outSeg;
 808            	DRMS_Array_t *outArray;
 809            	
 810            	float *data_prt[3] = {bx_map, by_map, bz_map};
 811            	char *segName[3] = {BP_SEG_CEA, BT_SEG_CEA, BR_SEG_CEA};
 812            	
 813            	for (int iSeg = 0; iSeg < 3; iSeg++) {
 814            		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
 815            		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
 816            		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
 817            		outArray->parent_segment = outSeg;
 818            		outArray->israw = 0;
 819 xudong 1.1 		outArray->bzero = outSeg->bzero;
 820            		outArray->bscale = outSeg->bscale;
 821            		status = drms_segment_write(outSeg, outArray, 0);
 822            		if (status) return 1;
 823            		drms_free_array(outArray);
 824            	}
 825            	
 826            	//
 827            	
 828            	return 0;
 829            	
 830            }
 831            
 832            
 833            /* 
 834             * Mapping vector magnetogram errors
 835             *
 836             */
 837            
 838            int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
 839            {
 840 xudong 1.1 	
 841            	int status = 0;
 842            	
 843            	int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
 844            	int dims[2] = {nx, ny};
 845            	
 846            	// Compute propogated errors, using nearest neighbour interpolation
 847            	
 848            	float *bx_err = (float *) (malloc(nxny * sizeof(float)));
 849            	float *by_err = (float *) (malloc(nxny * sizeof(float)));
 850            	float *bz_err = (float *) (malloc(nxny * sizeof(float)));
 851            	
 852            	if (getBErr(bx_err, by_err, bz_err, bharpRec, mInfo)) {
 853            		free(bx_err); free(by_err); free(bz_err);
 854            		return 1;
 855            	}
 856            
 857            	// Write out
 858            	
 859            	DRMS_Segment_t *outSeg;
 860            	DRMS_Array_t *outArray;
 861 xudong 1.1 	
 862            	float *data_prt[3] = {bx_err, by_err, bz_err};
 863            	char *segName[3] = {BP_ERR_SEG_CEA, BT_ERR_SEG_CEA, BR_ERR_SEG_CEA};
 864            	
 865            	for (int iSeg = 0; iSeg < 3; iSeg++) {
 866            		outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
 867            		outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
 868            		outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
 869            		outArray->parent_segment = outSeg;
 870            		outArray->israw = 0;
 871            		outArray->bzero = outSeg->bzero;
 872            		outArray->bscale = outSeg->bscale;
 873            		status = drms_segment_write(outSeg, outArray, 0);
 874            		if (status) return 1;
 875            		drms_free_array(outArray);
 876            	}
 877            	
 878            	//
 879            	
 880            	return 0;
 881            	
 882 xudong 1.1 }
 883            
 884            
 885            
 886            /* 
 887             * Determine reference point coordinate and patch size according to keywords
 888             * xc, yc are the coordinate of patch center, in degrees
 889             * ncol and nrow are the final size
 890             *
 891             */
 892            
 893            int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo)
 894            {
 895            	
 896            	int status = 0;
 897            	int harpnum = drms_getkey_int(inRec, "HARPNUM", &status);
 898            	TIME trec = drms_getkey_time(inRec, "T_REC", &status);
 899            	float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
 900            	
 901            	/* Center coord */
 902            	
 903 xudong 1.1 	float minlon = drms_getkey_float(inRec, "LONDTMIN", &status); if (status) return 1;		// Stonyhurst lon
 904            	float maxlon = drms_getkey_float(inRec, "LONDTMAX", &status); if (status) return 1;
 905            	float minlat = drms_getkey_float(inRec, "LATDTMIN", &status); if (status) return 1;
 906            	float maxlat = drms_getkey_float(inRec, "LATDTMAX", &status); if (status) return 1;
 907            	
 908            	// A bug fixer for HARP (per M. Turmon)
 909            	// When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
 910            	// Also keywords such as "SIZE" will be NaN
 911            	// We compute minlon & minlat then by
 912            	// LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
 913            	
 914            	float psize = drms_getkey_float(inRec, "SIZE", &status);
 915            	if (psize != psize) {
 916            		TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1;
 917            		double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
 918            		char firstRecQuery[100], t0_str[100];
 919            		sprint_time(t0_str, t0, "TAI", 0);
 920            		snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str);
 921            		DRMS_RecordSet_t *tmpRS = drms_open_records(drms_env, firstRecQuery, &status);
 922            		if (status || tmpRS->n != 1) return 1;
 923            		DRMS_Record_t *tmpRec = tmpRS->records[0];
 924 xudong 1.1 		double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1;
 925            		double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1;
 926            		minlon = minlon0 + (trec - t0) * omega / SECINDAY;
 927            		maxlon = maxlon0 + (trec - t0) * omega / SECINDAY;
 928            		printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon);
 929            	}
 930            	
 931            	mInfo->xc = (maxlon + minlon) / 2. + disk_lonc;
 932            	mInfo->yc = (maxlat + minlat) / 2.;
 933            	
 934            	/* Size */
 935            	
 936            	mInfo->ncol = round((maxlon - minlon) / mInfo->xscale);
 937            	mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
 938            	
 939            	return 0;
 940            	
 941            }
 942            
 943            
 944            /*
 945 xudong 1.1  * Fetch ephemeris info from a DRMS record
 946             * No error checking for now
 947             *
 948             */
 949            
 950            int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem)
 951            {
 952            	
 953            	int status = 0;
 954            	
 955            	float crota2 = drms_getkey_float(inRec, "CROTA2", &status);	// rotation
 956            	double sina = sin(crota2 * RADSINDEG); 
 957            	double cosa = cos(crota2 * RADSINDEG);
 958            	
 959            	ephem->pa = - crota2 * RADSINDEG;
 960            	ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG;
 961            	ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG;
 962            	
 963            	float crvalx = 0.0;
 964            	float crvaly = 0.0;
 965            	float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
 966 xudong 1.1 	float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
 967            	float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
 968                    printf("cdelt=%f\n",cdelt);
 969            	ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;		// Center of disk in pixel, starting at 0
 970            	ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
 971            	
 972            	float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status);
 973            	float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status);
 974            	if (status) rSun_ref = 6.96e8;
 975            	
 976            	ephem->asd = asin(rSun_ref/dSun);
 977            	ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt;
 978            	
 979            	return 0;
 980            	
 981            }
 982            
 983            
 984            /*
 985             * Compute the coordinates to be sampled on full disk image
 986             * mInfo->xi_out & mInfo->zeta_out
 987 xudong 1.1  * This is oversampled, its size is ncol0 & nrow0 as shown below
 988             *
 989             *
 990             */
 991            
 992            void findCoord(struct mapInfo *mInfo)
 993            {
 994            	
 995            	int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;	// pad with nbin/2 on edge to avoid NAN
 996            	int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
 997            	
 998            	float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG;		// oversampling resolution
 999            	float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG;		// in rad
1000            	
1001            	double lonc = mInfo->xc * RADSINDEG;	// in rad
1002            	double latc = mInfo->yc * RADSINDEG;
1003            	
1004            	double disk_lonc = (mInfo->ephem).disk_lonc;
1005            	double disk_latc = (mInfo->ephem).disk_latc;
1006            	
1007            	double rSun = (mInfo->ephem).rSun;
1008 xudong 1.1 	double disk_xc = (mInfo->ephem).disk_xc / rSun;
1009            	double disk_yc = (mInfo->ephem).disk_yc / rSun;
1010            	double pa = (mInfo->ephem).pa;
1011            	
1012            	// Temp pointers
1013            	
1014            	float *xi_out = mInfo->xi_out;
1015            	float *zeta_out = mInfo->zeta_out;
1016            	
1017            	// start
1018            	
1019            	double x, y;		// map coord
1020            	double lat, lon;	// helio coord
1021            	double xi, zeta;	// image coord (for one point)
1022            	
1023            	int ind_map;
1024            	
1025            	for (int row0 = 0; row0 < nrow0; row0++) {
1026            		for (int col0 = 0; col0 < ncol0; col0++) {
1027            			
1028            			ind_map = row0 * ncol0 + col0;
1029 xudong 1.1 			
1030            			x = (col0 + 0.5 - ncol0/2.) * xscale0;		// in rad
1031            			y = (row0 + 0.5 - nrow0/2.) * yscale0;
1032            			
1033            			/* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates. 
1034            			 * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
1035            			 * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
1036            			 * is the heliographic longitude and latitude of the map center. Both are in degree.    
1037            			 */
1038            			
1039            			if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1040            				xi_out[ind_map] = -1;
1041            				zeta_out[ind_map] = -1;
1042            				continue;
1043            			}
1044            			
1045            			/* map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the
1046            			 * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
1047            			 */
1048            			
1049            			if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta, 
1050 xudong 1.1 							disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1051            				xi_out[ind_map] = -1;
1052            				zeta_out[ind_map] = -1;
1053            				continue;
1054            			}
1055            			
1056            			xi_out[ind_map] = xi * rSun;
1057            			zeta_out[ind_map] = zeta * rSun;
1058            			
1059            		}
1060            	}
1061            	
1062            }
1063            
1064            
1065            /* 
1066             * Sampling function
1067             * oversampling by nbin, then binning using a Gaussian
1068             * save results in outData, always of float type
1069             *
1070             */
1071 xudong 1.1 
1072            int performSampling(float *outData, float *inData, struct mapInfo *mInfo)
1073            {
1074            	
1075            	int status = 0;
1076            	
1077            	int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;	// pad with nbin/2 on edge to avoid NAN
1078            	int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
1079            	
1080            	float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1081            	
1082            	float *xi_out = mInfo->xi_out;
1083            	float *zeta_out = mInfo->zeta_out;
1084            
1085            	// Interpolation
1086            	
1087            	struct fint_struct pars;
1088            	int interpOpt = INTERP;		// Use Wiener by default, 6 order, 1 constraint
1089            	
1090            	switch (interpOpt) {
1091            		case 0:			// Wiener, 6 order, 1 constraint
1092 xudong 1.1 			init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath);
1093            			break;
1094            		case 1:			// Cubic convolution
1095            			init_finterpolate_cubic_conv(&pars, 1., 3.);
1096            			break;
1097            		case 2:			// Bilinear
1098            			init_finterpolate_linear(&pars, 1.);
1099            			break;
1100            		default:
1101            			return 1;
1102            	}
1103            	
1104            	finterpolate(&pars, inData, xi_out, zeta_out, outData0, 
1105            				 FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1106            	
1107            	// Rebinning, smoothing
1108            	
1109            	frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);		// Gaussian
1110            	
1111            	//
1112            	
1113 xudong 1.1 	return 0;
1114            	
1115            }
1116            
1117            
1118            /* 
1119             * Performing local vector transformation
1120             *  xyz: z refers to vertical (radial) component, x EW (phi), y NS
1121             *
1122             */
1123            
1124            void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo)
1125            {
1126            	
1127            	int ncol = mInfo->ncol;
1128            	int nrow = mInfo->nrow;
1129            	
1130            	float xscale = mInfo->xscale * RADSINDEG;		// in rad
1131            	float yscale = mInfo->yscale * RADSINDEG;
1132            	
1133            	double lonc = mInfo->xc * RADSINDEG;	// in rad
1134 xudong 1.1 	double latc = mInfo->yc * RADSINDEG;
1135            	
1136            	double disk_lonc = (mInfo->ephem).disk_lonc;
1137            	double disk_latc = (mInfo->ephem).disk_latc;
1138            	
1139            	double rSun = (mInfo->ephem).rSun;
1140            	double disk_xc = (mInfo->ephem).disk_xc / rSun;
1141            	double disk_yc = (mInfo->ephem).disk_yc / rSun;
1142            	double pa = (mInfo->ephem).pa;
1143            	
1144            	int ind_map;
1145            	double x, y;
1146            	double lat, lon;	// lat / lon for current point
1147            	
1148            	double bx_tmp, by_tmp, bz_tmp;
1149            	
1150            	//
1151            	
1152            	for (int row = 0; row < mInfo->nrow; row++) {
1153            		for (int col = 0; col < mInfo->ncol; col++) {
1154            			
1155 xudong 1.1 			ind_map = row * mInfo->ncol + col;
1156            			
1157            			x = (col + 0.5 - ncol / 2.) * xscale;
1158            			y = (row + 0.5 - nrow / 2.) * yscale;
1159            			
1160            			if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1161            				bx_map[ind_map] = DRMS_MISSING_FLOAT;
1162            				by_map[ind_map] = DRMS_MISSING_FLOAT;
1163            				bz_map[ind_map] = DRMS_MISSING_FLOAT;
1164            				continue;
1165            			}
1166            			
1167            			bx_tmp = by_tmp = bz_tmp = 0;
1168            			
1169            			img2helioVector (bx_map[ind_map], by_map[ind_map], bz_map[ind_map],
1170            							 &bx_tmp, &by_tmp, &bz_tmp,
1171            							 lon, lat, disk_lonc, disk_latc, pa);
1172            			
1173            			bx_map[ind_map] = bx_tmp;
1174            			by_map[ind_map] = by_tmp;
1175            			bz_map[ind_map] = bz_tmp;
1176 xudong 1.1 			
1177            		}
1178            	}
1179            
1180            }
1181            
1182            
1183            
1184            /* 
1185             * Map and propogate vector field errors
1186             *
1187             */
1188            
1189            int getBErr(float *bx_err, float *by_err, float *bz_err,
1190            			 DRMS_Record_t *inRec, struct mapInfo *mInfo)
1191            {
1192            	
1193            	int status = 0;
1194            	
1195            	// Get variances and covariances, filling factor assume to be 1
1196            	
1197 xudong 1.1 	float *bT = (float *) (malloc(FOURK2 * sizeof(float)));	// field
1198            	float *bI = (float *) (malloc(FOURK2 * sizeof(float)));	// inclination
1199            	float *bA = (float *) (malloc(FOURK2 * sizeof(float)));	// azimuth
1200            	
1201            	float *errbT = (float *) (malloc(FOURK2 * sizeof(float)));
1202            	float *errbI = (float *) (malloc(FOURK2 * sizeof(float)));
1203            	float *errbA = (float *) (malloc(FOURK2 * sizeof(float)));
1204            	
1205            	float *errbTbI = (float *) (malloc(FOURK2 * sizeof(float)));
1206            	float *errbTbA = (float *) (malloc(FOURK2 * sizeof(float)));
1207            	float *errbIbA = (float *) (malloc(FOURK2 * sizeof(float)));
1208            	
1209            	if (readVectorBErr(inRec, 
1210            					   bT, bI, bA,
1211            					   errbT, errbI, errbA, 
1212            					   errbTbI, errbTbA, errbIbA)) {
1213            		printf("Read full disk variances & covariances error\n");
1214            		free(bT); free(bI); free(bA);
1215            		free(errbT); free(errbI); free(errbA);
1216            		free(errbTbI); free(errbTbA); free(errbIbA);
1217            		return 1;
1218 xudong 1.1 	}
1219            	
1220            	// Size
1221            	
1222            	int ncol = mInfo->ncol;
1223            	int nrow = mInfo->nrow;
1224            	
1225            	float xscale = mInfo->xscale * RADSINDEG;		// in rad
1226            	float yscale = mInfo->yscale * RADSINDEG;
1227            	
1228            	double lonc = mInfo->xc * RADSINDEG;	// in rad
1229            	double latc = mInfo->yc * RADSINDEG;
1230            	
1231            	double disk_lonc = (mInfo->ephem).disk_lonc;
1232            	double disk_latc = (mInfo->ephem).disk_latc;
1233            	
1234            	double rSun = (mInfo->ephem).rSun;
1235            	double disk_xc = (mInfo->ephem).disk_xc / rSun;
1236            	double disk_yc = (mInfo->ephem).disk_yc / rSun;
1237            	double pa = (mInfo->ephem).pa;
1238            	
1239 xudong 1.1 	// Start
1240            	
1241            	double x, y;          // map coord
1242            	double lat, lon;      // spherical coord
1243            	double xi, zeta;      // image coord, round to full pixel
1244            	
1245            	int ind_map, ind_img;
1246            	
1247            	double bpSigma2, btSigma2, brSigma2;		// variances after prop
1248            
1249            	for (int row = 0; row < mInfo->nrow; row++) {
1250            		for (int col = 0; col < mInfo->ncol; col++) {
1251            			
1252            			ind_map = row * mInfo->ncol + col;
1253            			
1254            			x = (col + 0.5 - ncol / 2.) * xscale;
1255            			y = (row + 0.5 - nrow / 2.) * yscale;
1256            			
1257            			if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1258            				bx_err[ind_map] = DRMS_MISSING_FLOAT;
1259            				by_err[ind_map] = DRMS_MISSING_FLOAT;
1260 xudong 1.1 				bz_err[ind_map] = DRMS_MISSING_FLOAT;
1261            				continue;
1262            			}
1263            			
1264            			if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta, 
1265            							disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1266            				bx_err[ind_map] = DRMS_MISSING_FLOAT;
1267            				bx_err[ind_map] = DRMS_MISSING_FLOAT;
1268            				bx_err[ind_map] = DRMS_MISSING_FLOAT;
1269            				continue;
1270            			}
1271            			
1272            			xi *= rSun; xi = round(xi);
1273            			zeta *= rSun; zeta = round(zeta);     // nearest neighbor
1274            			
1275            			ind_img = round(zeta * FOURK + xi);
1276            			
1277            			if (errorprop(bT, bA, bI, 
1278            						  errbT, errbA, errbI, errbTbA, errbTbI, errbIbA, 
1279            						  lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta, 
1280            						  &btSigma2, &bpSigma2, &brSigma2)) {
1281 xudong 1.1 				bx_err[ind_map] = DRMS_MISSING_FLOAT;
1282            				by_err[ind_map] = DRMS_MISSING_FLOAT;
1283            				bz_err[ind_map] = DRMS_MISSING_FLOAT;
1284            				continue;
1285            			}
1286            			
1287            			bx_err[ind_map] = sqrt(bpSigma2);
1288            			by_err[ind_map] = sqrt(btSigma2);
1289            			bz_err[ind_map] = sqrt(brSigma2);
1290            			
1291            		}
1292            	}
1293            	
1294            	//
1295            	
1296            	free(bT); free(bI); free(bA);
1297            	free(errbT); free(errbI); free(errbA);
1298            	free(errbTbI); free(errbTbA); free(errbIbA);
1299            	return 0;
1300            	
1301            }
1302 xudong 1.1 
1303            
1304            
1305            /*
1306             * Read full disk vector magnetograms
1307             * Fill factor is 1, use default disambiguity resolution
1308             *
1309             */
1310            
1311            int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img)
1312            {
1313            	
1314            	int status = 0;
1315            	
1316            	DRMS_Segment_t *inSeg;
1317            	DRMS_Array_t *inArray_ambig;
1318                    DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1319            	
1320            	char *ambig;
1321            	float *bTotal, *bAzim, *bIncl;
1322            	
1323 xudong 1.1 	inSeg = drms_segment_lookup(inRec, "disambig");
1324            	inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
1325            	if (status) return 1;
1326            	ambig = (char *)inArray_ambig->data;
1327            	
1328            	inSeg = drms_segment_lookup(inRec, "field");
1329            	inArray_bTotal = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1330            	if (status) return 1;
1331            	bTotal = (float *)inArray_bTotal->data;
1332            	
1333            	inSeg = drms_segment_lookup(inRec, "azimuth");
1334            	inArray_bAzim = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1335            	if (status) return 1;
1336            	bAzim = (float *)inArray_bAzim->data;
1337            	
1338            	inSeg = drms_segment_lookup(inRec, "inclination");
1339            	inArray_bIncl = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1340            	if (status) return 1;
1341            	bIncl = (float *)inArray_bIncl->data;
1342            	
1343            	// Convert CCD xyz
1344 xudong 1.1 	
1345            	int llx, lly;		// lower-left corner
1346            	int bmx, bmy;		// bitmap size
1347            	
1348            	llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1349            	lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1350            	
1351            	bmx = inArray_ambig->axis[0];
1352            	bmy = inArray_ambig->axis[1];
1353            	
1354            	int kx, ky, kOff;
1355            	int ix = 0, jy = 0, yOff = 0, iData = 0;
1356            	int xDim = FOURK, yDim = FOURK;
1357            	
1358            	for (jy = 0; jy < yDim; jy++)
1359            	{
1360            		ix = 0;
1361            		yOff = jy * xDim;
1362            		ky = jy - lly;
1363            		for (ix = 0; ix < xDim; ix++)
1364            		{
1365 xudong 1.1 			iData = yOff + ix;
1366            			kx = ix - llx;
1367            			
1368            			// zero azi pointing up, zero incl pointing out from sun
1369            			bx_img[iData] = - bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * sin(bAzim[iData] * RADSINDEG);
1370            			by_img[iData] = bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * cos(bAzim[iData] * RADSINDEG);
1371            			bz_img[iData] = bTotal[iData] * cos(bIncl[iData] * RADSINDEG);
1372                        
1373            			// Disambiguation
1374            			
1375            			if (kx < 0 || kx >= bmx || ky < 0 || ky >= bmy) {
1376            				continue;
1377            			} else {
1378            				kOff = ky * bmx + kx;
1379            				if (ambig[kOff] % 2) {		// 180
1380            					bx_img[iData] *= -1.; by_img[iData] *= -1.;
1381            				} 
1382            			}
1383            		}
1384            	}
1385            	
1386 xudong 1.1 	// Clean up
1387            	
1388            	drms_free_array(inArray_ambig);
1389            	drms_free_array(inArray_bTotal);
1390            	drms_free_array(inArray_bAzim);
1391            	drms_free_array(inArray_bIncl);
1392            	
1393            	return 0;
1394            	
1395            }
1396            
1397            
1398            /*
1399             * Read variances and covariances of vector magnetograms
1400             *
1401             */
1402            
1403            int readVectorBErr(DRMS_Record_t *inRec, 
1404            				   float *bT, float *bI, float *bA,
1405            				   float *errbT, float *errbI, float *errbA, 
1406            				   float *errbTbI, float *errbTbA, float *errbIbA)
1407 xudong 1.1 {
1408            	
1409            	int status = 0;
1410            	
1411            	float *data_ptr[9];
1412            	char *segName[9] = {"field", "inclination", "azimuth",
1413            						"field_err", "inclination_err", "azimuth_err",
1414            						"field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1415            	DRMS_Segment_t *inSegs[9];
1416            	DRMS_Array_t *inArrays[9];
1417            	
1418            	// Read full disk images
1419            	
1420            	for (int iSeg = 0; iSeg < 9; iSeg++) {
1421            		
1422            		inSegs[iSeg] = drms_segment_lookup(inRec, segName[iSeg]);
1423            		inArrays[iSeg] = drms_segment_read(inSegs[iSeg], DRMS_TYPE_FLOAT, &status);
1424            		data_ptr[iSeg] = (float *) inArrays[iSeg]->data;
1425            		
1426            	}
1427            	
1428 xudong 1.1 	float *bT0 = data_ptr[0], *bI0 = data_ptr[1], *bA0 = data_ptr[2];
1429            	float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];
1430            	float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];
1431            	
1432            	// Convert errors to variances, correlation coefficients to covariances
1433            	
1434            	for (int i = 0; i < FOURK2; i++) {
1435            		
1436            		if (fabs(errbI0[i]) > 180.) errbI0[i] = 180.;
1437            		if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;
1438            		
1439            		bT[i] = bT0[i];
1440            		bI[i] = bI0[i];
1441            		bA[i] = bA0[i];
1442            		
1443            		errbT[i] = errbT0[i] * errbT0[i];
1444            		errbI[i] = errbI0[i] * errbI0[i] * RADSINDEG * RADSINDEG;
1445            		errbA[i] = errbA0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1446            		
1447            		errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG;
1448                    errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG;
1449 xudong 1.1         errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1450            	
1451            	}
1452            	
1453            	//
1454            	
1455            	for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
1456            
1457            	return 0;
1458            	
1459            }
1460            
1461            
1462            /*
1463             * Create Cutout record: top level subroutine
1464             * Do the loops on segments and set the keywords here
1465             * Work is done in writeCutout routine below
1466             * 
1467             */
1468            
1469            int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec, 
1470 xudong 1.1 					DRMS_Record_t *dopRec, DRMS_Record_t *contRec, 
1471            					DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
1472            {
1473            	
1474            	int status = 0;
1475            	
1476            	int iHarpSeg;
1477            	int nMharpSegs = ARRLENGTH(MharpSegs), nBharpSegs = ARRLENGTH(BharpSegs);
1478            	
1479            	// Cutout Mharp
1480            	
1481            	for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) {
1482            		if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) {
1483            			printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]);
1484            			break;
1485            		}
1486            	}
1487 xudong 1.3 	if (iHarpSeg != nMharpSegs) {
1488            		SHOW("Cutout: segment number unmatch\n");
1489            		return 1;		// if failed
1490            	}
1491 xudong 1.1 	printf("Magnetogram cutout done.\n");
1492            	
1493            	// Cutout Doppler
1494            	
1495            	if (writeCutout(sharpRec, dopRec, mharpRec, "Dopplergram")) {
1496            		printf("Doppler cutout failed\n");
1497            		return 1;
1498            	}
1499            	printf("Dopplergram cutout done.\n");
1500            	
1501            	// Cutout Continuum
1502            	
1503            	if (writeCutout(sharpRec, contRec, mharpRec, "continuum")) {
1504            		printf("Continuum cutout failed\n");
1505            		return 1;
1506            	}
1507            	printf("Intensitygram cutout done.\n");
1508            	
1509            	// Coutout Bharp
1510            	
1511            	for (iHarpSeg = 0; iHarpSeg < nBharpSegs; iHarpSeg++) {
1512 xudong 1.1 		if (writeCutout(sharpRec, bharpRec, mharpRec, BharpSegs[iHarpSeg])) {
1513            			printf("Bharp cutout fails for %s\n", BharpSegs[iHarpSeg]);
1514            			break;
1515            		}
1516            	}
1517            	if (iHarpSeg != nBharpSegs) return 1;		// if failed
1518            	printf("Vector B cutout done.\n");
1519            	
1520            	// Keywords & Links
1521            	
1522            	drms_copykey(sharpRec, mharpRec, "T_REC");
1523            	drms_copykey(sharpRec, mharpRec, "HARPNUM");
1524            	
1525            	DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
1526            	if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
1527            	DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
1528            	if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
1529            	
1530            	setSWIndex(sharpRec, swKeys_ptr);	// Set space weather indices
1531            	setKeys(sharpRec, bharpRec);		// Set all other keywords
1532            
1533 xudong 1.1 	// Stats
1534            
1535            	int nCutSegs = ARRLENGTH(CutSegs);
1536            	for (int iSeg = 0; iSeg < nCutSegs; iSeg++) {
1537            		DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
1538            		DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
1539            		set_statistics(outSeg, outArray, 1);
1540            		drms_free_array(outArray);
1541            	}
1542            		
1543            	return 0;
1544            	
1545            }            
1546            
1547            
1548            /* 
1549             * Get cutout and write segment
1550             * Change DISAMB_AZI to apply disambiguation to azimuth
1551             *
1552             */
1553            
1554 xudong 1.1 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName)
1555            {
1556            	
1557            	int status = 0;
1558            	
1559            	DRMS_Segment_t *inSeg = NULL, *outSeg = NULL;
1560            	DRMS_Array_t *cutoutArray = NULL;
1561            	//	DRMS_Type_t arrayType;
1562            	
1563            	int ll[2], ur[2], nx, ny, nxny;		// lower-left and upper right coords
1564            	
1565            	/* Info */
1566            	
1567            	inSeg = drms_segment_lookup(inRec, SegName);
1568            	if (!inSeg) return 1;
1569            	
1570            	nx = (int) drms_getkey_float(harpRec, "CRSIZE1", &status);
1571            	ny = (int) drms_getkey_float(harpRec, "CRSIZE2", &status);
1572            	nxny = nx * ny;
1573            	ll[0] = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1; if (status) return 1;
1574            	ll[1] = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1; if (status) return 1;
1575 xudong 1.1 	ur[0] = ll[0] + nx - 1; if (status) return 1;
1576            	ur[1] = ll[1] + ny - 1; if (status) return 1;
1577            	
1578            	if (inSeg->axis[0] == nx && inSeg->axis[1] == ny) {			// for bitmaps, infomaps, etc.
1579            		cutoutArray = drms_segment_read(inSeg, DRMS_TYPE_DOUBLE, &status);
1580            		if (status) return 1;
1581            	} else if (inSeg->axis[0] == FOURK && inSeg->axis[1] == FOURK) {		// for full disk ones
1582            		cutoutArray = drms_segment_readslice(inSeg, DRMS_TYPE_DOUBLE, ll, ur, &status);
1583            		if (status) return 1;
1584            	} else {
1585            		return 1;
1586            	}
1587            	 
1588            	/* Adding disambiguation resolution to cutout azimuth? */
1589            
1590            #if DISAMB_AZI
1591            	if (!strcmp(SegName, "azimuth")) {
1592            		DRMS_Segment_t *disambSeg = drms_segment_lookup(inRec, "disambig");
1593            		if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
1594            		DRMS_Array_t *disambArray;
1595            		if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1596 xudong 1.1 			disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1597            			if (status) {drms_free_array(cutoutArray); return 1;}
1598            		} else {
1599            			return 1;
1600            		}
1601            		double *azimuth = (double *) cutoutArray->data;
1602            		char *disamb = (char *) disambArray->data;
1603            		for (int n = 0; n < nxny; n++) {
1604            			if (disamb[n]) azimuth[n] += 180.;
1605            		}
1606            		drms_free_array(disambArray);
1607            	}
1608            #endif
1609            
1610            	/* Write out */
1611            	
1612            	outSeg = drms_segment_lookup(outRec, SegName);
1613            	if (!outSeg) return 1;
1614            	drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray);
1615            	outSeg->axis[0] = cutoutArray->axis[0];
1616            	outSeg->axis[1] = cutoutArray->axis[1];
1617 xudong 1.1 	cutoutArray->parent_segment = outSeg;
1618            	cutoutArray->israw = 0;		// always compressed
1619                cutoutArray->bzero = outSeg->bzero;
1620                cutoutArray->bscale = outSeg->bscale;		// Same as inArray's
1621            	status = drms_segment_write(outSeg, cutoutArray, 0);
1622            	drms_free_array(cutoutArray);
1623            	if (status) return 1;
1624            	
1625            	return 0;
1626            	
1627            }
1628            
1629            
1630            /* 
1631             * Compute space weather indices, no error checking for now
1632             * Based on M. Bobra's swharp_vectorB.c
1633             * No error checking for now
1634             *
1635             */
1636            
1637            void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo)
1638 xudong 1.1 {
1639            	
1640            	int status = 0;
1641            	int nx = mInfo->ncol, ny = mInfo->nrow;
1642            	int nxny = nx * ny;
1643            	int dims[2] = {nx, ny};
1644            	// Get bx, by, bz, mask
1645            	
1646                    // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1647            	//DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "bitmap");
1648            	//DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1649            	//int *mask = (int *) maskArray->data;		// get the previously made mask array
1650            
1651                    //Use conf_disambig map as a threshold on spaceweather quantities 
1652 xudong 1.2 	DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");         
1653 xudong 1.1 	DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1654            	int *mask = (int *) maskArray->data;		// get the previously made mask array
1655            
1656            	DRMS_Segment_t *bxSeg = drms_segment_lookup(inRec, BP_SEG_CEA);
1657            	DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
1658            	float *bx = (float *) bxArray->data;		// bx
1659            	
1660            	DRMS_Segment_t *bySeg = drms_segment_lookup(inRec, BT_SEG_CEA);
1661            	DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);
1662            	float *by = (float *) byArray->data;		// by
1663            	for (int i = 0; i < nxny; i++) by[i] *= -1;
1664            	
1665            	DRMS_Segment_t *bzSeg = drms_segment_lookup(inRec, BR_SEG_CEA);
1666            	DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
1667            	float *bz = (float *) bzArray->data;		// bz
1668            	
1669            	// Get emphemeris
1670            	
1671            	//float cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
1672            	float cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);
1673            	float dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
1674 xudong 1.1 	double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);
1675            	double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);
1676            	float imcrpix1    = drms_getkey_float(inRec, "IMCRPIX1", &status);
1677            	float imcrpix2    = drms_getkey_float(inRec, "IMCRPIX2", &status);
1678            	float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);
1679            	float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);
1680            	
1681            	//float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)
1682            
1683                    printf("cdelt1=%f\n",cdelt1);
1684                    printf("rsun_ref=%f\n",rsun_ref);
1685                    printf("rsun_obs=%f\n",rsun_obs);
1686                    printf("dsun_obs=%f\n",dsun_obs);
1687            
1688 xudong 1.2 	// Temp arrays                
1689 xudong 1.1 	
1690            	float *bh = (float *) (malloc(nxny * sizeof(float)));
1691            	float *bt = (float *) (malloc(nxny * sizeof(float)));
1692            	float *jz = (float *) (malloc(nxny * sizeof(float)));
1693            	float *bpx = (float *) (malloc(nxny * sizeof(float)));
1694            	float *bpy = (float *) (malloc(nxny * sizeof(float)));
1695            	float *bpz = (float *) (malloc(nxny * sizeof(float)));
1696            	float *derx = (float *) (malloc(nxny * sizeof(float)));
1697            	float *dery = (float *) (malloc(nxny * sizeof(float)));
1698            	float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
1699            	float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
1700            	float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
1701            	float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
1702            	float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1703            	float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1704            	                                   
1705            	// Compute      
1706            	
1707            	if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf), 
1708            					   mask, cdelt1, rsun_ref, rsun_obs)){
1709            		swKeys_ptr->absFlux = DRMS_MISSING_FLOAT;		// If fail, fill in NaN
1710 xudong 1.1 		swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1711            	}
1712            	
1713            	for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
1714                   	greenpot(bpx, bpy, bpz, nx, ny);                      
1715            	
1716            	computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask);
1717            	
1718            	if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask))
1719            		swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
1720            	
1721            	computeB_total(bx, by, bz, bt, dims, mask);
1722            	
1723            	if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt))
1724            		swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
1725            	
1726            	if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh))
1727            		swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
1728            	
1729            	if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, derx_bz, dery_bz))
1730            		swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1731 xudong 1.1 	
1732            
1733            
1734            	if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, 
1735                                 cdelt1, rsun_ref, rsun_obs, derx, dery)) {
1736            		swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
1737            		swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
1738            	}
1739            	
1740                                    printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
1741            
1742                    if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, cdelt1, rsun_ref, rsun_obs))
1743            		swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
1744            	
1745            	if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih), 
1746            						&(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih), 
1747            						mask, cdelt1, rsun_ref, rsun_obs)) {  
1748            		swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT; 
1749            		swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
1750            		swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
1751            	}
1752 xudong 1.1 	
1753            	if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz), 
1754            								 mask, cdelt1, rsun_ref, rsun_obs))
1755            		swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
1756            
1757            	
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;
1763            	}
1764            	
1765            	if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims, 
1766            						  &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45), 
1767            						  &(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h), 
1768            						  mask)) {
1769            		swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1770            		swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
1771            		swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1772            		swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;
1773 xudong 1.1 	}
1774             
1775            	// Clean up
1776            	
1777            	drms_free_array(maskArray);
1778            	drms_free_array(bxArray);
1779            	drms_free_array(byArray);
1780            	drms_free_array(bzArray);
1781            	
1782            	free(bh); free(bt); free(jz);
1783            	free(bpx); free(bpy); free(bpz);
1784            	free(derx); free(dery);	
1785            	free(derx_bt); free(dery_bt);	
1786            	free(derx_bz); free(dery_bz);	
1787            	free(derx_bh); free(dery_bh);
1788            	
1789            }
1790            
1791            
1792            /* 
1793             * Set space weather indices, no error checking for now
1794 xudong 1.1  *
1795             */
1796            
1797            void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
1798            {
1799            	drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
1800            	drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
1801            	drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
1802            	drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
1803            	drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
1804            	drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
1805            	drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
1806            	drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
1807            	drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
1808            	drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
1809            	drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
1810            	drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
1811            	drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
1812            	drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);
1813            	drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
1814            	drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
1815 xudong 1.1 };
1816            
1817            
1818            /* 
1819             * Set all keywords, no error checking for now
1820             *
1821             */
1822            
1823            void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)
1824            {
1825               copy_me_keys(inRec, outRec);
1826               copy_patch_keys(inRec, outRec);
1827               copy_geo_keys(inRec, outRec);
1828               copy_ambig_keys(inRec, outRec);
1829 xudong 1.3 
1830               char timebuf[1024];
1831               float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1832               double val;
1833               int status = DRMS_SUCCESS;
1834            
1835               val = drms_getkey_double(inRec, "DATE",&status); 
1836               drms_setkey_double(outRec, "DATE_B", val);
1837               sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
1838               drms_setkey_string(outRec, "DATE", timebuf);
1839            
1840               // set cvs commit version into keyword HEADER
1841               char *cvsinfo = strdup("$Header: $");
1842               status = drms_setkey_string(outRec, "HEADER", cvsinfo);
1843            
1844 xudong 1.1 };
1845            
1846            //
1847            //
1848            
1849            /* ############# Nearest neighbour interpolation ############### */
1850            
1851            float nnb (float *f, int nx, int ny, double x, double y)
1852            {
1853            	
1854            	if (x <= -0.5 || y <= -0.5 || x > nx - 0.5 || y > ny - 0.5)
1855            		return DRMS_MISSING_FLOAT;
1856            	int ilow = floor (x);
1857            	int jlow = floor (y);
1858            	int i = ((x - ilow) > 0.5) ? ilow + 1 : ilow;
1859            	int j = ((y - jlow) > 0.5) ? jlow + 1 : jlow;
1860            	return f[j * nx + i];
1861            	
1862            }
1863            
1864            /* ################## Wrapper for Jesper's rebin code ################## */
1865 xudong 1.1 
1866            void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)
1867            {
1868            	
1869            	struct fresize_struct fresizes;
1870            	int nxout, nyout, xoff, yoff;
1871            	int nlead = nx;
1872            	
1873            	nxout = nx / nbin; nyout = ny / nbin;
1874            	if (gauss && nbin != 1)
1875            		init_fresize_gaussian(&fresizes, (nbin / 2), (nbin / 2 * 2), nbin);		// for nbin=3, sigma=1, half truncate width=2
1876            	else
1877            		init_fresize_bin(&fresizes, nbin);
1878            	xoff = nbin / 2 + nbin / 2;
1879            	yoff = nbin / 2 + nbin / 2;
1880            	fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);
1881            	
1882            }

Karen Tian
Powered by
ViewCVS 0.9.4