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

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

Karen Tian
Powered by
ViewCVS 0.9.4