107 mbobra 1.1 // Mapping method
108 enum projection {
109 carree,
110 cassini,
111 mercator,
112 cyleqa,
113 sineqa,
114 gnomonic,
115 postel,
116 stereographic,
117 orthographic,
118 lambert
119 };
120
121 // WSC code
122 char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
123 "SIN", "ZEA"};
124
125 // Ephemeris information
126 struct ephemeris {
127 double disk_lonc, disk_latc;
128 mbobra 1.1 double disk_xc, disk_yc;
129 double rSun, asd, pa;
130 };
131
132 // Mapping information
133 struct mapInfo {
134 float xc, yc; // reference point: center
135 int nrow, ncol; // size
136 float xscale, yscale; // scale
137 int nbin;
138 enum projection proj; // projection method
139 struct ephemeris ephem; // ephemeris info
140 float *xi_out, *zeta_out; // coordinate on full disk image to sample at
141 };
142
143 /* ========================================================================================================== */
144
145 /* Get all input data series */
146 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, char *mharpQuery);
147
148 /* Get other data series */
149 mbobra 1.1 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS);
150
151 /* Find record from record set with given T_rec */
152 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
153
154 /* Create CEA record */
155 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
156
157 /* Mapping single segment, wrapper */
158 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
159 struct mapInfo *mInfo, char *segName);
160
161 /* Determine reference point coordinate and patch size according to input */
162 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo);
163
164 /* Get ephemeris information */
165 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem);
166
167 /* Compute the coordinates at which the full disk image is sampled */
168 void findCoord(struct mapInfo *mInfo);
169
170 mbobra 1.1 /* Mapping function */
171 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
172
173 // ===================
174
175 /* Create Cutout record */
176 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
177
178 /* Get cutout and write segment */
179 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName);
180
181 // ===================
182
183 /* Compute space weather indices */
184 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo);
185
186 /* Set space weather indices */
187 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
188
189 /* Set all keywords */
190 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, struct mapInfo *mInfo);
191 mbobra 1.1
192 // ===================
193
194 /* Nearest neighbor interpolation */
195 float nnb (float *f, int nx, int ny, double x, double y);
196
197 /* Wrapper for Jesper's rebin code */
198 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss);
199
200 /* ========================================================================================================== */
201
202 /* Cutout segment names, input identical to output */
203 char *MharpSegs[] = {"magnetogram", "bitmap"};
204 char *CutSegs[] = {"magnetogram", "bitmap", "continuum"};
205 char *CEASegs[] = {"magnetogram", "bitmap", "continuum"};
206 // For BUNIT
207 char *CutBunits[] = {"Mx/cm^2", " ", "DN/s"};
208 char *CEABunits[] = {"Mx/cm^2", " ", "DN/s"};
209 /* ========================================================================================================== */
210
211 char *module_name = "smarp";
212 mbobra 1.1
213 ModuleArgs_t module_args[] =
214 {
215 {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
216 {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
217 {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
218 {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
219 {ARG_END}
220 };
221
222 int DoIt(void)
223 {
224 int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
225 int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
226 int status = DRMS_SUCCESS;
227 int nrecs, irec;
228 char *mharpQuery;
229 char *contQuery;
230 char *sharpCeaQuery, *sharpCutQuery;
231 DRMS_RecordSet_t *mharpRS = NULL;
232 DRMS_RecordSet_t *contRS = NULL;
233 mbobra 1.1
234 /* Get parameters */
235
236 mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
237 sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
238 sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
239 contQuery = (char *) params_get_str(&cmdparams, "cont");
240
241 /* Get input data, check everything */
242 if (getInputRS(&mharpRS, mharpQuery))
243 DIE("Input harp data error.");
244 nrecs = mharpRS->n;
245
246 if (getInputRS_aux(&contRS, contQuery, mharpRS))
247 DIE("Input continuum data error.");
248
249 /* Start */
250
251 printf("==============\nStart. %d image(s) in total.\n", nrecs);
252
253 for (irec = 0; irec < nrecs; irec++) {
254 mbobra 1.1
255 /* Records in work */
256
257 DRMS_Record_t *mharpRec = NULL;
258 DRMS_Record_t *contRec = NULL;
259
260 mharpRec = mharpRS->records[irec];
261 TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
262
263 struct swIndex swKeys;
264
265 if (getInputRec_aux(&contRec, contRS, trec)) {
266 printf("Fetching Continuum failed, image #%d skipped.\n", irec);
267 continue;
268 }
269
270 printf("Obtained all the data \n");
271
272 /* Create CEA record */
273
274 DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
275 mbobra 1.1 if (status) { // if failed
276 printf("Creating CEA failed, image #%d skipped.\n", irec);
277 continue;
278 }
279 if (createCeaRecord(mharpRec, contRec, sharpCeaRec, &swKeys)) { // do the work
280 printf("Creating CEA failed, image #%d skipped.\n", irec);
281 drms_close_record(sharpCeaRec, DRMS_FREE_RECORD);
282 continue;
283 } // swKeys updated here
284
285 drms_close_record(sharpCeaRec, DRMS_INSERT_RECORD);
286
287 printf("Created CEA record \n");
288
289 /* Create Cutout record */
290
291 DRMS_Record_t *sharpCutRec = drms_create_record(drms_env, sharpCutQuery, DRMS_PERMANENT, &status);
292 if (status) { // if failed
293 printf("Creating cutout failed, image #%d skipped.\n", irec);
294 continue;
295 }
296 mbobra 1.1
297 if (createCutRecord(mharpRec, contRec, sharpCutRec, &swKeys)) { // do the work
298 printf("Creating cutout failed, image #%d skipped.\n", irec);
299 drms_close_record(sharpCutRec, DRMS_FREE_RECORD);
300 continue;
301 } // swKeys used here
302 drms_close_record(sharpCutRec, DRMS_INSERT_RECORD);
303 printf("Created CUT record \n");
304 /* Done */
305
306 printf("Image #%d done.\n", irec);
307
308 } // irec
309
310
311 drms_close_records(mharpRS, DRMS_FREE_RECORD);
312 drms_close_records(contRS, DRMS_FREE_RECORD);
313
314 return 0;
315
316 } // DoIt
317 mbobra 1.1
318 // ===================================================================
319 // ===================================================================
320 // ===================================================================
321
322 /*
323 * Get input data series, including mHarp and bharp
324 * Need all records to match, otherwise quit
325 *
326 */
327
328 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, char *mharpQuery)
329 {
330 int status = 0;
331 *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
332 if (status || (*mharpRS_ptr)->n == 0) return 1;
333 return 0;
334 }
335
336 /*
337 * Get other data series, check all T_REC are available
338 mbobra 1.1 */
339
340 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS)
341 {
342
343 int status = 0;
344
345 *inRS_ptr = drms_open_records(drms_env, inQuery, &status);
346 if (status || (*inRS_ptr)->n == 0) return status;
347
348 // Check if all T_rec are available, need to match both ways
349 int n = harpRS->n, n0 = (*inRS_ptr)->n;
350
351 for (int i0 = 0; i0 < n0; i0++) {
352 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
353 TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
354 TIME trec = 0;
355 for (int i = 0; i < n; i++) {
356 DRMS_Record_t *harpRec = harpRS->records[i];
357 trec = drms_getkey_time(harpRec, "T_REC", &status);
358 if (fabs(trec0 - trec) < 10) break;
359 mbobra 1.1 }
360 if (fabs(trec0 - trec) >= 10) return 1;
361 }
362
363 for (int i = 0; i < n; i++) {
364 DRMS_Record_t *harpRec = harpRS->records[i];
365 TIME trec = drms_getkey_time(harpRec, "T_REC", &status);
366 TIME trec0 = 0;
367 for (int i0 = 0; i0 < n0; i0++) {
368 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
369 trec0 = drms_getkey_time(inRec, "T_REC", &status);
370 if (fabs(trec0 - trec) < 10) break;
371 }
372 if (fabs(trec0 - trec) >= 10) return 1;
373 }
374
375 return 0;
376
377 }
378
379 /*
380 mbobra 1.1 * Find record from record set with given T_rec
381 */
382
383 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec)
384 {
385
386 int status = 0;
387
388 int n = inRS->n;
389 for (int i = 0; i < n; i++) {
390 *inRec_ptr = inRS->records[i];
391 TIME trec0 = drms_getkey_time((*inRec_ptr), "T_REC", &status);
392 if (fabs(trec0 - trec) < 10) return 0;
393 }
394
395 return 1;
396
397 }
398
399
400
401 mbobra 1.1
402 /*
403 * Create CEA record: top level subroutine
404 * Also compute all the space weather keywords here
405 */
406
407 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
408 {
409
410 int status = 0;
411 DRMS_Segment_t *inSeg;
412 DRMS_Array_t *inArray;
413 int val;
414
415 struct mapInfo mInfo;
416 mInfo.proj = (enum projection) cyleqa; // projection method
417 mInfo.xscale = XSCALE;
418 mInfo.yscale = YSCALE;
419
420 int ncol0, nrow0; // oversampled map size
421
422 mbobra 1.1 // Get ephemeris
423
424 if (getEphemeris(mharpRec, &(mInfo.ephem))) {
425 SHOW("CEA: get ephemeris error\n");
426 return 1;
427 }
428
429 // Find position
430
431 if (findPosition(mharpRec, &mInfo)) {
432 SHOW("CEA: find position error\n");
433 return 1;
434 }
435
436 // ========================================
437 // Do this for all bitmaps, Aug 12 2013 XS
438 // ========================================
439
440 mInfo.nbin = 1; // for bitmaps. suppress anti-aliasing
441 ncol0 = mInfo.ncol;
442 nrow0 = mInfo.nrow;
443 mbobra 1.1
444 mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
445 mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
446
447 findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions
448
449 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
450 SHOW("CEA: mapping bitmap error\n");
451 return 1;
452 }
453 printf("Bitmap mapping done.\n");
454
455 free(mInfo.xi_out);
456 free(mInfo.zeta_out);
457
458 // ========================================
459 // Do this again for floats, Aug 12 2013 XS
460 // ========================================
461 // Create xi_out, zeta_out array in mInfo:
462 // Coordinates to sample in original full disk image
463
464 mbobra 1.1 mInfo.nbin = NBIN;
465 ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
466 nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
467
468 mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
469 mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
470
471 findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions
472
473 // Mapping single segment: Mharp, etc.
474
475 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
476 SHOW("CEA: mapping magnetogram error\n");
477 return 1;
478 }
479 printf("Magnetogram mapping done.\n");
480
481 if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
482 SHOW("CEA: mapping continuum error\n");
483 return 1;
484 }
485 mbobra 1.1 printf("Intensitygram mapping done.\n");
486
487 // Keywords & Links
488 copy_patch_keys(mharpRec, sharpRec);
489 copy_geo_keys(mharpRec, sharpRec);
490 // rename HARPNUM to TARPNUM
491 val = drms_getkey_double(mharpRec, "HARPNUM", &status);
492 drms_setkey_double(sharpRec, "TARPNUM", val);
493 // copy everything else
494 drms_copykey(sharpRec, mharpRec, "T_REC");
495 drms_copykey(sharpRec, mharpRec, "CDELT1");
496 drms_copykey(sharpRec, mharpRec, "RSUN_OBS");
497 drms_copykey(sharpRec, mharpRec, "DSUN_OBS");
498 drms_copykey(sharpRec, mharpRec, "OBS_VR");
499 drms_copykey(sharpRec, mharpRec, "OBS_VW");
500 drms_copykey(sharpRec, mharpRec, "OBS_VN");
501 drms_copykey(sharpRec, mharpRec, "CRLN_OBS");
502 drms_copykey(sharpRec, mharpRec, "CRLT_OBS");
503 drms_copykey(sharpRec, mharpRec, "CAR_ROT");
504 drms_copykey(sharpRec, mharpRec, "SIZE_SPT");
505 drms_copykey(sharpRec, mharpRec, "AREA_SPT");
506 mbobra 1.1 drms_copykey(sharpRec, mharpRec, "DATE__OBS");
507 drms_copykey(sharpRec, mharpRec, "T_OBS");
508 drms_copykey(sharpRec, mharpRec, "T_MAXPIX");
509 drms_copykey(sharpRec, mharpRec, "QUALITY");
510 drms_copykey(sharpRec, mharpRec, "NPIX_SPT");
511 drms_copykey(sharpRec, mharpRec, "ARS_NCLN");
512 drms_copykey(sharpRec, mharpRec, "ARS_MODL");
513 drms_copykey(sharpRec, mharpRec, "ARS_EDGE");
514 drms_copykey(sharpRec, mharpRec, "ARS_BETA");
515 drms_copykey(sharpRec, mharpRec, "T_MID1");
516 drms_copykey(sharpRec, mharpRec, "T_CMPASS");
517
518 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MTARP");
519 if (mHarpLink) {
520 drms_link_set("MTARP", sharpRec, mharpRec);
521 }
522
523 // set other keywords
524 setKeys(sharpRec, mharpRec, &mInfo);
525
526 // set space weather keywords
527 mbobra 1.1 computeSWIndex(swKeys_ptr, sharpRec, &mInfo);
528 printf("Space weather indices done.\n");
529 setSWIndex(sharpRec, swKeys_ptr);
530
531 // set statistical keywords (e.g. DATAMIN, DATAMAX, etc.)
532 //int nCEASegs = ARRLENGTH(CEASegs);
533 int nCEASegs = 3;
534 for (int iSeg = 0; iSeg < 3; iSeg++) {
535 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
536 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
537 int stat = set_statistics(outSeg, outArray, 1);
538 //printf("%d => %d\n", iSeg, stat);
539 drms_free_array(outArray);
540 }
541
542 free(mInfo.xi_out);
543 free(mInfo.zeta_out);
544 return 0;
545
546 }
547
548 mbobra 1.1 /*
549 * Mapping a single segment
550 * Read in full disk image, utilize mapImage for mapping
551 * then write the segment out, segName same in in/out Rec
552 */
553
554 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
555 struct mapInfo *mInfo, char *segName)
556 {
557
558 int status = 0;
559 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
560 int dims[2] = {nx, ny};
561 int interpOpt = INTERP; // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
562
563 // Input full disk array
564
565 DRMS_Segment_t *inSeg = NULL;
566 inSeg = drms_segment_lookup(inRec, segName);
567 if (!inSeg) return 1;
568
569 mbobra 1.1 DRMS_Array_t *inArray = NULL;
570 inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
571 if (!inArray) return 1;
|
578 mbobra 1.1
579 float *inData;
580 int xsz = inArray->axis[0], ysz = inArray->axis[1];
581 if ((xsz != FOURK) || (ysz != FOURK)) { // for bitmap, make tmp full disk
582 float *inData0 = (float *) inArray->data;
583 inData = (float *) (calloc(FOURK2, sizeof(float)));
584 int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
585 int y0 = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1;
586 int ind_map;
587 for (int row = 0; row < ysz; row++) {
588 for (int col = 0; col < xsz; col++) {
589 ind_map = (row + y0) * FOURK + (col + x0);
590 inData[ind_map] = inData0[row * xsz + col];
591 }
592 }
593 drms_free_array(inArray); inArray = NULL;
594 } else {
595 inData = (float *) inArray->data;
596 }
597
598 // Mapping
599 mbobra 1.1
600 float *map = (float *) (malloc(nxny * sizeof(float)));
601 if (performSampling(map, inData, mInfo, interpOpt)) // Add interpOpt for different types, Aug 12 XS
602 {if (inArray) drms_free_array(inArray); free(map); return 1;}
603
604 // Write out
605
606 DRMS_Segment_t *outSeg = NULL;
607 outSeg = drms_segment_lookup(sharpRec, segName);
608 if (!outSeg) return 1;
609
610 // DRMS_Type_t arrayType = outSeg->info->type;
611 DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
612 if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
613
614 // convert to needed data type
615
616 // drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray); // Jan 02 2013
617
618 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
619 // outArray->parent_segment = outSeg;
620 mbobra 1.1 outArray->israw = 0; // always compressed
621 outArray->bzero = outSeg->bzero;
622 outArray->bscale = outSeg->bscale;
623
624 status = drms_segment_write(outSeg, outArray, 0);
625 if (status) return 0;
626
627 if (inArray) drms_free_array(inArray);
628 if ((xsz != FOURK) || (ysz != FOURK)) free(inData); // Dec 18 2012
629 if (outArray) drms_free_array(outArray);
630 return 0;
631
632 }
633
634 /*
635 * Determine reference point coordinate and patch size according to keywords
636 * xc, yc are the coordinate of patch center, in degrees
637 * ncol and nrow are the final size
638 */
639
640 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo)
641 mbobra 1.1 {
642
643 int status = 0;
644 int harpnum = drms_getkey_int(inRec, "TARPNUM", &status);
645 TIME trec = drms_getkey_time(inRec, "T_REC", &status);
646 float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
647
648 /* Center coord */
649 // Changed into double Jun 16 2014 XS
650
651 double minlon = drms_getkey_double(inRec, "LONDTMIN", &status); if (status) return 1; // Stonyhurst lon
652 double maxlon = drms_getkey_double(inRec, "LONDTMAX", &status); if (status) return 1;
653 double minlat = drms_getkey_double(inRec, "LATDTMIN", &status); if (status) return 1;
654 double maxlat = drms_getkey_double(inRec, "LATDTMAX", &status); if (status) return 1;
655
656 // A bug fixer for HARP (per M. Turmon)
657 // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
658 // Also keywords such as "SIZE" will be NaN
659 // We compute minlon & minlat then by
660 // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
661
662 mbobra 1.1 // float psize = drms_getkey_float(inRec, "SIZE", &status);
663 // if (psize != psize) {
664
665 if (minlon != minlon || maxlon != maxlon) { // check lons instead of SIZE
666 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
667 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
668 char firstRecQuery[100], t0_str[100];
669 sprint_time(t0_str, t0, "TAI", 0);
670 snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str);
671 DRMS_RecordSet_t *tmpRS = drms_open_records(drms_env, firstRecQuery, &status);
672 if (status || tmpRS->n != 1) return 1;
673 DRMS_Record_t *tmpRec = tmpRS->records[0];
674 double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1;
675 double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1;
676 minlon = minlon0 + (trec - t0) * omega / SECINDAY;
677 maxlon = maxlon0 + (trec - t0) * omega / SECINDAY;
678 printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon);
679 }
680
681 mInfo->xc = (maxlon + minlon) / 2. + disk_lonc;
682 mInfo->yc = (maxlat + minlat) / 2.;
683 mbobra 1.1
684 /* Size */
685 // Rounded to 1.d3 precision first. Jun 16 2014 XS
686 // The previous fix does not work. LONDTMAX-LONDTMIN varies from frame to frame
687 // Need to find out the maximum possible difference, MAXLONDIFF (1.2e-4)
688 // Now, ncol = (maxlon-minlon)/xscale, if the decimal part is outside 0.5 \pm (MAXLONDIFF/xscale)
689 // proceed as it is. else, all use floor on ncol
690
691 float dpix = (MAXLONDIFF / mInfo->xscale) * 1.5; // "danger zone"
692 float ncol = (maxlon - minlon) / mInfo->xscale;
693 float d_ncol = fabs(ncol - floor(ncol) - 0.5); // distance to 0.5
694 if (d_ncol < dpix) {
695 mInfo->ncol = floor(ncol);
696 } else {
697 mInfo->ncol = round(ncol);
698 }
699
700 mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
701
702 return 0;
703
704 mbobra 1.1 }
705
706
707 /*
708 * Fetch ephemeris info from a DRMS record
709 */
710
711 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem)
712 {
713
714 int status = 0;
715
716 float crota2 = drms_getkey_float(inRec, "CROTA2", &status); // rotation
717 double sina = sin(crota2 * RADSINDEG);
718 double cosa = cos(crota2 * RADSINDEG);
719
720 ephem->pa = - crota2 * RADSINDEG;
721 ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG;
722 ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG;
723
724 float crvalx = 0.0;
725 mbobra 1.1 float crvaly = 0.0;
726 float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
727 float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
728 float cdelt = drms_getkey_float(inRec, "CDELT1", &status); // in arcsec, assumimg dx=dy
729 ephem->disk_xc = PIX_X(0.0,0.0) - 1.0; // Center of disk in pixel, starting at 0
730 ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
731
732 float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status);
733 float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status);
734 if (status) rSun_ref = 6.96e8;
735
736 ephem->asd = asin(rSun_ref/dSun);
737 ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt;
738
739 return 0;
740
741 }
742
743
744 /*
745 * Compute the coordinates to be sampled on full disk image
746 mbobra 1.1 * mInfo->xi_out & mInfo->zeta_out
747 * This is oversampled, its size is ncol0 & nrow0 as shown below
748 */
749
750 void findCoord(struct mapInfo *mInfo)
751 {
752
753 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
754 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
755
756 float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG; // oversampling resolution
757 float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG; // in rad
758
759 double lonc = mInfo->xc * RADSINDEG; // in rad
760 double latc = mInfo->yc * RADSINDEG;
761
762 double disk_lonc = (mInfo->ephem).disk_lonc;
763 double disk_latc = (mInfo->ephem).disk_latc;
764
765 double rSun = (mInfo->ephem).rSun;
766 double disk_xc = (mInfo->ephem).disk_xc / rSun;
767 mbobra 1.1 double disk_yc = (mInfo->ephem).disk_yc / rSun;
768 double pa = (mInfo->ephem).pa;
769
770 // Temp pointers
771
772 float *xi_out = mInfo->xi_out;
773 float *zeta_out = mInfo->zeta_out;
774
775 // start
776
777 double x, y; // map coord
778 double lat, lon; // helio coord
779 double xi, zeta; // image coord (for one point)
780
781 int ind_map;
782
783 for (int row0 = 0; row0 < nrow0; row0++) {
784 for (int col0 = 0; col0 < ncol0; col0++) {
785
786 ind_map = row0 * ncol0 + col0;
787
788 mbobra 1.1 x = (col0 + 0.5 - ncol0/2.) * xscale0; // in rad
789 y = (row0 + 0.5 - nrow0/2.) * yscale0;
790
791 /* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
792 * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
793 * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
794 * is the heliographic longitude and latitude of the map center. Both are in degree.
795 */
796
797 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
798 xi_out[ind_map] = -1;
799 zeta_out[ind_map] = -1;
800 continue;
801 }
802
803 /* map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the
804 * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
805 */
806
807 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
808 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
809 mbobra 1.1 xi_out[ind_map] = -1;
810 zeta_out[ind_map] = -1;
811 continue;
812 }
813
814 xi_out[ind_map] = xi * rSun;
815 zeta_out[ind_map] = zeta * rSun;
816
817 }
818 }
819
820 }
821
822
823 /*
824 * Sampling function
825 * oversampling by nbin, then binning using a Gaussian
826 * save results in outData, always of float type
827 */
828
829 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
830 mbobra 1.1 {
831
832 int status = 0;
833 int ind_map;
834
835 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
836 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
837
838 // Changed Aug 12 2013, XS, for bitmaps
839 float *outData0;
840 if (interpOpt == 3 && mInfo->nbin == 1) {
841 outData0 = outData;
842 } else {
843 outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
844 }
845
846 float *xi_out = mInfo->xi_out;
847 float *zeta_out = mInfo->zeta_out;
848
849 // Interpolation
850
851 mbobra 1.1 struct fint_struct pars;
852 // Aug 12 2013, passed in as argument now
853
854 switch (interpOpt) {
855 case 0: // Wiener, 6 order, 1 constraint
856 init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath);
857 break;
858 case 1: // Cubic convolution
859 init_finterpolate_cubic_conv(&pars, 1., 3.);
860 break;
861 case 2: // Bilinear
862 init_finterpolate_linear(&pars, 1.);
863 break;
864 case 3: // Near neighbor
865 break;
866 default:
867 return 1;
868 }
869
870 //printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
871 if (interpOpt == 3) { // Aug 6 2013, Xudong
872 mbobra 1.1 for (int row0 = 0; row0 < nrow0; row0++) {
873 for (int col0 = 0; col0 < ncol0; col0++) {
874 ind_map = row0 * ncol0 + col0;
875 outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
876 }
877 }
878 } else {
879 finterpolate(&pars, inData, xi_out, zeta_out, outData0,
880 FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
881 }
882
883 // Rebinning, smoothing
884
885 if (interpOpt == 3 && mInfo->nbin == 1) {
886 return 0;
887 } else {
888 frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian
889 free(outData0);
890 }
891
892 return 0;
893 mbobra 1.1
894 }
895
896 /*
897 * Create Cutout record: top level subroutine
898 * Do the loops on segments and set the keywords here
899 * Work is done in writeCutout routine below
900 */
901
902 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
903 {
904
905 int status = 0;
906 int val;
907 int iHarpSeg;
908 int nMharpSegs = ARRLENGTH(MharpSegs);
909
910 // Cutout Mharp
911
912 for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) {
913 if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) {
914 mbobra 1.1 printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]);
915 printf("iHarpSeg nMharpSegs %d %d \n",iHarpSeg,nMharpSegs);
916 break;
917 }
918 }
919 if (iHarpSeg != nMharpSegs) {
920 SHOW("Cutout: segment number mismatch\n");
921 return 1; // if failed
922 }
923 printf("Magnetogram cutout done.\n");
924
925 // Cutout Continuum
926
927 if (writeCutout(sharpRec, contRec, mharpRec, "continuum")) {
928 printf("Continuum cutout failed\n");
929 return 1;
930 }
931 printf("Intensitygram cutout done.\n");
932
933 // Keywords & Links
934 copy_patch_keys(mharpRec, sharpRec);
935 mbobra 1.1 copy_geo_keys(mharpRec, sharpRec);
936 // rename HARPNUM to TARPNUM
937 val = drms_getkey_double(mharpRec, "HARPNUM", &status);
938 drms_setkey_double(sharpRec, "TARPNUM", val);
939 // copy everything else
940 drms_copykey(sharpRec, mharpRec, "T_REC");
941 drms_copykey(sharpRec, mharpRec, "CDELT1");
942 drms_copykey(sharpRec, mharpRec, "RSUN_OBS");
943 drms_copykey(sharpRec, mharpRec, "DSUN_OBS");
944 drms_copykey(sharpRec, mharpRec, "OBS_VR");
945 drms_copykey(sharpRec, mharpRec, "OBS_VW");
946 drms_copykey(sharpRec, mharpRec, "OBS_VN");
947 drms_copykey(sharpRec, mharpRec, "CRLN_OBS");
948 drms_copykey(sharpRec, mharpRec, "CRLT_OBS");
949 drms_copykey(sharpRec, mharpRec, "CAR_ROT");
950 drms_copykey(sharpRec, mharpRec, "SIZE_SPT");
951 drms_copykey(sharpRec, mharpRec, "AREA_SPT");
952 drms_copykey(sharpRec, mharpRec, "DATE__OBS");
953 drms_copykey(sharpRec, mharpRec, "T_OBS");
954 drms_copykey(sharpRec, mharpRec, "T_MAXPIX");
955 drms_copykey(sharpRec, mharpRec, "QUALITY");
956 mbobra 1.1 drms_copykey(sharpRec, mharpRec, "NPIX_SPT");
957 drms_copykey(sharpRec, mharpRec, "ARS_NCLN");
958 drms_copykey(sharpRec, mharpRec, "ARS_MODL");
959 drms_copykey(sharpRec, mharpRec, "ARS_EDGE");
960 drms_copykey(sharpRec, mharpRec, "ARS_BETA");
961 drms_copykey(sharpRec, mharpRec, "T_MID1");
962 drms_copykey(sharpRec, mharpRec, "T_CMPASS");
963
964 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MTARP");
965 if (mHarpLink) {
966 drms_link_set("MTARP", sharpRec, mharpRec);
967 }
968
969 setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices
970 setKeys(sharpRec, mharpRec, NULL); // Set all other keywords, NULL specifies cutout
971
972 // Stats
973
974 int nCutSegs = 3;
975 for (int iSeg = 0; iSeg < 3; iSeg++) {
976 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
977 mbobra 1.1 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
978 set_statistics(outSeg, outArray, 1);
979 drms_free_array(outArray);
980 }
981
982 return 0;
983
984 }
985
986 /*
987 * Get cutout and write segment
988 */
989
990 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName)
991 {
992
993 int status = 0;
994
995 DRMS_Segment_t *inSeg = NULL, *outSeg = NULL;
996 DRMS_Array_t *cutoutArray = NULL;
997 // DRMS_Type_t arrayType;
998 mbobra 1.1
999 int ll[2], ur[2], nx, ny, nxny; // lower-left and upper right coords
1000
1001 /* Info */
1002
1003 inSeg = drms_segment_lookup(inRec, SegName);
1004 if (!inSeg) return 1;
1005 //printf("SegName=%s\n",SegName); fflush(stdout);
1006 nx = (int) drms_getkey_float(harpRec, "CRSIZE1", &status);
1007 ny = (int) drms_getkey_float(harpRec, "CRSIZE2", &status);
1008 nxny = nx * ny;
1009 ll[0] = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1; if (status) return 1;
1010 ll[1] = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1; if (status) return 1;
1011 ur[0] = ll[0] + nx - 1; if (status) return 1;
1012 ur[1] = ll[1] + ny - 1; if (status) return 1;
1013 if (inSeg->axis[0] == nx && inSeg->axis[1] == ny) { // for bitmaps, infomaps, etc.
1014 cutoutArray = drms_segment_read(inSeg, DRMS_TYPE_DOUBLE, &status);
1015 if (status) return 1;
1016 } else if (inSeg->axis[0] == FOURK && inSeg->axis[1] == FOURK) { // for full disk ones
1017 cutoutArray = drms_segment_readslice(inSeg, DRMS_TYPE_DOUBLE, ll, ur, &status);
1018 if (status) return 1;
1019 mbobra 1.1 } else {
1020 return 1;
1021 }
1022 /* Write out */
1023 outSeg = drms_segment_lookup(outRec, SegName);
1024 if (!outSeg) return 1;
1025 outSeg->axis[0] = cutoutArray->axis[0];
1026 outSeg->axis[1] = cutoutArray->axis[1];
1027 cutoutArray->israw = 0; // always compressed
1028 cutoutArray->bzero = outSeg->bzero;
1029 cutoutArray->bscale = outSeg->bscale; // Same as inArray's
1030 status = drms_segment_write(outSeg, cutoutArray, 0);
1031 drms_free_array(cutoutArray);
1032 if (status) return 1;
1033 //printf("line1068\n"); fflush(stdout);
1034 return 0;
1035
1036 }
1037
1038
1039 /*
1040 mbobra 1.1 * Compute space weather indices
1041 */
1042
1043 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo)
1044 {
1045
1046 int status = 0;
1047 int nx = mInfo->ncol, ny = mInfo->nrow;
1048 int nxny = nx * ny;
1049 int dims[2] = {nx, ny};
1050
1051 // Get bx, by, bz, mask
1052
1053 // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1054 DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
1055 DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
1056 int *bitmask = (int *) bitmaskArray->data; // get the previously made mask array
1057
1058 //Use magnetogram map to compute R
1059 DRMS_Segment_t *losSeg = drms_segment_lookup(inRec, "magnetogram");
1060 DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
1061 mbobra 1.1 float *los = (float *) losArray->data; // los
1062
1063 // Get emphemeris
1064 float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);
1065 float dsun_obs = drms_getkey_float(inRec, "DSUN_OBS", &status);
1066 double rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status);
1067 double rsun_obs = drms_getkey_double(inRec, "RSUN_OBS", &status);
1068 float imcrpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
1069 float imcrpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
1070 float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status);
1071 float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status);
1072
1073 // convert cdelt1_orig from degrees to arcsec
1074 float cdelt1 = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
1075
1076 // Temp arrays
1077 float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1078 float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1079
1080 // define some values for the R calculation
1081 int scale = round(2.0/cdelt1);
1082 mbobra 1.1 int nx1 = nx/scale;
1083 int ny1 = ny/scale;
1084 int nxp = nx1+40; // same comment as above
1085 int nyp = ny1+40; // why is this a +40 pixel size? is this an MDI pixel?
1086 float *rim = (float *)malloc(nx1*ny1*sizeof(float));
1087 float *p1p0 = (float *)malloc(nx1*ny1*sizeof(float));
1088 float *p1n0 = (float *)malloc(nx1*ny1*sizeof(float));
1089 float *p1p = (float *)malloc(nx1*ny1*sizeof(float));
1090 float *p1n = (float *)malloc(nx1*ny1*sizeof(float));
1091 float *p1 = (float *)malloc(nx1*ny1*sizeof(float));
1092 float *pmap = (float *)malloc(nxp*nyp*sizeof(float));
1093 float *p1pad = (float *)malloc(nxp*nyp*sizeof(float));
1094 float *pmapn = (float *)malloc(nx1*ny1*sizeof(float));
1095
1096
1097 // THREE spaceweather quantities computed: USFLUX, MEANGBZ, R_VALUE
1098 if (computeAbsFlux(los, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1099 &(swKeys_ptr->count_mask), bitmask, cdelt1, rsun_ref, rsun_obs))
1100 {
1101 swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1102 swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1103 mbobra 1.1 swKeys_ptr->count_mask = DRMS_MISSING_INT;
1104 }
1105
1106
1107 if (computeBzderivative(los, dims, &(swKeys_ptr->mean_derivative_bz),
1108 bitmask, derx_bz, dery_bz))
1109 {
1110 swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1111 }
1112
1113
1114 if (computeR(los, dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,
1115 p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
1116 {
1117 swKeys_ptr->Rparam = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1118 }
1119
1120
1121 // Clean up the arrays
1122 drms_free_array(bitmaskArray);
1123 //drms_free_array(bzArray);
1124 mbobra 1.1 drms_free_array(losArray);
1125
1126 // free arrays related to Bz derivative
1127 free(derx_bz); free(dery_bz);
1128 // free the arrays that are related to the r calculation
1129 free(rim);
1130 free(p1p0);
1131 free(p1n0);
1132 free(p1p);
1133 free(p1n);
1134 free(p1);
1135 free(pmap);
1136 free(p1pad);
1137 free(pmapn);
1138 }
1139
1140 /*
1141 * Set space weather indices
1142 */
1143
1144 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
1145 mbobra 1.1 {
1146 drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
1147 drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
1148 drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);
1149 drms_setkey_float(outRec, "CMASK", swKeys_ptr->count_mask);
1150 };
1151
1152 /*
1153 * Set all keywords, no error checking for now
1154 */
1155
1156 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, struct mapInfo *mInfo)
1157 {
1158
1159 int status = 0;
1160
1161 // Change a few geometry keywords for CEA & cutout records
1162 if (mInfo != NULL) { // CEA
1163 printf("Calculating CEA keys\n");
1164 drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
1165 drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
1166 mbobra 1.1 drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
1167 drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
1168 drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
1169 drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
1170 drms_setkey_string(outRec, "CUNIT1", "degree");
1171 drms_setkey_string(outRec, "CUNIT2", "degree");
1172 char key[64];
1173 snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
1174 drms_setkey_string(outRec, "CTYPE1", key);
1175 snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
1176 drms_setkey_string(outRec, "CTYPE2", key);
1177 drms_setkey_float(outRec, "CROTA2", 0.0);
1178 // Set BUNIT for each segment
1179 int nSeg = 3;
1180 for (int iSeg = 0; iSeg < nSeg; iSeg++) {
1181 DRMS_Segment_t *outSeg = NULL;
1182 outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);
1183 if (!outSeg) continue;
1184 char bunit_xxx[20];
1185 sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
1186 //printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);
1187 mbobra 1.1 drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);
1188 }
1189
1190 } else { // Cutout
1191
1192 float disk_xc, disk_yc;
1193 disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);
1194 disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);
1195 float x_ll = drms_getkey_float(mharpRec, "CRPIX1", &status);
1196 float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &status);
1197 // Defined as disk center's pixel address wrt lower-left of cutout
1198 drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
1199 drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
1200 // Always 0.
1201 drms_setkey_float(outRec, "CRVAL1", 0);
1202 drms_setkey_float(outRec, "CRVAL2", 0);
1203
1204 // Jan 2 2014 XS
1205 int nSeg = ARRLENGTH(CutSegs);
1206 for (int iSeg = 0; iSeg < nSeg; iSeg++)
1207 {
1208 mbobra 1.1 DRMS_Segment_t *outSeg = NULL;
1209 outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);
1210 if (!outSeg) continue;
1211 // Set Bunit
1212 char bunit_xxx[20];
1213 sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
1214 //printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);
1215 drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);
1216 }
1217 }
1218
1219
1220 TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1221 tnow = (double)time(NULL);
1222 tnow += UNIX_epoch;
1223
1224 val = drms_getkey_time(mharpRec, "DATE", &status);
1225 drms_setkey_time(outRec, "DATE", tnow);
1226
1227 // set cvs commit version into keyword CODEVER7
1228 char *cvsinfo = strdup("$Id");
1229 mbobra 1.1 char *cvsinfo2 = smarp_functions_version();
1230 char cvsinfoall[2048];
1231 strcat(cvsinfoall,cvsinfo);
1232 strcat(cvsinfoall,"\n");
1233 strcat(cvsinfoall,cvsinfo2);
1234 status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
1235
1236 };
1237
1238 /* ############# Nearest neighbour interpolation ############### */
1239
1240 float nnb (float *f, int nx, int ny, double x, double y)
1241 {
1242
1243 if (x <= -0.5 || y <= -0.5 || x > nx - 0.5 || y > ny - 0.5)
1244 return DRMS_MISSING_FLOAT;
1245 int ilow = floor (x);
1246 int jlow = floor (y);
1247 int i = ((x - ilow) > 0.5) ? ilow + 1 : ilow;
1248 int j = ((y - jlow) > 0.5) ? jlow + 1 : jlow;
1249 return f[j * nx + i];
1250 mbobra 1.1
1251 }
1252
1253 /* ################## Wrapper for Jesper's rebin code ################## */
1254
1255 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)
1256 {
1257
1258 struct fresize_struct fresizes;
1259 int nxout, nyout, xoff, yoff;
1260 int nlead = nx;
1261
1262 nxout = nx / nbin; nyout = ny / nbin;
1263 if (gauss && nbin != 1)
1264 init_fresize_gaussian(&fresizes, (nbin / 2), (nbin / 2 * 2), nbin); // for nbin=3, sigma=1, half truncate width=2
1265 else
1266 init_fresize_bin(&fresizes, nbin);
1267 xoff = nbin / 2 + nbin / 2;
1268 yoff = nbin / 2 + nbin / 2;
1269 fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);
1270
1271 mbobra 1.1 }
|
Return to
smarp.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps