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

Karen Tian
Powered by
ViewCVS 0.9.4