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

Karen Tian
Powered by
ViewCVS 0.9.4