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