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

Karen Tian
Powered by
ViewCVS 0.9.4