1 tplarson 1.1 /* this module is based on jv2helio, but instead interpolates its input
2 * to the desired resolution, optionally correcting for distortion and
3 * errors in the p-angle and carrington inclination.
4 */
5
6
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <time.h>
10 #include <math.h>
11 #include <sys/time.h>
12 #include <sys/resource.h>
13 #include <mkl_lapack.h> //needed for dsecnd()
14 #include "jsoc_main.h"
15 #include "fstats.h"
16 #include "drms_dsdsapi.h"
|
17 tplarson 1.3 #include "projection.h"
|
18 tplarson 1.1
19 #define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0]))
20 #define QUAL_NODATA (0x80000000)
21
22 #define PI (M_PI)
23 #define RADSINDEG (PI/180)
24 #define ARCSECSINRAD (3600*180/PI)
25 #define DAYSINYEAR (365.2425)
26 #define SECSINDAY (86400)
27 #define TAU_A (499.004783806) // light travel time in seconds, = 1 AU/c
28 //#define TAU_A (499.004782) // this value used in old v2helio
29 #define TCARR (25.38) // days
30 #define RTRUE (6.96000000e8) // meters
31 #define AU (149597870691) // meters/au
32 //#define AU (1.49597870e11) // this value used in old v2helio
33 #define MAXLEN (256)
34 #define NO_DATASET (-1)
35 #define NO_IMAGE (-1)
36 #define kMAX_SKIPERRMSG 1024
37 #define kMAXROWS 65536
38
39 tplarson 1.1 #define kNOTSPECIFIED "not specified"
40
41 /* Must do decls */
42 double dsecnd();
43
44 typedef enum
45 {
46 V2HStatus_Success,
47 V2HStatus_MissingParameter,
48 V2HStatus_IllegalTimeFormat,
49 V2HStatus_TimeConvFailed,
50 V2HStatus_Unimplemented,
51 V2HStatus_IllegalOrientation
52 } V2HStatus_t;
53
54 char *module_name = "undistortmdi";
|
55 tplarson 1.4 char *cvsinfo_undistortmdi = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/undistortmdi.c,v 1.3 2013/05/03 20:43:25 tplarson Exp $";
|
56 tplarson 1.1
57 ModuleArgs_t module_args[] =
58 {
59 {ARG_STRING, "in", "", "input data records"},
60 {ARG_STRING, "out", "", "output data series"},
61 {ARG_STRING, "segin", kNOTSPECIFIED, "name of input segment if not using segment 0"},
62 {ARG_STRING, "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
63 {ARG_STRING, "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
64 {ARG_STRING, "srclink", "SRCDATA", "name of link to source data"},
65 {ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
66 {ARG_INT, "VERB", "1", "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", ""}, /* debug messages */
67 {ARG_INT, "OUTCOLS", "4096", "number of output columns"},
68 {ARG_INT, "OUTROWS", "4096", "number of output rows"},
|
69 tplarson 1.2 {ARG_FLOAT, "MAPRMAX", "1.1", "maximum image radius", ""},
|
70 tplarson 1.1 {ARG_INT, "RMAXFLAG", "1", "set to nonzero to set output to DRMS_MISSING outside MAPRMAX, otherwise uses 0.0 outside MAPRMAX"}, /* Non 0 sets data outside RMAX MISSING */
71 {ARG_INT, "INTERPO", "1", "option for interpolation: 0=bilinear; 1=cubic convolution", ""}, /* 2 methods - see soi_fun.h */
72 {ARG_FLOAT, "OUTSCALE", "1.0", "bscale to use for output", ""}, /* scale for output */
73 {ARG_FLOAT, "OUTBIAS", "0.0", "bzero to use for output", ""}, /* bias for scaled output */
74 {ARG_INT, "DISTORT", "0", "option for distortion correction: 0=none; 1=full disk(fd) data; 2=vector-weighted(vw) data", ""}, /* 0 for none, 1 for FD, 2 for vw */
75 {ARG_FLOAT, "CUBIC", "7.06E-9", "cubic distortion in fd units", ""}, /* Cubic distortion in FD units */
76 {ARG_FLOAT, "TILTALPHA","2.59", "tilt of CCD, degrees", ""}, /* TILT of CCD in degrees */
77 {ARG_FLOAT, "TILTBETA", "56.0", "direction of CCD tilt, degrees", ""}, /* Direction of TILT in degrees */
78 {ARG_FLOAT, "TILTFEFF", "12972.629", "effective focal length", ""}, /* Effective focal length */
79 {ARG_INT, "DATASIGN", "0", "value to multiply data; set to 0 to take DATASIGN from keyword, or 1.0 if not found", ""}, /* Non 0 forces datasign to value*/
80 {ARG_INT, "MAXMISSVALS", "0", "if >0, this becomes threshold on MISSVALS from keyword", ""}, /* max. allowed MISSING pixels */
81 {ARG_FLOAT, "SOLAR_P", "999.0", "P-angle; if unset, taken from keywords", ""}, /* can't use D_MISSING here */
82 {ARG_FLOAT, "PSIGN", "1.0", "sign of SOLAR_P", ""}, /* Sign of P. For MWO data. */
83 {ARG_FLOAT, "PERR", "0.0", "fixed P-angle error, likely -0.22", ""}, /* Fixed P-angle error. Maybe -0.22. */
84 {ARG_FLOAT, "IERR", "0.0", "error in Carrington inclination, likely -0.10", ""}, /* Error in Carrington inclination. Maybe -0.10. */
85 {ARG_TIME, "REF_TB0", "2001.06.06_06:57:22_TAI", "reference time for computing correction to P and B angles, roughly when B0=0", ""},
86 {ARG_END}
87 };
88
|
89 tplarson 1.3 #include "imageinterp.c"
|
90 tplarson 1.1 #include "saveparm.c"
|
91 tplarson 1.3 #include "timing.c"
92 #include "set_history.c"
|
93 tplarson 1.1
94 static inline void ParameterDef(int status,
95 char *pname,
96 double defaultp,
97 double *p,
98 int iRec,
99 int verbflag)
100 {
101 /* logs warning and sets parameter to default value */
102 if (status != 0)
103 {
104 *p = defaultp;
105 if (verbflag)
106 fprintf(stderr, "WARNING: default value %g used for %s, iRec = %d, status = %d\n", defaultp, pname, iRec, status);
107 }
108 }
109
110
111 #define PARAMETER_ERROR(PNAME) \
112 if (status != DRMS_SUCCESS) \
113 { \
114 tplarson 1.1 CreateBlankRecord(inrec, outrec, quality); \
115 fprintf(stderr, \
116 "SKIP: error getting keyword %s: iRec = %d, status = %d, T_REC = %s, recnum = %lld\n", \
117 PNAME, \
118 iRec, \
119 status, \
120 trecstr, \
121 inrec->recnum); \
122 if (inarr) drms_free_array(inarr); \
123 if (outarr) drms_free_array(outarr); \
124 if (orientation) free(orientation); \
125 continue; \
126 }
127
128
129 /* Segment will be empty, but there will be a record! */
130 static void CreateBlankRecord(DRMS_Record_t *inrec, DRMS_Record_t *outrec, int quality)
131 {
132 /* create 'blank' data */
133 // might insert 'quality = quality | MASK' here.
134 quality = quality | QUAL_NODATA;
135 tplarson 1.1 drms_copykey(outrec, inrec, "T_REC");
136 drms_setkey_int(outrec, "QUALITY", quality);
137 drms_close_record(outrec, DRMS_INSERT_RECORD);
138 }
139
140
141 int DoIt(void)
142 {
143
144 int newstat = 0;
145 int status = DRMS_SUCCESS;
146 V2HStatus_t vstat = V2HStatus_Success;
147
148 const char *orientationdef = "SESW ";
149 char *orientation = NULL;
150 int paramsign;
151 int longitude_shift, velocity_correction, interpolation, apodization;
152 int mag_correction;
153 int mag_offset;
154 int row;
155 int mapped_lmax, sinb_divisions, mapcols, maprows, nmax, nmin;
156 tplarson 1.1 int length[2];
157 int carrStretch = 0;
158 float diffrotA = 0.0;
159 float diffrotB = 0.0;
160 float diffrotC = 0.0;
161 double tobs, tmearth, tref, trefb0;
162 double smajor, sminor, sangle;
163 double xscale, yscale, imagescale;
164 int xpixels, ypixels, pixels;
165 double obs_vr, obs_vw, obs_vn;
166 double b0, bmax, bmin, sinBdelta;
167 double longmax, longmin, longmax_adjusted, longmin_adjusted, longinterval;
168 double p0, p, rmax;
169 double ierr, perr, psign;
170 double x0, y0;
171 double obsdist, longshift, obsl0, refl0, mapl0, longrate, rtrue, rsun, S;
172 double rsunDef;
173 int obsCR;
174 int apel;
175 double apinner, apwidth, apx, apy;
176 double scale, bias;
177 tplarson 1.1 double colsperdeg;
178
179 int quality, NaN_beyond_rmax, nRecs;
180 double satrot, instrot;
181 double dsignout, vsign;
182 int distsave;
183 double cubsave, tiltasave, tiltbsave, tiltfsave;
184 TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
185 char trecstr[100];
186
|
187 tplarson 1.3 LIBPROJECTION_Dist_t distP;
|
188 tplarson 1.1 DRMS_Segment_t *segin = NULL;
189 DRMS_Segment_t *segout = NULL;
190 DRMS_Array_t *inarr = NULL;
191 DRMS_Array_t *outarr = NULL;
192 DRMS_Record_t *inrec = NULL;
193 DRMS_Record_t *outrec = NULL;
194 DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
195 DRMS_RecLifetime_t lifetime;
196
197 long long histrecnum=-1;
198
199 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
200 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
201
202 struct timeval tv0, tv1, tv;
203 double ct0, ct1, ct;
204 double wt0, wt1, wt;
205 double ut0, ut1, ut;
206 double st0, st1, st;
207 double tt0, tt1, tt;
208
209 tplarson 1.1 wt0=getwalltime();
210 getcputime(&ut0, &st0);
211
212 gettimeofday(&tv0, NULL);
213 ct0=dsecnd();
214 tt0=1000*((double)clock())/CLOCKS_PER_SEC;
215
216 char *inRecQuery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
217 char *outSeries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
218
219 int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
220 int maxmissvals = cmdparams_save_int(&cmdparams, "MAXMISSVALS", &newstat);
221 int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
222 if (permflag)
223 lifetime = DRMS_PERMANENT;
224 else
225 lifetime = DRMS_TRANSIENT;
226
227
228 interpolation = cmdparams_save_int(&cmdparams, "INTERPO", &newstat);
229 paramsign = cmdparams_save_int(&cmdparams, "DATASIGN", &newstat);
230 tplarson 1.1
231 distsave = cmdparams_save_int(&cmdparams, "DISTORT", &newstat);
232 cubsave = cmdparams_save_double(&cmdparams, "CUBIC", &newstat);
233 tiltasave = cmdparams_save_double(&cmdparams, "TILTALPHA", &newstat);
234 tiltbsave = cmdparams_save_double(&cmdparams, "TILTBETA", &newstat);
235 tiltfsave = cmdparams_save_double(&cmdparams, "TILTFEFF", &newstat);
236
237 scale = cmdparams_save_double(&cmdparams, "OUTSCALE", &newstat);
238 bias = cmdparams_save_double(&cmdparams, "OUTBIAS", &newstat);
239 p0 = cmdparams_save_double(&cmdparams, "SOLAR_P", &newstat);
240 psign = cmdparams_save_double(&cmdparams, "PSIGN", &newstat);
241 perr = cmdparams_save_double(&cmdparams, "PERR", &newstat);
242 ierr = cmdparams_save_double(&cmdparams, "IERR", &newstat);
243 trefb0 = cmdparams_save_time(&cmdparams, "REF_TB0", &newstat);
244 rmax = cmdparams_save_double(&cmdparams, "MAPRMAX", &newstat);
245 NaN_beyond_rmax = cmdparams_save_int(&cmdparams, "RMAXFLAG", &newstat);
246
247 SetDistort(distsave, cubsave, tiltasave, tiltbsave, tiltfsave, &distP);
248
249 mapcols = cmdparams_save_int(&cmdparams, "OUTCOLS", &newstat);
250 maprows = cmdparams_save_int(&cmdparams, "OUTROWS", &newstat);
251 tplarson 1.1
252 length[0] = mapcols;
253 length[1] = maprows;
254
255 char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
256 char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
257 int seginflag = strcmp(kNOTSPECIFIED, segnamein);
258 int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
259
260 char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
261 char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
262
|
263 tplarson 1.3 if (newstat)
264 {
265 fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
266 cpsave_decode_error(newstat);
267 return 1;
268 }
269 else if (savestrlen != strlen(savestr))
270 {
271 fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
272 return 1;
273 }
274
|
275 tplarson 1.1 // set up ancillary dataseries for processing metadata
276 DRMS_Record_t *tempoutrec = drms_create_record(drms_env,
277 outSeries,
278 DRMS_TRANSIENT,
279 &status);
280
281 if (status != DRMS_SUCCESS)
282 {
283 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outSeries, status);
284 return 1;
285 }
286
287 DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
288 DRMS_Link_t *srclink = hcon_lookup_lower(&tempoutrec->links, srclinkname);
289
290 if (histlink != NULL)
291 {
|
292 tplarson 1.3 histrecnum=set_history(histlink);
293 if (histrecnum < 0)
294 {
295 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
296 return 1;
297 }
|
298 tplarson 1.1 }
299 else
300 {
|
301 tplarson 1.3 fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
|
302 tplarson 1.1 }
303
304 int itest;
305 // these must be present in the output dataseries and variable, not links or constants
306 /*
307 char *outchecklist[] = {"T_REC", "QUALITY", "DATASIGN",
308 "CRPIX1", "CRVAL1", "CROTA1", "CDELT1",
309 "CRPIX2", "CRVAL2", "CROTA2", "CDELT2" };
310
311 for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
312 {
313 DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
314 if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1)
315 {
316 fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
317 return 1;
318 }
319 }
320 */
321 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
322
323 tplarson 1.1
324 DRMS_RecordSet_t *inRecSet = drms_open_records(drms_env, inRecQuery, &status);
325 if (status != DRMS_SUCCESS || inRecSet == NULL)
326 {
327 fprintf(stderr, "ERROR: problem opening input recordset: status = %d\n", status);
328 return 1;
329 }
330
331 nRecs = inRecSet->n;
332 if (nRecs == 0)
333 {
334 printf("WARNING: input recordset contains no records\nmodule %s successful completion\n", cmdparams.argv[0]);
335 return 1;
336 }
337
338 if (verbflag)
339 printf("input recordset opened, nRecs = %d\n",nRecs);
340
341
342 // go ahead and check for the presence of these input keywords as well
343 char *inchecklist[] = {"T_REC", "QUALITY", "T_OBS", "CRLT_OBS", "CRLN_OBS",
344 tplarson 1.1 "CDELT1", "CDELT2"};
345
346 DRMS_Record_t *tempinrec = inRecSet->records[0];
347 for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
348 {
349 DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&tempinrec->keywords, inchecklist[itest]);
350 if (!inkeytest)
351 {
352 fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
353 return 1;
354 }
355 }
356
357
358 int iRec;
359 int error=0; // only set error before a break
360 int nsuccess=0;
361 for (iRec = 0; iRec < nRecs; iRec++)
362 {
363 if (verbflag > 1)
364 {
|
365 tplarson 1.3 wt1=getwalltime();
366 getcputime(&ut1, &st1);
|
367 tplarson 1.1 gettimeofday(&tv1, NULL);
368 ct1=dsecnd(); //((double)clock())/CLOCKS_PER_SEC;
|
369 tplarson 1.3 tt1=1000*((double)clock())/CLOCKS_PER_SEC;
|
370 tplarson 1.1 printf("processing record %d\n", iRec);
371 }
372 inrec = inRecSet->records[iRec];
373
374 // create an output record
375 outrec = drms_create_record(drms_env,
376 outSeries,
377 lifetime,
378 &status);
379
380 if (status != DRMS_SUCCESS || outrec==NULL)
381 {
382 fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status);
383 error=2;
384 break;
385 }
386
387 if (histlink)
388 drms_setlink_static(outrec, histlinkname, histrecnum);
389
390 if (srclink)
391 tplarson 1.1 drms_setlink_static(outrec, srclinkname, inrec->recnum);
392
393
394 if (segoutflag)
395 segout = drms_segment_lookup(outrec, segnameout);
396 else
397 segout = drms_segment_lookupnum(outrec, 0);
398
399 // create an output array
400 outarr = drms_array_create(usetype, 2, length, NULL, &status);
401
402 if (!segout || !outarr || status != DRMS_SUCCESS)
403 {
404 fprintf(stderr, "ERROR: problem with output segment or array: iRec = %d, status = %d\n", iRec, status);
405 if (outarr)
406 drms_free_array(outarr);
407 error=1;
408 break;
409 }
410
411 trec = drms_getkey_time(inrec, "T_REC", &status);
412 tplarson 1.1 if (status != DRMS_SUCCESS)
413 {
414 fprintf(stderr,
415 "SKIP: error getting prime keyword %s: iRec = %d, status = %d, recnum = %lld\n",
416 "T_REC",
417 iRec,
418 status,
419 inrec->recnum);
420 drms_free_array(outarr);
421 continue;
422 }
423 sprint_time(trecstr, trec, "TAI", 0);
424
425 quality = drms_getkey_int(inrec, "QUALITY", &status);
426 PARAMETER_ERROR("QUALITY")
427 // insert tests on quality here.
428 // if we encounter a keyword error causing a continue within the loop, we should set a bit in quality for this.
429 // in other words, CreateBlankRecord should always be preceded by setting quality, or could move this inside CreateBlankRecord.
430 if (quality & QUAL_NODATA)
431 {
432 CreateBlankRecord(inrec, outrec, quality);
433 tplarson 1.1 fprintf(stderr,"SKIP: record rejected based on quality = %d: iRec = %d, T_REC = %s, recnum = %lld\n", quality, iRec, trecstr, inrec->recnum);
434 drms_free_array(outarr);
435 continue;
436 }
437
438 if (seginflag)
439 segin = drms_segment_lookup(inrec, segnamein);
440 else
441 segin = drms_segment_lookupnum(inrec, 0);
442
443 if (segin)
444 inarr = drms_segment_read(segin, usetype, &status);
445
446 if (!segin || !inarr || status != DRMS_SUCCESS)
447 {
448 // CreateBlankRecord(inrec, outrec, quality);
449 fprintf(stderr, "ERROR: problem with input segment or array: iRec = %d, status = %d, T_REC = %s, recnum = %lld \n", iRec, status, trecstr, inrec->recnum);
450 drms_free_array(outarr);
451 if (inarr)
452 drms_free_array(inarr);
453 error=1;
454 tplarson 1.1 break;
455 }
456
457 if (maxmissvals > 0)
458 {
459 int missvals = drms_getkey_int(inrec, "MISSVALS", &status);
460 if (status == DRMS_ERROR_UNKNOWNKEYWORD)
461 {
462 fprintf(stderr, "ERROR: required keyword %s missing, iRec = %d, T_REC = %s, recnum = %lld\n", "MISSVALS", iRec, trecstr, inrec->recnum);
463 drms_free_array(inarr);
464 drms_free_array(outarr);
465 free(orientation);
466 error=1;
467 break;
468 }
469 PARAMETER_ERROR("MISSVALS")
470 if (missvals > maxmissvals)
471 {
472 CreateBlankRecord(inrec, outrec, quality);
473 fprintf(stderr, "SKIP: %d pixels MISSING, max allowed is %d: iRec = %d, status = %d, T_REC = %s, recnum = %lld\n", missvals, maxmissvals, iRec, status, trecstr, inrec->recnum);
474 drms_free_array(inarr);
475 tplarson 1.1 drms_free_array(outarr);
476 free(orientation);
477 continue;
478 }
479 }
480
481 tobs = drms_getkey_time(inrec, "T_OBS", &status);
482 PARAMETER_ERROR("T_OBS")
483
484 // MDI keyword was OBS_B0
485 b0 = drms_getkey_double(inrec, "CRLT_OBS", &status);
486 PARAMETER_ERROR("CRLT_OBS")
487
488 // MDI keyword was OBS_L0
489 obsl0 = drms_getkey_double(inrec, "CRLN_OBS", &status);
490 PARAMETER_ERROR("CRLN_OBS")
491
492 if (p0 == 999.0)
493 {
494 // MDI keyword was SOLAR_P = -(SAT_ROT + INST_ROT)
495 /*
496 tplarson 1.1 satrot = drms_getkey_double(inrec, "SAT_ROT", &status);
497 PARAMETER_ERROR("SAT_ROT")
498 instrot = drms_getkey_double(inrec, "INST_ROT", &status);
499 PARAMETER_ERROR("INST_ROT")
500 p=-(satrot+instrot);
501 */
502 double crota = drms_getkey_double(inrec, "CROTA2", &status);
503 PARAMETER_ERROR("CROTA2")
504 p=-crota;
505 }
506 else
507 {
508 p = p0;
509 }
510
511 // fix for 1988 MWO
512 // p = 180.0 - p ;
513
514 p = psign * p ;
515
516 // 991839442. corresponds to hour 73878 minute 57 second 22
517 tplarson 1.1 // or 73878.956 or day 3078.2898, roughly when B0 is 0
518 // b0=b0 * 0.986207; One way of correcting.
519 // The following is pretty good
520 // b0=b0-0.1*sin((tobs-991839442.)/31557600.*2*PI);
521 // b0 = b0 + ierr * sin((tobs - 991839442.) / 31557600. * 2 * PI);
522 b0 = b0 + ierr * sin((tobs - trefb0) / 31557600. * 2 * PI);
523
524 // p=p-0.2;
525 // p=p-0.2+0.1*cos((tobs-991839442.)/31557600.*2*PI);
526 // p = p + perr - ierr * cos((tobs - 991839442.) / 31557600. * 2 * PI);
527 p = p + perr - ierr * cos((tobs - trefb0) / 31557600. * 2 * PI);
528
529 if (paramsign != 0)
530 {
531 vsign = paramsign;
532 }
533 else
534 {
535 vsign = drms_getkey_double(inrec, "DATASIGN", &status);
536 ParameterDef(status, "DATASIGN", 1.0, &vsign, iRec, 1);
537 }
538 tplarson 1.1
539 xscale = drms_getkey_double(inrec, "CDELT1", &status);
540 PARAMETER_ERROR("CDELT1")
541 yscale = drms_getkey_double(inrec, "CDELT2", &status);
542 PARAMETER_ERROR("CDELT2")
543
544 if (xscale != yscale)
545 {
546 fprintf(stderr, "ERROR: %s != %s not supported, iRec = %d, T_REC = %s, recnum = %lld \n", "CDELT1", "CDELT2", iRec, trecstr, inrec->recnum);
547 drms_free_array(inarr);
548 error++;
549 continue;
550 }
551 imagescale=xscale;
552
553 // MDI keyword was OBS_DIST, in AU
554 obsdist = drms_getkey_double(inrec, "DSUN_OBS", &status);
555 obsdist /= AU; // note that an incorrect value of 1.49597892e11 has sometimes been used to convert between OBS_DIST and DSUN_OBS when porting data from DSDS to DRMS
556 ParameterDef(status, "DSUN_OBS", 1.0, &obsdist, iRec, 1);
557 S = rtrue / obsdist; // radians - approx. arcsin(rtrue/obsdist)
558
559 tplarson 1.1
560 rsun = drms_getkey_double(inrec, "R_SUN", &status);
561 if (status != DRMS_SUCCESS)
562 {
563 rsun = drms_getkey_double(inrec, "RSUN_OBS", &status);
564 if (status != DRMS_SUCCESS)
565 rsun = ARCSECSINRAD * S / sqrt(1.0 - S * S) / imagescale;
566 else
567 rsun /= imagescale;
568 }
569
570 xpixels = inarr->axis[0];
571 ypixels = inarr->axis[1];
572
573 // MDI keyword was X0
574 x0 = drms_getkey_double(inrec, "CRPIX1", &status);
575 ParameterDef(status, "CRPIX1", xpixels / 2, &x0, iRec, 1);
576 x0 -= 1.0;
577
578 // MDI keyword was Y0
579 y0 = drms_getkey_double(inrec, "CRPIX2", &status);
580 tplarson 1.1 ParameterDef(status, "CRPIX2", ypixels / 2, &y0, iRec, 1);
581 y0 -= 1.0;
582
583 if (status = imageinterp((float *)inarr->data,
584 (float *)outarr->data,
585 xpixels,
586 ypixels,
587 x0,
588 y0,
589 p * RADSINDEG,
590 rsun,
591 rmax,
592 NaN_beyond_rmax,
593 interpolation,
594 mapcols,
595 maprows,
596 vsign,
597 &distP))
598 {
599 CreateBlankRecord(inrec, outrec, quality);
600 fprintf(stderr, "failure in imageinterp: iRec = %d, status = %d, T_REC = %s, recnum = %lld\n", iRec, status, trecstr, inrec->recnum);
601 tplarson 1.1 drms_free_array(inarr);
602 drms_free_array(outarr);
603 free(orientation);
604 continue; // go to next image
605 }
606
|
607 tplarson 1.4 drms_setkey_int(outrec, "TOTVALS", maprows*mapcols);
|
608 tplarson 1.1 set_statistics(segout, outarr, 1);
609
610 // Write to the output record
611 outarr->bzero=bias;
612 outarr->bscale=scale;
613 segout->axis[0]=outarr->axis[0]; // required for vardim output
614 segout->axis[1]=outarr->axis[1];
615 drms_segment_write(segout, outarr, 0);
616
617 // drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
618 drms_copykey(outrec, inrec, "T_REC");
619 drms_setkey_int(outrec, "QUALITY", quality);
620
621 double xratio = (double)mapcols/xpixels;
622 double yratio = (double)maprows/ypixels;
623
624 drms_setkey_float(outrec, "CRPIX1", (drms_getkey_float(inrec, "CRPIX1", &status)-0.5)*xratio+0.5);
625 drms_setkey_float(outrec, "CRPIX2", (drms_getkey_float(inrec, "CRPIX2", &status)-0.5)*yratio+0.5);
626 drms_setkey_float(outrec, "CDELT1", drms_getkey_float(inrec, "CDELT1", &status)/xratio);
627 drms_setkey_float(outrec, "CDELT2", drms_getkey_float(inrec, "CDELT2", &status)/yratio);
628 drms_setkey_float(outrec, "CROTA2", 0.0);
629 tplarson 1.1 drms_setkey_float(outrec, "X0", (x0 + 0.5)*xratio - 0.5);
630 drms_setkey_float(outrec, "Y0", (y0 + 0.5)*yratio - 0.5);
631 drms_setkey_float(outrec, "IM_SCALE", imagescale/xratio);
632 drms_setkey_float(outrec, "R_SUN", rsun*xratio);
633 drms_setkey_float(outrec, "RSUN_OBS", drms_getkey_float(inrec, "RSUN_OBS", &status)*xratio);
634
635 drms_copykey(outrec, inrec, "CRVAL1");
636 drms_copykey(outrec, inrec, "CRVAL2");
637 drms_copykey(outrec, inrec, "T_OBS");
638 drms_copykey(outrec, inrec, "CADENCE");
639 drms_copykey(outrec, inrec, "EXPTIME");
640 drms_copykey(outrec, inrec, "CRLT_OBS");
641 drms_copykey(outrec, inrec, "CRLN_OBS");
642 drms_copykey(outrec, inrec, "CAR_ROT");
643 drms_copykey(outrec, inrec, "DSUN_OBS");
644 drms_copykey(outrec, inrec, "R_SUN_REF");
645 drms_copykey(outrec, inrec, "OBS_VR");
646 drms_copykey(outrec, inrec, "OBS_VW");
647 drms_copykey(outrec, inrec, "OBS_VN");
648
649 dsignout=vsign*drms_getkey_double(inrec, "DATASIGN", &status);
650 tplarson 1.1 if (status != DRMS_SUCCESS) dsignout=vsign;
651 dsignout/=fabs(dsignout);
652 drms_setkey_int(outrec, "DATASIGN", (int)dsignout);
653
654 tnow = (double)time(NULL);
655 tnow += UNIX_epoch;
656 drms_setkey_time(outrec, "DATE", tnow);
657
658 drms_close_record(outrec, DRMS_INSERT_RECORD);
659
660 if (verbflag > 1)
661 {
|
662 tplarson 1.3 wt=getwalltime();
663 getcputime(&ut, &st);
664 gettimeofday(&tv, NULL);
665 ct=dsecnd(); //((double)clock())/CLOCKS_PER_SEC;
666 tt=1000*((double)clock())/CLOCKS_PER_SEC;
667 fprintf(stdout,
|
668 tplarson 1.1 "record %d done, %f ms wall time, %f ms cpu time\n",
669 iRec,
670 (float)((tv.tv_sec * 1000000.0 + tv.tv_usec -
671 (tv1.tv_sec * 1000000.0 + tv1.tv_usec)) / 1000.0),
672 (ct-ct1)*1000);
673 printf("test: %f ms wall time, %f ms cpu time\n", wt-wt1, (ut+st)-(ut1+st1));
674 printf("clock: %f ms\n",tt-tt1);
675 }
676
677 drms_free_array(inarr);
678 drms_free_array(outarr);
679 free(orientation);
680 nsuccess++;
681
682 } // end loop
683
684 if (inRecSet)
685 drms_close_records(inRecSet, DRMS_FREE_RECORD);
686
687 if (error == 1)
688 drms_close_record(outrec, DRMS_FREE_RECORD);
689 tplarson 1.1
690 wt=getwalltime();
691 getcputime(&ut, &st);
692
693 gettimeofday(&tv, NULL);
694 ct=dsecnd(); //((double)clock())/CLOCKS_PER_SEC;
695 tt=1000*((double)clock())/CLOCKS_PER_SEC;
696 if (verbflag)
697 {
698 printf("number of records processed = %d\n", nsuccess);
699 fprintf(stdout, "total time spent: %f ms wall time, %f ms cpu time\n",
700 (float)((tv.tv_sec * 1000000.0 + tv.tv_usec - (tv0.tv_sec * 1000000.0 + tv0.tv_usec)) / 1000.0),
701 (ct-ct0)*1000);
702 printf("test: %f ms wall time, %f ms cpu time\n", wt-wt0, (ut+st)-(ut0+st0));
703 printf("clock: %f ms\n",tt-tt0);
704 if (!error)
705 printf("module %s successful completion\n", cmdparams.argv[0]);
706 }
707
708 return 0;
709
710 tplarson 1.1 } // end DoIt
|
Return to
undistortmdi.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps