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

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

Karen Tian
Powered by
ViewCVS 0.9.4