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

Karen Tian
Powered by
ViewCVS 0.9.4