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

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

Karen Tian
Powered by
ViewCVS 0.9.4