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

Karen Tian
Powered by
ViewCVS 0.9.4