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

Karen Tian
Powered by
ViewCVS 0.9.4