1 tplarson 1.4 /* this JSOC module is a combination of 3 SOI modules: v2helio, helio2mlat, qdotprod
|
2 tplarson 1.1 * ported by Tim Larson
3 * removes check on mapbmin, mapbmax, and maprows, these keywords are no longer carried
4 */
5
6 /*
|
7 tplarson 1.9 * v2helio.c ~soi/(version)/src/modules/v2helio.c
8 *
9 * This module interpolates CCD velocity data to estimate values that
10 * would be obtained at specified equal increments of heliographic
11 * longitude and sine of latitude. Apodization and corrections for
12 * solar rotation and limbshift are included.
13 *
14 * Responsible: Kay Leibrand KLeibrand@solar.Stanford.EDU
15 *
16 * Bugs:
17 * This module is under development. Look for ??? in comments.
18 *
19 * Planned updates:
20 * Decide when and where and if default values should be used.
21 *
22 */
23
24 /*
25 * helio2mlat.c ~soi/(version)/src/modules/helio2mlat.c
26 *
27 * Description:
28 tplarson 1.9 * Adapted by Kay Leibrand from pipeLNU shtfft
29 * fft by map_rows. Use FORTRAN library Cmlib for float version.
30 * extends data to 360 degrees prior to transform but only saves
31 * cols up to lmax. Performs transpose needed prior to dotprod.
32 *
33 * Responsible: Kay Leibrand KLeibrand@solar.Stanford.EDU
34 *
35 * Bugs:
36 * This module is under development. Look for ??? in comments.
37 *
38 * Planned updates:
39 * Restructure to use functions?
40 * Refine parameter definitions and names to be consistent with
41 * new keywords and ancillary data flow.
42 * Fix phase
43 *
44 */
45
46 /*
|
47 tplarson 1.1 * qdotprod.c ~soi/(version)/src/modules/qdotprod.c
48 *
49 * Description:
50 * Conversion of Jesper Schou's FORTRAN q(uick)dotprod to C
51 * from file ~schou/testdot/testdot4c.f
52 * uses FORTRAN functions from blas library
53 * contains optimizations, calculates masks, allows "chunking" in l
54 *
55 * Responsible: Kay Leibrand KLeibrand@solar.Stanford.EDU
56 *
57 * Bugs:
58 * Inadequate checks for consistency between data specifications,
59 * i.e. dataset names, and LMIN, LMAX, LCHUNK parameters.
60 * This module is under development. Look for ??? and TBD in comments.
61 * Normalization in plm's is different from PHS pipeLNU
62 *
63 * Planned updates:
64 * Fix known bugs.
65 * Use a consistent pointer style, i.e. pointer arith vs []
66 * Restructure to use functions.
67 * Refine parameter definitions and names to be consistent with
68 tplarson 1.1 * new keywords and ancillary data flow.
69 *
70 * Revision history is at end of file.
71 */
72
73 /* Normalization of the resulting time-series when norm != 0:
74
75 Let the observed surface behaviour of an oscillation be given by
76 Re(exp(-i\omega t) Y_l^m (\theta,\phi)) with (Y_l^m)^2 normalized
77 to have an average of 1 over the unit sphere (this is 4\pi times
78 the usual definition where the integral is 1). Assume that the
79 whole (4\pi) Sun is observed with no velocity projection factor.
80 Then the resulting time series is given by exp(-i\omega t).
81
82 This is equivalent to preserving the rms value of the real parts
83 of the oscillations.
84
85 And maybe I got the implementation straight.
86
87 */
88
89 tplarson 1.1 #include <stdio.h>
90 #include <math.h>
91 #include <stdlib.h>
92 #include <time.h>
93 #include <sys/time.h>
94 #include <sys/resource.h>
95 #include <fftw3.h>
96 #include "jsoc_main.h"
|
97 tplarson 1.9 #include "fstats.h"
|
98 tplarson 1.1 #include "drms_dsdsapi.h"
99 #include "errlog.h"
|
100 tplarson 1.21 #include "projection.h"
|
101 tplarson 1.23 #include "atoinc.h"
|
102 tplarson 1.1
103 #define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0]))
104 #define PI (M_PI)
105
106 #define absval(x) (((x) < 0) ? (-(x)) : (x))
107 #define minval(x,y) (((x) < (y)) ? (x) : (y))
108 #define maxval(x,y) (((x) < (y)) ? (y) : (x))
109 #define very_small (1.0e-30)
110 #define is_very_small(x) (absval(x) < very_small)
111 #define is_even(x) (!((x)%2))
112 #define is_odd(x) ((x)%2)
113
114 #define RADSINDEG (PI/180)
115 #define ARCSECSINRAD (3600*180/PI)
116 #define DAYSINYEAR (365.2425)
117 #define SECSINDAY (86400)
118 #define TAU_A (499.004783806) // light travel time in seconds, = 1 AU/c
119 //#define TAU_A (499.004782) // this value used in old v2helio
120 #define TCARR (25.38) // days
121 #define RTRUE (6.96000000e8) // meters
122 #define AU (149597870691) // meters/au
123 tplarson 1.1 //#define AU (1.49597870e11) // this value used in old v2helio
124
|
125 tplarson 1.20 #define QUAL_NODATA (0x80000000)
126 #define QUAL_MIXEDCALVER (0x00000001)
|
127 tplarson 1.1 #define kNOTSPECIFIED "not specified"
128
129 char *module_name = "jv2ts";
|
130 tplarson 1.24 char *cvsinfo_jv2ts = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jv2ts.c,v 1.23 2013/07/02 20:32:03 tplarson Exp $";
|
131 tplarson 1.1
132 ModuleArgs_t module_args[] =
133 {
134 // these inputs are common to all modules
135 {ARG_STRING, "in", NULL, "input data records"},
|
136 tplarson 1.9 {ARG_STRING, "tsout", kNOTSPECIFIED, "output data series"},
|
137 tplarson 1.1 {ARG_STRING, "segin", kNOTSPECIFIED, "input data segment"},
138 {ARG_STRING, "segout", kNOTSPECIFIED, "output data segment"},
|
139 tplarson 1.21 {ARG_STRING, "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata. specify \"none\" to suppress warning messages."},
|
140 tplarson 1.1 {ARG_STRING, "srclink", "SRCDATA", "name of link to source data"}, // used only for jv2helio and jhelio2mlat output
141 {ARG_STRING, "v2hout", kNOTSPECIFIED, "output data series for jv2helio"},
142 {ARG_STRING, "h2mout", kNOTSPECIFIED, "output data series for jhelio2mlat"},
143 {ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
|
144 tplarson 1.8 {ARG_INT, "VERB", "1", "option for level of verbosity: 0 gives only error and warning messages; >0 prints messages outside of loop; >1 prints messages inside loop; >2 for debugging output"},
145 {ARG_INT, "FORCEOUTPUT", "0", "option to specify behavior on missing inputs; 0 gives an error on missing or duplicate inputs, >0 makes outputs regardless"},
146 {ARG_STRING, "TAG", "none", "this parameter sets a keyword of the same name to the same value, usually used as an additional primekey"},
|
147 tplarson 1.23 {ARG_STRING, "VERSION", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, useful for selecting obsolete records"},
|
148 tplarson 1.20 {ARG_INT, "CALVERKEY", "2", "short integer bit mask determining which 4-bit fields of CALVER64 are permitted to change on input. set the bit to disallow change of the corresponding nybble."},
|
149 tplarson 1.1 // these are from jqdotprod
150 {ARG_INT, "LMIN", "0", "minimum l in output"},
151 // {ARG_INT, "LMAX", "0", "maximum l in output, take from input if 0"}, /* if 0, default is LMAX of in_sds */
152 {ARG_INT, "LCHUNK", "0", "increment in l on output, default is lmax+1"}, /* if 0, is LMAX+1 */
153 {ARG_INT, "NORM", "1", "set to nonzero to use cnorm=sinbdelta*sqrt(2) in sgemm call, otherwise use cnorm=1"}, /* Uses old norm if =0, see below */
|
154 tplarson 1.8 {ARG_TIME, "TSTART", NULL, "start time of first output record"},
155 {ARG_STRING, "TTOTAL", NULL, "total length of time in output"},
|
156 tplarson 1.1 {ARG_STRING, "TCHUNK", kNOTSPECIFIED, "length of output timeseries"},
157 // these are from jhelio2mlat
158 {ARG_INT, "LMAX", "300", "maximum l (maximum m) in the output, cannot be greater than MAPMMAX", NULL},
159 {ARG_INT, "SUBMEAN", "0", "nonzero subtracts mean of input image", NULL},
160 {ARG_INT, "NORMLIZE", "0", "nonzero multiplies by sqrt((fraction nonmissing)/2) for each row", NULL},
161 {ARG_INT, "CENTLONG", "1", "nonzero centers the longitude Fourier transform on the center of the remapped image", NULL},
|
162 tplarson 1.10 {ARG_INT, "ZEROMISS", "0", "zero sets any row containing a NaN to DRMS_MISSING", NULL},
|
163 tplarson 1.1 {ARG_INT, "LGAPOD", "0", "nonzero apodizes in longitude", NULL},
164 {ARG_DOUBLE, "LGAPMIN", "60.0", "start of longitude apodization, degrees", NULL},
165 {ARG_DOUBLE, "LGAPWIDT", "10.0", "width of longitude apodization, degrees", NULL},
166 // these are from jv2helio
167 {ARG_INT, "MAPMMAX", "1536", "determines mapcols", ""}, /* determines mapcols, default value is 3*512 */
168 {ARG_INT, "SINBDIVS", "512", "number of increments in sin latitude from 0 to 1", ""}, /* # of = increments in sinB from sin(0) to sin(PI/2) */
169 {ARG_FLOAT, "MAPRMAX", "0.95", "maximum image radius", ""},
|
170 tplarson 1.24 {ARG_INT, "NAN_BEYOND_RMAX", "0", "set to nonzero to set output to DRMS_MISSING outside MAPRMAX, otherwise uses 0.0 outside MAPRMAX"}, /* Non 0 sets data outside RMAX MISSING */
|
171 tplarson 1.1 {ARG_FLOAT, "MAPLGMAX", "72.0", "longitude maximum, degrees", ""}, /* degrees */
172 {ARG_FLOAT, "MAPLGMIN", "-72.0", "longitude minimum, degrees", ""},
173 {ARG_FLOAT, "MAPBMAX", "72.0", "latitude maximum, degrees, also used for minimum", ""},
174 {ARG_INT, "LGSHIFT", "0", "option for longitude shift: 0=none; 1=fixed rate; 2=nearest degree; 3=nearest tenth of a degree", ""}, /* 0=none; 1=fixed rate; 2=nearest Degree */
175 {ARG_TIME, "REF_T0", "1987.01.03_17:31:12_TAI", "reference time for computing fixed rate longitude shift", ""},
176 {ARG_FLOAT, "REF_L0", "0.0", "reference longitude for computing fixed rate longitude shift ", ""},
177 {ARG_FLOAT, "SOLAR_P", "999.0", "P-angle; if unset, taken from keywords", ""}, /* can't use D_MISSING here */
178 {ARG_FLOAT, "PSIGN", "1.0", "sign of SOLAR_P", ""}, /* Sign of P. For MWO data. */
179 {ARG_FLOAT, "PERR", "0.0", "fixed P-angle error, likely -0.22", ""}, /* Fixed P-angle error. Maybe -0.22. */
180 {ARG_FLOAT, "IERR", "0.0", "error in Carrington inclination, likely -0.10", ""}, /* Error in Carrington inclination. Maybe -0.10. */
181 {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", ""},
182 {ARG_INT, "INTERPO", "1", "option for interpolation: 0=bilinear; 1=cubic convolution", ""}, /* 2 methods - see soi_fun.h */
183 {ARG_INT, "APODIZE", "0", "option for apodization: 0=none; 1=use true solar coordinates; 2=use ideal solar coordinates (b0=0)", ""}, /* see soi_fun.h or apodize.c */
184 {ARG_FLOAT, "APINNER", "0.90", "start of apodization in fractional image radius", ""}, /* start of apodization */
185 {ARG_FLOAT, "APWIDTH", "0.05", "width of apodization in fractional image radius", ""}, /* width of apodization */
186 {ARG_INT, "APEL", "0", "set to nonzero for elliptical apodization described by APX and APY", ""}, /* do elliptical apodization described by apx and apy */
187 {ARG_FLOAT, "APX", "1.00", "divide the x position by this before applying apodization", ""}, /* divide the x position by this before applying apodization */
188 {ARG_FLOAT, "APY", "1.00", "divide the y position by this before applying apodization", ""}, /* divide the y position by this before applying apodization */
189 {ARG_INT, "VCORLEV", "2", "option for velocity correction: 0=none; 1=subtract a model of differential rotation; 2=also divide by line of sight projection factor for purely radial velocities", ""}, /* 3 levels - see soi_fun.h*/
|
190 tplarson 1.24 {ARG_INT, "MCORLEV", "0", "option for magnetic correction: 0=none; 1=line of sight; 2=radial", ""}, /* 2 levels - see soi_fun.h*/
|
191 tplarson 1.1 {ARG_INT, "MOFFSETFLAG", "0", "set to nonzero to get BFITZERO from input record and subtract from data before interpolating", ""}, /* 1=apply BFITZERO correction*/
192 {ARG_FLOAT, "OUTSCALE", "1.0", "bscale to use for output", ""}, /* scale for output */
193 {ARG_FLOAT, "OUTBIAS", "0.0", "bzero to use for output", ""}, /* bias for scaled output */
194 {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 */
195 {ARG_FLOAT, "CUBIC", "7.06E-9", "cubic distortion in fd units", ""}, /* Cubic distortion in FD units */
196 {ARG_FLOAT, "TILTALPHA", "2.59", "tilt of CCD, degrees", ""}, /* TILT of CCD in degrees */
197 {ARG_FLOAT, "TILTBETA", "56.0", "direction of CCD tilt, degrees", ""}, /* Direction of TILT in degrees */
198 {ARG_FLOAT, "TILTFEFF", "12972.629", "effective focal length", ""}, /* Effective focal length */
199 {ARG_INT, "OFLAG", "0", "set to nonzero to force reading orientation from keyword, otherwise \"SESW\" is assumed)", ""}, /* Non 0 skips checko (SESW assumed) */
200 {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*/
201 {ARG_INT, "MAXMISSVALS", "0", "if >0, this becomes threshold on MISSVALS from keyword", ""}, /* max. allowed MISSING pixels */
202 {ARG_INT, "CARRSTRETCH", "0", "set to nonzero to correct for differential rotation according to DIFROT_[ABC]", ""}, /* 0 - don't correct for diff rot, 1 - correct */
203 {ARG_FLOAT, "DIFROT_A", "13.562", "A coefficient in differential rotation adjustment (offset)", ""}, /* A coefficient in diff rot adj (offset) */
204 {ARG_FLOAT, "DIFROT_B", "-2.04", "B coefficient (to sin(lat) ^ 2)", ""}, /* B coefficient (to sin(lat) ^ 2) */
205 {ARG_FLOAT, "DIFROT_C", "-1.4875", "C coefficient (to sin(lat) ^ 4)", ""}, /* C coefficient (to sin(lat) ^ 4) */
206
207 {ARG_END}
208 };
209
210 #include "saveparm.c"
211 #include "timing.c"
|
212 tplarson 1.8 #include "set_history.c"
|
213 tplarson 1.18 #include "calversfunctions.c"
|
214 tplarson 1.1
|
215 tplarson 1.21 int SetDistort(int dist, double cubic, double alpha, double beta, double feff, LIBPROJECTION_Dist_t *dOut);
|
216 tplarson 1.19
|
217 tplarson 1.21 int obs2helio(float *V, float *U, int xpixels, int ypixels, double x0, double y0, double BZero, double P,
|
218 tplarson 1.19 double S, double rsun, double Rmax, int interpolation, int cols, int rows, double Lmin,
219 double Ldelta, double Ladjust, double sinBdelta, double smajor, double sminor, double sangle,
220 double xscale, double yscale, const char *orientation, int mag_correction, int velocity_correction,
221 double obs_vr, double obs_vw, double obs_vn, double vsign, int NaN_beyond_rmax, int carrStretch,
|
222 tplarson 1.21 const LIBPROJECTION_Dist_t *distpars, float diffrotA, float diffrotB, float diffrotC,
223 LIBPROJECTION_RotRate_t *rRates, int size);
|
224 tplarson 1.19
225 int apodize(float *data, double b0, int cols, int rows, double Lmin, double Ldelta, double sinBdelta,
226 int apodlevel, double apinner, double apwidth, int apel, double apx, double apy);
227
228 void setplm2(int lmin,int lmax,int m,long nx,int *indx,double *x,long nplm,double *plm,double *dplm);
229
|
230 tplarson 1.21 //char *getshtversion(void);
|
231 tplarson 1.19
|
232 arta 1.3 /* forward decls for fortran functions */
233 void scopy_(int *, float *, int *, float *, int *);
234 void setplm_(int *, int *, int *, int *, int *, double *, int *, double *);
235 void sgemm_(); /* give up */
236
|
237 tplarson 1.1 typedef enum
238 {
239 V2HStatus_Success,
240 V2HStatus_MissingParameter,
241 V2HStatus_IllegalTimeFormat,
242 V2HStatus_TimeConvFailed,
243 V2HStatus_Unimplemented,
244 V2HStatus_IllegalOrientation
245 } V2HStatus_t;
246
247 static void CheckO(const char *orientation, V2HStatus_t *stat)
248 {
249 /* check for legal orientation string */
250 static char o[16];
251 char *c = o;
252
253 if(!orientation)
254 {
255 *stat = V2HStatus_MissingParameter;
256 }
257 else if (4 != sscanf(orientation, "%[NS]%[EW]%[NS]%[EW]", c++, c++, c++, c++))
258 tplarson 1.1 {
259 *stat = V2HStatus_IllegalOrientation;
260 }
261 else if ((o[0] == o[2]) && (o[1] == o[3]))
262 {
263 *stat = V2HStatus_IllegalOrientation;
264 }
265 else if ((o[0] !=o [2]) && (o[1] != o[3]))
266 {
267 *stat = V2HStatus_IllegalOrientation;
268 }
269 else
270 {
271 *stat = V2HStatus_Success;
272 }
273 }
274
275 static inline void ParameterDef(int status,
276 char *pname,
277 double defaultp,
278 double *p,
|
279 tplarson 1.8 long long recnum,
|
280 tplarson 1.1 int verbflag)
281 {
282 /* logs warning and sets parameter to default value */
283 if (status != 0)
284 {
285 *p = defaultp;
286 if (verbflag)
|
287 tplarson 1.8 fprintf(stderr, "WARNING: default value %g used for %s, status = %d, recnum = %lld\n", defaultp, pname, status, recnum);
|
288 tplarson 1.1 }
289 }
290
291 #define PARAMETER_ERROR(PNAME) \
292 if (status != DRMS_SUCCESS) \
293 { \
|
294 tplarson 1.8 fprintf(stderr, "ERROR: problem with required keyword %s, status = %d, T_REC = %s, recnum = %lld, histrecnum = %lld\n", PNAME, status, trecstr, inrec->recnum, histrecnum); \
|
295 tplarson 1.1 drms_free_array(inarr); \
|
296 tplarson 1.8 return 0; \
|
297 tplarson 1.1 }
298
299
300 int DoIt(void)
301 {
302 int newstat=0;
303 DRMS_RecChunking_t chunkstat = kRecChunking_None;
304 int fetchstat = DRMS_SUCCESS;
305 int status = DRMS_SUCCESS;
|
306 tplarson 1.8 char *inrecquery = NULL;
|
307 tplarson 1.1 char *outseries = NULL;
308 char *segnamein = NULL;
309 char *segnameout = NULL;
|
310 tplarson 1.8 DRMS_RecordSet_t *inrecset = NULL;
|
311 tplarson 1.1 DRMS_Record_t *inrec = NULL;
312 DRMS_Record_t *outrec = NULL;
313 DRMS_Segment_t *segin = NULL;
314 DRMS_Segment_t *segout = NULL;
315 DRMS_Array_t *inarr = NULL;
316 DRMS_Array_t *outarr = NULL;
317 DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
318 DRMS_RecLifetime_t lifetime;
319 long long histrecnum=-1;
320 int length[2];
321
322 double wt0, wt1, wt2, wt3, wt;
323 double ut0, ut1, ut2, ut3, ut;
324 double st0, st1, st2, st3, st;
325 double ct0, ct1, ct2, ct3, ct;
326
327 TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
328 char trecstr[100], tstartstr[100];
329
|
330 tplarson 1.6 int quality;
|
331 tplarson 1.1
332 float *v2hptr, *h2mptr;
333
334 // from jqdotprod, norm changed to normflag
335 float *oddpart, *evenpart, *inptr, *workptr, *mptr, *outptr;
336 float *folded, *masks, *real4evenl, *real4oddl, *imag4evenl, *imag4oddl, *outx;
337
338 /* used for setting up plm's */
339 double *plm, *saveplm, *latgrid;
340 int *indx;
341
342 double sinBdelta;
343 double *plmptr;
344 float *maskptr;
345
346 int nsn, fournsn, snx, maxnsn;
347 int lchunksize, lchunkfirst, lchunklast, lchunk, l;
348 int lmin, lmax;
349 int msize, mx, foldedsize;
350 int maprows, nlat, imagesize, latx, poslatx, neglatx, moffset, nlatx;
351 int i, m, modd, meven;
352 tplarson 1.1 int lfirst, llast, ifirst, ilast, lstart, ldone;
353 int lfirsteven, llasteven, nevenl, lfirstodd, llastodd, noddl;
354 int fromoffset, tooffset, imageoffset;
355
356 int increment1 = 1; /* for scopy call */
357 int increment2 = 2; /* for scopy call */
358
359 /* arguments for sgemm call */
360 char transpose[] = "t";
361 char normal[] = "n";
362 float one = 1.0;
363 float zero = 0.0;
364 int normflag;
365 float cnorm; /* Constant to get proper normalization. */
366
|
367 tplarson 1.9 double tstart, tepoch, tstep, tround, cadence, nseconds, chunksecs, cadence0;
|
368 tplarson 1.1 char *ttotal, *tchunk;
369 int nrecs, irec, trecind, bc, iset, ntimechunks;
370 int *bad;
|
371 tplarson 1.6 int ndt;
|
372 tplarson 1.1
373 // from jhelio2mlat, i, lmax duplicated.
374 double mean, norm=1.0, normx;
375 int subtract_mean, normalize, cent_long, zero_miss, lgapod;
376 double lgapmin, lgapwidt, lgapmax, lon;
377
378 int row, col;
379 int lfft, mapped_lmax;
380 int mapcols2;
381 int nfft, nmean, nok, nout;
382
383 float *buf, *bp, *ip, *inp, *op, *outp, *weight, *wp;
384 float *wbuf;
385 fftwf_plan fftwp;
386
387 // from jv2helio, row, mapped_lmax, maprows, sinBdelta duplicated
388 V2HStatus_t vstat = V2HStatus_Success;
389 const char *orientationdef = "SESW ";
390 char *orientation = NULL;
391 int paramsign;
392 int longitude_shift, velocity_correction, interpolation, apodization;
393 tplarson 1.1 int mag_correction;
394 int mag_offset;
395 int sinb_divisions, mapcols, nmax, nmin;
396 int carrStretch = 0;
397 float diffrotA = 0.0;
398 float diffrotB = 0.0;
399 float diffrotC = 0.0;
400 double tobs, tmearth, tref, trefb0;
401 double smajor, sminor, sangle;
402 double xscale, yscale, imagescale;
403 int xpixels, ypixels, pixels;
404 double obs_vr, obs_vw, obs_vn;
405 double b0, bmax, bmin;
406 double longmax, longmin, longmax_adjusted, longmin_adjusted, longinterval;
407 double p0, p, rmax;
408 double ierr, perr, psign;
409 double x0, y0;
410 double obsdist, longshift, obsl0, refl0, mapl0, longrate, rtrue, rsun, S;
|
411 tplarson 1.21 double rsunDef, rsunobs;
|
412 tplarson 1.1 int obsCR;
413 int apel;
414 double apinner, apwidth, apx, apy;
415 double scale, bias;
416 double colsperdeg;
417
418 double satrot, instrot;
419 double dsignout, vsign;
420 int distsave;
421 double cubsave, tiltasave, tiltbsave, tiltfsave;
|
422 tplarson 1.21 LIBPROJECTION_Dist_t distP;
|
423 tplarson 1.1
|
424 tplarson 1.2 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
425 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
426
|
427 tplarson 1.1 wt0=getwalltime();
428 ct0=getcputime(&ut0, &st0);
429
|
430 tplarson 1.8 inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
|
431 tplarson 1.9 outseries = (char *)cmdparams_save_str(&cmdparams, "tsout", &newstat);
432 int tsflag = strcmp(kNOTSPECIFIED, outseries);
|
433 tplarson 1.8 segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
434 segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
|
435 tplarson 1.1 int seginflag = strcmp(kNOTSPECIFIED, segnamein);
436 int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
437 int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
438 int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
439 if (permflag)
440 lifetime = DRMS_PERMANENT;
441 else
442 lifetime = DRMS_TRANSIENT;
|
443 tplarson 1.8 int forceoutput = cmdparams_save_int(&cmdparams, "FORCEOUTPUT", &newstat);
444 char *tag = (char *)cmdparams_save_str(&cmdparams, "TAG", &newstat);
|
445 tplarson 1.23 char *version = (char *)cmdparams_save_str(&cmdparams, "VERSION", &newstat);
446 int verflag = strcmp(kNOTSPECIFIED, version);
|
447 tplarson 1.20 unsigned short calverkey = (unsigned short)cmdparams_save_int(&cmdparams, "CALVERKEY", &newstat);
|
448 tplarson 1.1
|
449 tplarson 1.8 char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
450 char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
|
451 tplarson 1.1
|
452 tplarson 1.8 char *v2hout = (char *)cmdparams_save_str(&cmdparams, "v2hout", &newstat);
453 char *h2mout = (char *)cmdparams_save_str(&cmdparams, "h2mout", &newstat);
|
454 tplarson 1.1 int v2hflag = strcmp(kNOTSPECIFIED, v2hout);
455 int h2mflag = strcmp(kNOTSPECIFIED, h2mout);
|
456 tplarson 1.14 int histflag = strncasecmp("none", histlinkname, 4);
|
457 tplarson 1.1
|
458 tplarson 1.9 if (!v2hflag && !h2mflag && !tsflag)
459 {
460 fprintf(stderr, "ERROR: no outputs specified.\n");
461 return 1;
462 }
463
|
464 tplarson 1.1 lmin=cmdparams_save_int(&cmdparams, "LMIN", &newstat);
465 lmax=cmdparams_save_int(&cmdparams, "LMAX", &newstat);
466 lchunksize=cmdparams_save_int(&cmdparams, "LCHUNK", &newstat);
467 normflag=cmdparams_save_int(&cmdparams, "NORM", &newstat);
468
469 tstart=cmdparams_save_time(&cmdparams, "TSTART", &newstat);
|
470 tplarson 1.14 sprint_time(tstartstr, tstart, "TAI", 0);
|
471 tplarson 1.8 ttotal=(char *)cmdparams_save_str(&cmdparams, "TTOTAL", &newstat);
|
472 tplarson 1.23 nseconds=atoinc(ttotal);
473
474 /*
475 if (strcmp(kNOTSPECIFIED, ttotal))
476 {
477 nseconds=atoinc(ttotal);
478 }
479 else
480 nseconds=0.0;
481 */
482 /*
|
483 arta 1.7 status=drms_names_parseduration(&ttotal, &nseconds, 1);
|
484 tplarson 1.16 if (status != DRMS_SUCCESS)
|
485 tplarson 1.8 {
486 fprintf(stderr, "ERROR: problem parsing TTOTAL, = %s\n", ttotal);
487 return 1;
488 }
|
489 tplarson 1.23 */
|
490 tplarson 1.8 tchunk=(char *)cmdparams_save_str(&cmdparams, "TCHUNK", &newstat);
|
491 tplarson 1.1 if (strcmp(kNOTSPECIFIED, tchunk))
492 {
|
493 tplarson 1.23 chunksecs=atoinc(tchunk);
494 /*
|
495 tplarson 1.8 status=drms_names_parseduration(&tchunk, &chunksecs, 1);
|
496 tplarson 1.16 if (status != DRMS_SUCCESS)
|
497 tplarson 1.1 newstat = newstat | CPSAVE_UNKNOWN_ERROR;
|
498 tplarson 1.23 */
|
499 tplarson 1.1 }
|
500 tplarson 1.9 else if (!tsflag)
501 {
502 fprintf(stderr, "ERROR: TCHUNK must be specified if no tsout is given.\n");
503 return 1;
504 }
|
505 tplarson 1.1 else
|
506 tplarson 1.23 chunksecs=0.0;
|
507 tplarson 1.1
508
509 subtract_mean = cmdparams_save_int(&cmdparams, "SUBMEAN", &newstat);
510 normalize = cmdparams_save_int(&cmdparams, "NORMLIZE", &newstat);
511 /* CENTLONG=1 centers the longitude Fourier transform on the center
512 of the remapped image */
513 cent_long = cmdparams_save_int(&cmdparams, "CENTLONG", &newstat);
514 /* ZEROMISS=1 sets missing data to 0,
515 ZEROMISS=0 fills the output row with missing */
516 zero_miss = cmdparams_save_int(&cmdparams, "ZEROMISS", &newstat);
517 lgapod = cmdparams_save_int(&cmdparams, "LGAPOD", &newstat);
518 lgapmin = cmdparams_save_double(&cmdparams, "LGAPMIN", &newstat);
519 lgapwidt = cmdparams_save_double(&cmdparams, "LGAPWIDT", &newstat);
520 lgapmax = lgapmin+lgapwidt;
521
522 int checko = cmdparams_save_int(&cmdparams, "OFLAG", &newstat);
|
523 tplarson 1.10 int NaN_beyond_rmax = cmdparams_save_int(&cmdparams, "NAN_BEYOND_RMAX", &newstat);
|
524 tplarson 1.1 int maxmissvals = cmdparams_save_int(&cmdparams, "MAXMISSVALS", &newstat);
|
525 tplarson 1.10 // float beyondrmax = cmdparams_save_float(&cmdparams, "BEYONDRMAX", &newstat);
|
526 tplarson 1.1
527 carrStretch = cmdparams_save_int(&cmdparams, "CARRSTRETCH", &newstat);
528 diffrotA = cmdparams_save_float(&cmdparams, "DIFROT_A", &newstat);
529 diffrotB = cmdparams_save_float(&cmdparams, "DIFROT_B", &newstat);
530 diffrotC = cmdparams_save_float(&cmdparams, "DIFROT_C", &newstat);
531
532 longrate = 360.0 / TCARR - 360.0 / DAYSINYEAR; // degrees per day
533 longrate /= SECSINDAY; // degrees per sec
534
535 apodization = cmdparams_save_int(&cmdparams, "APODIZE", &newstat);
536 apinner = cmdparams_save_double(&cmdparams, "APINNER", &newstat);
537 apwidth = cmdparams_save_double(&cmdparams, "APWIDTH", &newstat);
538 apel = cmdparams_save_int(&cmdparams, "APEL", &newstat);
539 apx = cmdparams_save_double(&cmdparams, "APX", &newstat);
540 apy = cmdparams_save_double(&cmdparams, "APY", &newstat);
541 longitude_shift = cmdparams_save_int(&cmdparams, "LGSHIFT", &newstat);
542 mag_correction = cmdparams_save_int(&cmdparams, "MCORLEV", &newstat);
543 mag_offset = cmdparams_save_int(&cmdparams, "MOFFSETFLAG", &newstat);
544 velocity_correction = cmdparams_save_int(&cmdparams, "VCORLEV", &newstat);
545 interpolation = cmdparams_save_int(&cmdparams, "INTERPO", &newstat);
546 paramsign = cmdparams_save_int(&cmdparams, "DATASIGN", &newstat);
547 tplarson 1.1 rmax = cmdparams_save_double(&cmdparams, "MAPRMAX", &newstat);
548 refl0 = cmdparams_save_double(&cmdparams, "REF_L0", &newstat);
549
550 distsave = cmdparams_save_int(&cmdparams, "DISTORT", &newstat);
551 cubsave = cmdparams_save_double(&cmdparams, "CUBIC", &newstat);
552 tiltasave = cmdparams_save_double(&cmdparams, "TILTALPHA", &newstat);
553 tiltbsave = cmdparams_save_double(&cmdparams, "TILTBETA", &newstat);
554 tiltfsave = cmdparams_save_double(&cmdparams, "TILTFEFF", &newstat);
555
556 scale = cmdparams_save_double(&cmdparams, "OUTSCALE", &newstat);
557 bias = cmdparams_save_double(&cmdparams, "OUTBIAS", &newstat);
558 p0 = cmdparams_save_double(&cmdparams, "SOLAR_P", &newstat);
559 psign = cmdparams_save_double(&cmdparams, "PSIGN", &newstat);
560 perr = cmdparams_save_double(&cmdparams, "PERR", &newstat);
561 ierr = cmdparams_save_double(&cmdparams, "IERR", &newstat);
562 trefb0 = cmdparams_save_time(&cmdparams, "REF_TB0", &newstat);
563
564 SetDistort(distsave, cubsave, tiltasave, tiltbsave, tiltfsave, &distP);
565
566 tref = cmdparams_save_time(&cmdparams, "REF_T0", &newstat);
567
568 tplarson 1.1 // determine mapcols and adjust longmin and longmax */
569 mapped_lmax = cmdparams_save_int(&cmdparams, "MAPMMAX", &newstat);
570 longmax = cmdparams_save_double(&cmdparams, "MAPLGMAX", &newstat); /* degrees */
571 longmin = cmdparams_save_double(&cmdparams, "MAPLGMIN", &newstat); /* degrees */
572 longinterval = (180.0) / mapped_lmax; /* degrees */
573
574 // This does not always handle the case where 1/longinterval is an integer correctly.
575
576 // the next two statement do nothing, right?
577 // why do nmin and max keep getting set with different RHSs?
578 nmin = (int)(longmin / longinterval); // round towards 0
579 nmax = (int)(longmax / longinterval); // round towards 0
580 colsperdeg = mapped_lmax / 180.0;
581 nmin = (int)(longmin * colsperdeg); // round towards 0
582 nmax = (int)(longmax * colsperdeg); // round towards 0
583 mapcols = nmax - nmin + 1;
584 longmin_adjusted = nmin * longinterval;
585 longmax_adjusted = nmax * longinterval;
586
587 // determine maprows, bmax, bmin, and sinBdelta
588 sinb_divisions = cmdparams_save_int(&cmdparams, "SINBDIVS", &newstat);
589 tplarson 1.1 sinBdelta = 1.0/sinb_divisions;
590 bmax = cmdparams_save_double(&cmdparams, "MAPBMAX", &newstat); // degrees
591 bmin = -bmax;
592 nmax = (int)(sin(RADSINDEG*bmax)*sinb_divisions); // round towards 0
593 maprows = 2*nmax;
594
|
595 tplarson 1.8 if (normflag == 0)
596 cnorm=1.0;
597 else
598 cnorm = sqrt(2.)*sinBdelta;
599
|
600 tplarson 1.1 if (lmax > mapped_lmax || lmin > lmax)
601 {
602 fprintf(stderr, "ERROR: must have MAPMMAX >= LMAX >= LMIN, MAPMMAX = %d, LMAX= %d, LMIN = %d\n", mapped_lmax, lmax, lmin);
603 return 1;
604 }
605
|
606 tplarson 1.8 if (newstat)
607 {
608 fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
609 cpsave_decode_error(newstat);
610 return 1;
611 }
612 else if (savestrlen != strlen(savestr))
613 {
614 fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
615 return 1;
616 }
617
|
618 tplarson 1.9 DRMS_Record_t *tempoutrec;
619 DRMS_Link_t *histlink;
620 int itest;
|
621 tplarson 1.21 /*
622 // cvsinfo used to be passed in the call to set_history. now this information is encoded in CVSTAG, which is defined by a compiler flag in the make.
|
623 tplarson 1.17 char *cvsinfo;
624 cvsinfo = (char *)malloc(1024);
|
625 tplarson 1.24 strcpy(cvsinfo,"$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jv2ts.c,v 1.23 2013/07/02 20:32:03 tplarson Exp $");
|
626 tplarson 1.17 strcat(cvsinfo,"\n");
627 strcat(cvsinfo,getshtversion());
|
628 tplarson 1.21 */
|
629 tplarson 1.9 // assume all output dataseries link to the same dataseries for HISTORY
630 if (tsflag)
631 {
632 tempoutrec = drms_create_record(drms_env, outseries, DRMS_TRANSIENT, &status);
633 if (status != DRMS_SUCCESS)
634 {
635 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
636 return 1;
637 }
|
638 tplarson 1.1
|
639 tplarson 1.14 DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, "MAPMMAX");
640 if (outkeytest != NULL && outkeytest->info->recscope == 1)
641 {
642 int mapmmaxout=drms_getkey_int(tempoutrec, "MAPMMAX", &status);
643 if (mapmmaxout != mapped_lmax)
644 {
645 fprintf(stderr,"ERROR: output MAPMMAX=%d does not match input parameter MAPMMAX=%d, status = %d\n", mapmmaxout, mapped_lmax, status);
646 return 1;
647 }
648 }
649
650 outkeytest = hcon_lookup_lower(&tempoutrec->keywords, "SINBDIVS");
651 if (outkeytest != NULL && outkeytest->info->recscope == 1)
|
652 tplarson 1.9 {
|
653 tplarson 1.14 int sinbdivsout=drms_getkey_int(tempoutrec, "SINBDIVS", &status);
654 if (sinbdivsout != sinb_divisions)
|
655 tplarson 1.9 {
|
656 tplarson 1.14 fprintf(stderr,"ERROR: output SINBDIVS=%d does not match input parameter SINBDIVS=%d, status = %d\n", sinbdivsout, sinb_divisions, status);
657 return 1;
|
658 tplarson 1.9 }
659 }
|
660 tplarson 1.14
661 // set up ancillary dataseries for processing metadata
662 if (histflag)
|
663 tplarson 1.9 {
|
664 tplarson 1.14 histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
665 if (histlink != NULL)
666 {
|
667 tplarson 1.21 histrecnum=set_history(histlink);
|
668 tplarson 1.14 if (histrecnum < 0)
669 {
670 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
671 return 1;
672 }
673 }
674 else
675 {
676 fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
677 }
|
678 tplarson 1.9 }
|
679 tplarson 1.1
|
680 tplarson 1.9 // these must be present in the output dataseries and variable, not links or constants
|
681 tplarson 1.14 char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"};
|
682 tplarson 1.9 for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
683 {
684 DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
685 if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1)
686 {
687 fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
688 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
689 return 1;
690 }
691 }
|
692 tplarson 1.1
|
693 tplarson 1.9 cadence0=drms_getkey_float(tempoutrec, "T_STEP", &status);
694 tepoch=drms_getkey_time(tempoutrec, "T_START_epoch", &status);
695 tstep=drms_getkey_float(tempoutrec, "T_START_step", &status);
696 tround=drms_getkey_float(tempoutrec, "T_START_round", &status);
697 if (fmod(tstart-tepoch,tstep) > tround/2)
698 {
|
699 tplarson 1.14 sprint_time(trecstr, tepoch, "TAI", 0);
700 fprintf(stderr, "ERROR: output dataseries seems incompatible with input parameters (tstep must divide tstart-tepoch): TSTART = %s, T_START_epoch = %s, tstep = %f\n",
701 tstartstr, trecstr, tstep);
|
702 tplarson 1.9 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
703 return 1;
704 }
705 if (chunksecs == 0.0)
706 chunksecs = tstep;
707 else if (fmod(chunksecs,tstep))
|
708 tplarson 1.1 {
|
709 tplarson 1.9 fprintf(stderr, "ERROR: output dataseries seems incompatible with input parameters (tstep must divide chunksecs): chunksecs = %f, tstep = %f\n", chunksecs, tstep);
|
710 tplarson 1.4 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
|
711 tplarson 1.1 return 1;
712 }
|
713 tplarson 1.9
|
714 tplarson 1.14 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
|
715 tplarson 1.1 }
|
716 tplarson 1.14
717 if (v2hflag)
|
718 tplarson 1.1 {
|
719 tplarson 1.9 tempoutrec = drms_create_record(drms_env, v2hout, DRMS_TRANSIENT, &status);
720 if (status != DRMS_SUCCESS)
721 {
722 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", v2hout, status);
723 return 1;
724 }
725
726 // set up ancillary dataseries for processing metadata
|
727 tplarson 1.14 if (histflag && histrecnum < 0)
|
728 tplarson 1.9 {
|
729 tplarson 1.14 histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
730 if (histlink != NULL)
|
731 tplarson 1.9 {
|
732 tplarson 1.21 histrecnum=set_history(histlink);
|
733 tplarson 1.14 if (histrecnum < 0)
734 {
735 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
736 return 1;
737 }
738 }
739 else
740 {
741 fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
|
742 tplarson 1.9 }
743 }
|
744 tplarson 1.1
745 // these must be present in the output dataseries and variable, not links or constants
|
746 tplarson 1.9 char *outchecklist[] = {"T_REC", "QUALITY", "CRPIX1", "CRVAL1", "CDELT1", "CRPIX2", "CROTA2", "CDELT2" };
747 for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
|
748 tplarson 1.1 {
|
749 tplarson 1.9 DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
750 if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1)
751 {
752 fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
753 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
754 return 1;
755 }
|
756 tplarson 1.1 }
757
|
758 tplarson 1.14 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
|
759 tplarson 1.4 }
|
760 tplarson 1.14
761 if (h2mflag)
|
762 tplarson 1.4 {
|
763 tplarson 1.9 tempoutrec = drms_create_record(drms_env, h2mout, DRMS_TRANSIENT, &status);
764 if (status != DRMS_SUCCESS)
765 {
766 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", h2mout, status);
767 return 1;
768 }
769
770 // set up ancillary dataseries for processing metadata
|
771 tplarson 1.14 if (histflag && histrecnum < 0)
|
772 tplarson 1.9 {
|
773 tplarson 1.14 histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
774 if (histlink != NULL)
|
775 tplarson 1.9 {
|
776 tplarson 1.21 histrecnum=set_history(histlink);
|
777 tplarson 1.14 if (histrecnum < 0)
778 {
779 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
780 return 1;
781 }
782 }
783 else
784 {
785 fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
|
786 tplarson 1.9 }
787 }
788
789 // these must be present in the output dataseries and variable, not links or constants
790 char *outchecklist[] = {"T_REC", "QUALITY", "CRPIX1", "CDELT1", "CRPIX2", "CDELT2" };
791 for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
792 {
793 DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
794 if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1)
795 {
796 fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
797 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
798 return 1;
799 }
800 }
801
|
802 tplarson 1.14 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
|
803 tplarson 1.4 }
|
804 tplarson 1.14
|
805 tplarson 1.9
|
806 tplarson 1.4 if (fmod(nseconds,chunksecs) != 0.0)
807 {
|
808 tplarson 1.14 fprintf(stderr, "ERROR: input parameters seem incompatible (chunksecs must divide totalsecs): totalsecs = %f, chunksecs = %f\n", nseconds, chunksecs);
|
809 tplarson 1.8 return 1;
|
810 tplarson 1.4 }
811 ntimechunks=nseconds/chunksecs;
812
|
813 tplarson 1.8 inrecset = drms_open_recordset(drms_env, inrecquery, &status);
814 if (status != DRMS_SUCCESS || inrecset == NULL)
815 {
816 fprintf(stderr, "ERROR: problem opening input recordset: status = %d\n", status);
817 return 1;
818 }
|
819 tplarson 1.1
|
820 tplarson 1.20 int nrecsin = drms_count_records(drms_env, inrecquery, &status);
821 if (status != DRMS_SUCCESS)
822 {
823 fprintf(stderr, "ERROR: problem counting input records: status = %d, nrecs = %d\n", status, nrecsin);
824 drms_close_records(inrecset, DRMS_FREE_RECORD);
825 return 1;
826 }
827 if (nrecsin == 0)
828 {
829 fprintf(stderr, "ERROR: input recordset contains no records. if such was intended use jretile instead.\n");
830 drms_close_records(inrecset, DRMS_FREE_RECORD);
831 return 1;
832 }
833
834 //the above replaces the following. drms_open_recordset() no longer fills in the number of records.
835 /*
|
836 tplarson 1.8 if (inrecset->n == 0)
|
837 tplarson 1.1 {
|
838 tplarson 1.8 fprintf(stderr, "ERROR: input recordset contains no records. if such was intended use jretile instead.\n");
839 drms_close_records(inrecset, DRMS_FREE_RECORD);
840 return 1;
|
841 tplarson 1.1 }
|
842 tplarson 1.20 */
|
843 tplarson 1.1
844 if (verbflag)
|
845 tplarson 1.20 printf("input recordset opened, nrecs = %d\n", nrecsin);
|
846 tplarson 1.1
|
847 tplarson 1.8 inrec = drms_recordset_fetchnext(drms_env, inrecset, &fetchstat, &chunkstat, NULL);
|
848 tplarson 1.1
|
849 tplarson 1.8 // these must be present in the input dataseries
|
850 tplarson 1.12 char *inchecklist[] = {"T_REC", "QUALITY", "T_OBS", "CRLT_OBS", "CRLN_OBS", "CADENCE",
|
851 tplarson 1.14 // "SAT_ROT", "INST_ROT", "IM_SCALE",
|
852 tplarson 1.6 "CDELT1", "CDELT2"};
|
853 tplarson 1.1
|
854 tplarson 1.21 DRMS_Keyword_t *inkeytest;
|
855 tplarson 1.1 for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
856 {
|
857 tplarson 1.21 inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
|
858 tplarson 1.8 if (inkeytest == NULL)
|
859 tplarson 1.1 {
860 fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
|
861 tplarson 1.8 drms_close_records(inrecset, DRMS_FREE_RECORD);
862 return 1;
|
863 tplarson 1.1 }
864 }
865
|
866 tplarson 1.21 int readrsunref=0;
867 double rsunref;
868 inkeytest = hcon_lookup_lower(&inrec->keywords, "RSUN_REF");
869 if (inkeytest == NULL)
870 rtrue = RTRUE/AU;
871 else if (inkeytest->info->recscope == 1)
872 {
873 rsunref = drms_getkey_double(inrec, "RSUN_REF", &status);
874 ParameterDef(status, "RSUN_REF", RTRUE, &rsunref, inrec->recnum, 1);
875 rtrue=rsunref/AU;
876 }
877 else
878 readrsunref=1;
|
879 tplarson 1.1
880 trec = drms_getkey_time(inrec, "T_REC", &status);
|
881 tplarson 1.8 if (status != DRMS_SUCCESS)
|
882 tplarson 1.1 {
883 fprintf(stderr, "ERROR: problem with required parameter T_REC: status = %d, recnum = %lld\n", status, inrec->recnum);
|
884 tplarson 1.8 drms_close_records(inrecset, DRMS_FREE_RECORD);
885 return 1;
|
886 tplarson 1.1 }
887 sprint_time(trecstr, trec, "TAI", 0);
888
889 cadence=drms_getkey_float(inrec, "CADENCE", &status);
|
890 tplarson 1.8 if (status != DRMS_SUCCESS)
|
891 tplarson 1.1 {
|
892 tplarson 1.14 fprintf(stderr, "ERROR: problem with required parameter CADENCE: status = %d, T_REC = %s, recnum = %lld\n", status, trecstr, inrec->recnum);
|
893 tplarson 1.8 drms_close_records(inrecset, DRMS_FREE_RECORD);
894 return 1;
|
895 tplarson 1.1 }
|
896 tplarson 1.11
|
897 tplarson 1.9 if (!forceoutput)
898 {
|
899 tplarson 1.20 if (nrecsin != nseconds/cadence)
|
900 tplarson 1.9 {
901 fprintf(stderr, "ERROR: input recordset does not contain a record for every slot.\n");
902 drms_close_records(inrecset, DRMS_FREE_RECORD);
903 return 1;
904 }
905 }
|
906 tplarson 1.6
|
907 tplarson 1.9 if (tsflag && cadence != cadence0)
|
908 tplarson 1.8 {
|
909 tplarson 1.14 fprintf(stderr, "ERROR: input CADENCE does not match output T_STEP: status = %d, T_REC = %s, recnum = %lld\n", status, trecstr, inrec->recnum);
|
910 tplarson 1.8 drms_close_records(inrecset, DRMS_FREE_RECORD);
911 return 1;
912 }
|
913 tplarson 1.1
914 nrecs=chunksecs/cadence;
915 maxnsn=nsn=nrecs;
916 fournsn=4*nsn;
917
918 if (verbflag)
919 printf("ntimechunks = %d, recs per timechunk = %d\n", ntimechunks, nrecs);
920
921
922 if (trec >= tstart + nseconds)
923 {
|
924 tplarson 1.4 fprintf(stderr, "ERROR: no records processed: first input record is after last output record: T_REC = %s\n", trecstr);
|
925 tplarson 1.8 drms_close_records(inrecset, DRMS_FREE_RECORD);
926 return 1;
|
927 tplarson 1.1 }
928
929 while (trec < tstart && chunkstat != kRecChunking_NoMoreRecs)
930 {
|
931 tplarson 1.14 fprintf(stderr, "WARNING: input record will not be included in output: T_REC = %s, TSTART = %s \n", trecstr, tstartstr);
|
932 tplarson 1.8 inrec = drms_recordset_fetchnext(drms_env, inrecset, &fetchstat, &chunkstat, NULL);
933 if (inrec != NULL)
|
934 tplarson 1.1 {
935 trec = drms_getkey_time(inrec, "T_REC", &status);
936 sprint_time(trecstr, trec, "TAI", 0);
937 }
938 }
939 if (chunkstat == kRecChunking_NoMoreRecs)
940 {
|
941 tplarson 1.4 fprintf(stderr,"ERROR: no records processed: last input record is before first output record: T_REC = %s\n", trecstr);
|
942 tplarson 1.8 drms_close_records(inrecset, DRMS_FREE_RECORD);
943 return 1;
|
944 tplarson 1.1 }
945
946 msize = lmax+1;
947 if (lchunksize == 0) lchunksize = msize;
948
949 nlat = maprows/2;
950 imagesize = maprows*2*msize; /* out image could be smaller than in */
951
952 nout = 2 * (lmax + 1);
953
954 lfft = 2 * mapped_lmax;
955 nfft = lfft + 2;
956 mapcols2 = mapcols/2;
957
958 /* Let's try this since SGI's like odd leading dimensions of the first
959 array in sgemm */
|
960 tplarson 1.14 // nlatx=2*(nlat/2)+1;
|
961 tplarson 1.1
962 /* make nlatx divisible by 4 on linux systems */
|
963 tplarson 1.14 //#ifdef __linux__
|
964 tplarson 1.1 if (nlat % 4)
965 nlatx=4*(nlat/4+1);
966 else
967 nlatx=nlat;
|
968 tplarson 1.14 //#endif
|
969 tplarson 1.1
|
970 tplarson 1.14 DRMS_RecordSet_t *outrecset, *v2hrecset, *h2mrecset;
|
971 tplarson 1.9 if (tsflag)
972 {
973 real4evenl = (float *)(malloc (nlatx*maxnsn*sizeof(float)));
974 real4oddl = (float *)(malloc (nlatx*maxnsn*sizeof(float)));
975 imag4evenl = (float *)(malloc (nlatx*maxnsn*sizeof(float)));
976 imag4oddl = (float *)(malloc (nlatx*maxnsn*sizeof(float)));
977 outx = (float *)(malloc (maxnsn*2*lchunksize*sizeof(float)));
978
979 plm = (double *)(malloc ((lmax+1)*nlat*sizeof(double)));
980 saveplm = (double *)(malloc ((lmax+1)*nlat*2*sizeof(double)));
981 latgrid = (double *)(malloc (nlat*sizeof(double)));
982 for (i=0; i<nlat; i++) latgrid[i] = (i+0.5)*sinBdelta;
983
984 indx = (int *)(malloc ((lmax+1)*sizeof(int)));
985 for (l=0; l<=lmax; l++) indx[l]=l;
986
987 masks = (float *)(malloc (nlat*lchunksize*sizeof(float)));
988 foldedsize = 4*nlat*(lmax+1)*maxnsn;
989 folded = (float *)(malloc (foldedsize*sizeof(float)));
990 oddpart = (float *)(malloc (nlat*sizeof(float)));
991 evenpart = (float *)(malloc (nlat*sizeof(float)));
|
992 tplarson 1.14
993 lchunkfirst = lmin/lchunksize;
994 lchunklast = lmax/lchunksize;
995 int nlchunks = (lchunklast - lchunkfirst) + 1;
996 int nrecsout = nlchunks*ntimechunks;
997 outrecset = drms_create_records(drms_env, nrecsout, outseries, lifetime, &status);
998 if (status != DRMS_SUCCESS || outrecset == NULL)
999 {
1000 fprintf(stderr,"ERROR: unable to create records in output dataseries %s, status = %d\n", outseries, status);
1001 drms_close_records(inrecset, DRMS_FREE_RECORD);
1002 return 1;
1003 }
1004
1005 }
1006
1007 if (v2hflag)
1008 {
|
1009 tplarson 1.20 v2hrecset = drms_create_records(drms_env, nrecsin, v2hout, lifetime, &status);
|
1010 tplarson 1.14 if (status != DRMS_SUCCESS || v2hrecset == NULL)
1011 {
1012 fprintf(stderr,"ERROR: unable to create records in output dataseries %s, status = %d\n", v2hout, status);
1013 drms_close_records(inrecset, DRMS_FREE_RECORD);
1014 return 1;
1015 }
1016 }
1017
1018 if (h2mflag)
1019 {
|
1020 tplarson 1.20 h2mrecset = drms_create_records(drms_env, nrecsin, h2mout, lifetime, &status);
|
1021 tplarson 1.14 if (status != DRMS_SUCCESS || h2mrecset == NULL)
1022 {
1023 fprintf(stderr,"ERROR: unable to create records in output dataseries %s, status = %d\n", h2mout, status);
1024 drms_close_records(inrecset, DRMS_FREE_RECORD);
1025 return 1;
1026 }
|
1027 tplarson 1.9 }
|
1028 tplarson 1.1
|
1029 tplarson 1.14
|
1030 tplarson 1.1 bad = (int *)(malloc (nrecs*sizeof(int)));
1031 v2hptr = (float *)(malloc(maprows*mapcols*sizeof(float)));
1032 h2mptr = (float *)(malloc(maprows*nout*sizeof(float)));
1033
1034 /* get working buffer */
1035 buf = (float *)malloc(nfft * sizeof(float));
1036
1037 wbuf = (float *)malloc(nfft * sizeof(float));
1038 fftwp = fftwf_plan_r2r_1d(lfft, buf, wbuf, FFTW_R2HC, FFTW_ESTIMATE);
1039
1040 /* get weight array for apodizing */
1041 weight = (float *)malloc(nfft * sizeof(float));
1042
1043 wp = weight;
1044 for (col=0; col<mapcols; col++)
1045 {
1046 if (lgapod)
1047 {
1048 lon=abs(col-mapcols2)*360.0/(2*mapped_lmax);
1049 if (lon < lgapmin)
1050 *wp++=1.0;
1051 tplarson 1.1 else if (lon < lgapmax)
1052 *wp++=0.5+0.5*cos(PI*(lon-lgapmin)/lgapwidt);
1053 else
1054 *wp++=0.0;
1055 }
1056 else
1057 *wp++ = 1.0;
1058 }
1059
|
1060 tplarson 1.8 status=drms_stage_records(inrecset, 1, 0);
1061 if (status != DRMS_SUCCESS)
1062 {
1063 fprintf(stderr, "ERROR: drms_stage_records returned status = %d\n", status);
1064 return 1;
1065 }
|
1066 tplarson 1.1
|
1067 tplarson 1.20 unsigned long long calversout, calvers;
|
1068 tplarson 1.15 int calversunset=1;
|
1069 tplarson 1.20 unsigned int nybblearrout[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1070 int fixflagarr[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1071 for (i=0;i<16;i++)
1072 {
1073 if (getbits(calverkey,i,1))
1074 fixflagarr[i]=1;
1075 }
1076
1077 int mixflag=0;
|
1078 tplarson 1.15
|
1079 tplarson 1.4 int nrecords=0;
1080 int nsegments=0;
1081 int error=0;
1082 int nodata=0;
|
1083 tplarson 1.14 int irecout=0;
1084 int iv2hrec=0;
1085 int ih2mrec=0;
|
1086 tplarson 1.1 for (iset=0; iset < ntimechunks && chunkstat != kRecChunking_NoMoreRecs; iset++)
1087 {
1088 sprint_time(tstartstr, tstart, "TAI", 0);
1089 if (verbflag)
1090 {
1091 wt1=getwalltime();
1092 ct1=getcputime(&ut1, &st1);
1093 printf("processing timechunk %d, tstart = %s\n", iset, tstartstr);
1094 }
1095
1096 if (trec >= tstart+chunksecs)
1097 {
|
1098 tplarson 1.8 if (forceoutput)
1099 {
1100 nodata=1;
1101 goto skip_norecs;
1102 }
1103 else
1104 {
1105 fprintf(stderr, "ERROR: no data for timechunk beginning at %s\n", tstartstr);
1106 error++;
1107 tstart+=chunksecs;
1108 continue;
1109 }
|
1110 tplarson 1.1 }
1111
|
1112 tplarson 1.4 while (trec < tstart && chunkstat != kRecChunking_NoMoreRecs)
1113 {
|
1114 tplarson 1.8 inrec = drms_recordset_fetchnext(drms_env, inrecset, &fetchstat, &chunkstat, NULL);
|
1115 tplarson 1.6 if (inrec != NULL)
|
1116 tplarson 1.4 {
1117 trec = drms_getkey_time(inrec, "T_REC", &status);
1118 sprint_time(trecstr, trec, "TAI", 0);
1119 }
1120 }
1121
|
1122 tplarson 1.1 bc=0;
|
1123 tplarson 1.4 nodata=0;
|
1124 tplarson 1.20 mixflag=0;
1125 calversunset=1;
|
1126 tplarson 1.8 trecind=(trec-tstart+cadence/2)/cadence;
|
1127 tplarson 1.1
|
1128 tplarson 1.8 for (irec=0; irec < nrecs && chunkstat != kRecChunking_NoMoreRecs; irec++)
|
1129 tplarson 1.1 {
|
1130 tplarson 1.8
1131 if (trecind > irec)
1132 //some inputs were missing
|
1133 tplarson 1.1 {
|
1134 tplarson 1.8 if (forceoutput)
1135 {
1136 bad[bc++]=irec;
1137 continue;
1138 }
1139 else
1140 {
1141 fprintf(stderr, "ERROR: some input records missing, T_START = %s, T_REC = %s, irec = %d\n", tstartstr, trecstr, irec);
1142 error++;
1143 goto continue_outer_loop;
1144 }
1145 }
1146
1147 while (trecind < irec)
1148 //T_REC is duplicated in input
1149 {
1150 if (forceoutput)
1151 {
1152 inrec = drms_recordset_fetchnext(drms_env, inrecset, &fetchstat, &chunkstat, NULL);
1153 if (inrec != NULL)
1154 {
1155 tplarson 1.8 trec = drms_getkey_time(inrec, "T_REC", &status);
1156 sprint_time(trecstr, trec, "TAI", 0);
1157 trecind=(trec-tstart+cadence/2)/cadence;
1158 }
1159 }
1160 else
1161 {
1162 fprintf(stderr, "ERROR: some input records have duplicate T_REC, T_START = %s, T_REC = %s, irec = %d\n", tstartstr, trecstr, irec);
1163 error++;
1164 goto continue_outer_loop;
1165 }
|
1166 tplarson 1.1 }
1167
1168 if (verbflag > 1)
1169 {
1170 wt2=getwalltime();
1171 ct2=getcputime(&ut2, &st2);
1172 printf(" processing record %d\n", irec);
1173 }
1174
1175 quality=drms_getkey_int(inrec, "QUALITY", &status);
1176 if (status != DRMS_SUCCESS || (quality & QUAL_NODATA)) //may want stricter test on quality here
1177 {
1178 bad[bc++]=irec;
|
1179 tplarson 1.18 if (verbflag > 2)
|
1180 tplarson 1.8 fprintf(stderr, "SKIP: image rejected based on quality: T_REC = %s, quality = %08x\n", trecstr, quality);
|
1181 tplarson 1.1 goto skip;
1182 }
1183
|
1184 tplarson 1.20 if (tsflag)
|
1185 tplarson 1.15 {
|
1186 tplarson 1.20 if (calversunset)
1187 {
1188 calversout=drms_getkey_longlong(inrec, "CALVER64", &status);
1189 if (status != DRMS_SUCCESS)
1190 calversout = 0;
1191 else
1192 calversout = fixcalver64(calversout);
1193 calversunset=0;
1194
1195 for (i=0;i<16;i++)
1196 nybblearrout[i]=getbits(calversout,4*i+3,4);
1197
1198 }
1199
1200 calvers=drms_getkey_longlong(inrec, "CALVER64", &status);
|
1201 tplarson 1.15 if (status != DRMS_SUCCESS)
|
1202 tplarson 1.20 calvers = 0;
|
1203 tplarson 1.18 else
|
1204 tplarson 1.20 calvers = fixcalver64(calvers);
|
1205 tplarson 1.15
|
1206 tplarson 1.20 for (i=0;i<16;i++)
1207 {
1208 int nybble=getbits(calvers,4*i+3,4);
1209 if (fixflagarr[i])
1210 {
1211 if (nybble != nybblearrout[i])
1212 {
1213 fprintf(stderr, "ERROR: input data has mixed values for field %d of CALVER64: %d and %d, recnum = %lld, histrecnum = %lld\n", i, nybblearrout[i], nybble, inrec->recnum, histrecnum);
1214 error++;
1215 goto continue_outer_loop;
1216 }
1217 }
1218 else
1219 {
1220 if (nybble < nybblearrout[i])
1221 nybblearrout[i]=nybble;
1222 }
1223 }
|
1224 tplarson 1.15
|
1225 tplarson 1.20 if (!mixflag && calvers != calversout)
1226 mixflag=1;
|
1227 tplarson 1.15 }
1228
|
1229 tplarson 1.20
|
1230 tplarson 1.1 if (seginflag)
1231 segin = drms_segment_lookup(inrec, segnamein);
1232 else
1233 segin = drms_segment_lookupnum(inrec, 0);
|
1234 tplarson 1.16 if (segin != NULL)
|
1235 tplarson 1.1 inarr = drms_segment_read(segin, usetype, &status);
1236
|
1237 tplarson 1.8 if (segin == NULL || inarr == NULL || status != DRMS_SUCCESS)
|
1238 tplarson 1.1 {
|
1239 tplarson 1.8 fprintf(stderr, "ERROR: problem with input segment or array: status = %d, T_REC = %s, recnum = %lld, histrecnum = %lld\n", status, trecstr, inrec->recnum, histrecnum);
1240 return 0;
|
1241 tplarson 1.1 }
1242
1243 if (maxmissvals > 0)
1244 {
1245 int missvals = drms_getkey_int(inrec, "MISSVALS", &status);
1246 PARAMETER_ERROR("MISSVALS")
1247 if (missvals > maxmissvals)
1248 {
1249 bad[bc++]=irec;
|
1250 tplarson 1.8 if (verbflag > 1)
1251 fprintf(stderr, "SKIP: %d pixels MISSING, max allowed is %d: T_REC = %s, recnum = %lld\n", missvals, maxmissvals, trecstr, inrec->recnum);
|
1252 tplarson 1.1 drms_free_array(inarr);
1253 goto skip;
1254 }
1255 }
1256
1257 tobs = drms_getkey_time(inrec, "T_OBS", &status);
1258 PARAMETER_ERROR("T_OBS")
1259
1260 // MDI keyword was OBS_B0
1261 b0 = drms_getkey_double(inrec, "CRLT_OBS", &status);
1262 PARAMETER_ERROR("CRLT_OBS")
1263
1264 // MDI keyword was OBS_L0
1265 obsl0 = drms_getkey_double(inrec, "CRLN_OBS", &status);
1266 PARAMETER_ERROR("CRLN_OBS")
1267
1268 if (p0 == 999.0)
1269 {
1270 // MDI keyword was SOLAR_P = -(SAT_ROT + INST_ROT)
|
1271 tplarson 1.6 /*
|
1272 tplarson 1.1 satrot = drms_getkey_double(inrec, "SAT_ROT", &status);
1273 PARAMETER_ERROR("SAT_ROT")
1274 instrot = drms_getkey_double(inrec, "INST_ROT", &status);
1275 PARAMETER_ERROR("INST_ROT")
1276 p=-(satrot+instrot);
|
1277 tplarson 1.6 */
1278 double crota = drms_getkey_double(inrec, "CROTA2", &status);
1279 PARAMETER_ERROR("CROTA2")
1280 p=-crota;
|
1281 tplarson 1.1 }
1282 else
1283 {
1284 p = p0;
1285 }
1286
1287 p = psign * p ;
1288 b0 = b0 + ierr * sin((tobs - trefb0) / 31557600. * 2 * PI);
1289 p = p + perr - ierr * cos((tobs - trefb0) / 31557600. * 2 * PI);
1290
1291 // S_MAJOR, S_MINOR, S_ANGLE, X_SCALE, Y_SCALE were MDI keywords, not used for HMI, but still necessary for GONG data
1292 smajor = drms_getkey_double(inrec, "S_MAJOR", &status);
|
1293 tplarson 1.8 ParameterDef(status, "S_MAJOR", 1.0, &smajor, inrec->recnum, 0);
|
1294 tplarson 1.1
1295 sminor = drms_getkey_double(inrec, "S_MINOR", &status);
|
1296 tplarson 1.8 ParameterDef(status, "S_MINOR", 1.0, &sminor, inrec->recnum, 0);
|
1297 tplarson 1.1
1298 sangle = drms_getkey_double(inrec, "S_ANGLE", &status);
|
1299 tplarson 1.8 ParameterDef(status, "S_ANGLE", 0.0, &sangle, inrec->recnum, 0);
|
1300 tplarson 1.1
|
1301 tplarson 1.21 /*
1302 our calculation of CDELTi does not follow WCS conventions. it should be the plate scale
1303 at the reference pixel (disk center), but instead we use the average between center and limb.
1304 this is taken into account in the calculation of rsun below.
1305 */
|
1306 tplarson 1.1 xscale = drms_getkey_double(inrec, "CDELT1", &status);
1307 PARAMETER_ERROR("CDELT1")
1308 yscale = drms_getkey_double(inrec, "CDELT2", &status);
1309 PARAMETER_ERROR("CDELT2")
1310
1311 // use xscale and yscale for the following check, then set to 1.0 for the call to obs2helio
1312 if (xscale != yscale)
1313 {
|
1314 tplarson 1.8 fprintf(stderr, "ERROR: CDELT1 != CDELT2 not supported, CDELT1 = %f, CDELT2 = %f: T_REC = %s, recnum = %lld \n", xscale, yscale, trecstr, inrec->recnum);
|
1315 tplarson 1.1 drms_free_array(inarr);
|
1316 tplarson 1.4 error++;
|
1317 tplarson 1.8 goto continue_outer_loop;
|
1318 tplarson 1.1 }
|
1319 tplarson 1.13 imagescale=xscale;
|
1320 tplarson 1.1 xscale=1.0;
1321 yscale=1.0;
1322
|
1323 tplarson 1.6 /*
|
1324 tplarson 1.1 imagescale = drms_getkey_double(inrec, "IM_SCALE", &status);
1325 PARAMETER_ERROR("IM_SCALE")
|
1326 tplarson 1.6 */
|
1327 tplarson 1.13
|
1328 tplarson 1.1 if (paramsign != 0)
1329 {
1330 vsign = paramsign;
1331 }
1332 else
1333 {
1334 vsign = drms_getkey_double(inrec, "DATASIGN", &status);
|
1335 tplarson 1.8 ParameterDef(status, "DATASIGN", 1.0, &vsign, inrec->recnum, 1);
|
1336 tplarson 1.1 }
1337
1338 if (velocity_correction)
1339 {
1340 obs_vr = drms_getkey_double(inrec, "OBS_VR", &status);
|
1341 tplarson 1.8 ParameterDef(status, "OBS_VR", 0.0, &obs_vr, inrec->recnum, 1);
|
1342 tplarson 1.1
1343 obs_vw = drms_getkey_double(inrec, "OBS_VW", &status);
|
1344 tplarson 1.8 ParameterDef(status, "OBS_VW", 0.0, &obs_vw, inrec->recnum, 1);
|
1345 tplarson 1.1
1346 obs_vn = drms_getkey_double(inrec, "OBS_VN", &status);
|
1347 tplarson 1.8 ParameterDef(status, "OBS_VN", 0.0, &obs_vn, inrec->recnum, 1);
|
1348 tplarson 1.1 }
1349
1350 // MDI keyword was OBS_DIST, in AU
|
1351 tplarson 1.21 obsdist = drms_getkey_double(inrec, "DSUN_OBS", &status) / AU;
1352 // 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
1353 ParameterDef(status, "OBS_DIST", 1.0, &obsdist, inrec->recnum, 1);
1354 if (readrsunref)
1355 {
1356 rsunref = drms_getkey_double(inrec, "RSUN_REF", &status);
1357 ParameterDef(status, "RSUN_REF", RTRUE, &rsunref, inrec->recnum, 1);
1358 rtrue=rsunref/AU;
1359 }
1360 S = rtrue / obsdist; // radians - approx. arcsin(rtrue/obsdist), but don't undo this approximation, because it is assumed in obs2helio
|
1361 tplarson 1.1
|
1362 tplarson 1.21 rsunobs = drms_getkey_double(inrec, "RSUN_OBS", &status);
|
1363 tplarson 1.9 if (status == DRMS_SUCCESS)
|
1364 tplarson 1.21 rsun = rsunobs/imagescale; //this calculation of rsun assumes approximation of imagescale mentioned in comment above
|
1365 tplarson 1.9 else
1366 {
1367 rsun = drms_getkey_double(inrec, "R_SUN", &status);
1368 if (status != DRMS_SUCCESS)
1369 rsun = ARCSECSINRAD * S / sqrt(1.0 - S * S) / imagescale;
1370 }
1371
|
1372 tplarson 1.1 if (longitude_shift == 1)
1373 {
1374 tmearth = tobs+TAU_A*(1.0-obsdist);
1375 longshift = (obsl0-refl0)+longrate*(tmearth-tref); // degrees
1376 while (longshift > 180.0) longshift-=360.0;
1377 while (longshift < -180.0) longshift+=360.0;
1378 }
1379 else if (longitude_shift == 2) // Shift center to nearest Carrington Degree
1380 {
1381 longshift = obsl0 - (int)(obsl0);
1382 if (longshift > 0.5) longshift -= 1.0;
1383 }
1384 else if (longitude_shift == 3) // Shift center to nearest tenth of a degree
1385 {
1386 longshift = (obsl0 * 10 - (int)(obsl0 * 10)) / 10;
1387 if (longshift > 0.5) longshift -= 1.0;
1388 }
1389 else
1390 {
1391 longshift = 0.0;
1392 }
1393 tplarson 1.1
1394 mapl0 = obsl0 - longshift;
1395
1396 xpixels = inarr->axis[0];
1397 ypixels = inarr->axis[1];
1398 pixels = xpixels * ypixels;
1399
1400 // MDI keyword was X0
1401 x0 = drms_getkey_double(inrec, "CRPIX1", &status);
|
1402 tplarson 1.8 ParameterDef(status, "CRPIX1", xpixels / 2, &x0, inrec->recnum, 1);
|
1403 tplarson 1.1 x0 -= 1.0;
1404
1405 // MDI keyword was Y0
1406 y0 = drms_getkey_double(inrec, "CRPIX2", &status);
|
1407 tplarson 1.8 ParameterDef(status, "CRPIX2", ypixels / 2, &y0, inrec->recnum, 1);
|
1408 tplarson 1.1 y0 -= 1.0;
1409
1410 if (mag_offset)
1411 {
1412 float *dat = (float *)inarr->data;
1413 double bfitzero = drms_getkey_double(inrec, "BFITZERO", &status);
1414 PARAMETER_ERROR("BFITZERO")
1415 int i;
1416
1417 if (!isnan(bfitzero))
1418 {
1419 for (i = 0; i < pixels; ++i)
1420 {
1421 dat[i] -= (float)bfitzero;
1422 }
1423
1424 }
1425 }
1426
1427 if (checko)
1428 {
1429 tplarson 1.1 orientation = drms_getkey_string(inrec, "ORIENT", &status);
1430 PARAMETER_ERROR("ORIENT")
1431 CheckO(orientation, &vstat);
1432 if (vstat != V2HStatus_Success)
1433 {
|
1434 tplarson 1.8 fprintf(stderr,"ERROR: illegal ORIENT: T_REC = %s, recnum = %lld\n", trecstr, inrec->recnum);
|
1435 tplarson 1.1 drms_free_array(inarr);
1436 free(orientation);
|
1437 tplarson 1.4 error++;
|
1438 tplarson 1.8 goto continue_outer_loop;
|
1439 tplarson 1.1 }
1440 }
1441 else
1442 {
1443 orientation = strdup(orientationdef);
1444 }
1445
|
1446 tplarson 1.21 if (status = obs2helio((float *)inarr->data,
|
1447 tplarson 1.1 v2hptr,
1448 xpixels,
1449 ypixels,
1450 x0,
1451 y0,
1452 b0 * RADSINDEG,
1453 p * RADSINDEG,
1454 S,
1455 rsun,
1456 rmax,
1457 interpolation,
1458 mapcols,
1459 maprows,
1460 longmin_adjusted * RADSINDEG,
1461 longinterval * RADSINDEG,
1462 longshift * RADSINDEG,
1463 sinBdelta,
1464 smajor,
1465 sminor,
1466 sangle * RADSINDEG,
1467 xscale,
1468 tplarson 1.1 yscale,
1469 orientation,
1470 mag_correction,
1471 velocity_correction,
1472 obs_vr,
1473 obs_vw,
1474 obs_vn,
1475 vsign,
1476 NaN_beyond_rmax,
1477 carrStretch,
1478 &distP,
1479 diffrotA,
1480 diffrotB,
1481 diffrotC,
1482 NULL,
1483 0))
1484 {
|
1485 tplarson 1.14 fprintf(stderr, "ERROR: failure in obs2helio: status = %d, T_REC = %s, recnum = %lld\n", status, trecstr, inrec->recnum);
|
1486 tplarson 1.1 drms_free_array(inarr);
1487 free(orientation);
|
1488 tplarson 1.4 error++;
|
1489 tplarson 1.8 goto continue_outer_loop;
|
1490 tplarson 1.1 }
1491 drms_free_array(inarr);
1492 free(orientation);
1493
1494 if (status = apodize(v2hptr,
1495 b0 * RADSINDEG,
1496 mapcols,
1497 maprows,
1498 longmin_adjusted * RADSINDEG,
1499 longinterval * RADSINDEG,
1500 sinBdelta,
1501 apodization,
1502 apinner,
1503 apwidth,
1504 apel,
1505 apx,
1506 apy))
1507 {
|
1508 tplarson 1.14 fprintf(stderr, "ERROR: failure in apodize: status = %d, T_REC = %s, recnum = %lld\n", status, trecstr, inrec->recnum);
|
1509 tplarson 1.4 error++;
|
1510 tplarson 1.8 goto continue_outer_loop;
|
1511 tplarson 1.1 }
1512
1513 if (v2hflag)
1514 {
|
1515 tplarson 1.14 // outrec = drms_create_record(drms_env, v2hout, lifetime, &status);
1516 outrec=v2hrecset->records[iv2hrec++];
|
1517 tplarson 1.20 drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
|
1518 tplarson 1.1 DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
1519 DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
|
1520 tplarson 1.14 if (histlink != NULL)
|
1521 tplarson 1.1 drms_setlink_static(outrec, histlinkname, histrecnum);
|
1522 tplarson 1.14 if (srclink != NULL)
|
1523 tplarson 1.1 drms_setlink_static(outrec, srclinkname, inrec->recnum);
1524 if (segoutflag)
1525 segout = drms_segment_lookup(outrec, segnameout);
1526 else
1527 segout = drms_segment_lookupnum(outrec, 0);
1528 length[0]=mapcols;
1529 length[1]=maprows;
|
1530 tplarson 1.8 outarr = drms_array_create(usetype, 2, length, v2hptr, &status);
|
1531 tplarson 1.22 drms_setkey_int(outrec, "TOTVALS", maprows*mapcols);
|
1532 tplarson 1.9 set_statistics(segout, outarr, 1);
|
1533 tplarson 1.16 outarr->bzero=segout->bzero;
1534 outarr->bscale=segout->bscale;
|
1535 tplarson 1.8 status=drms_segment_write(segout, outarr, 0);
1536 free(outarr);
1537
1538 if (status != DRMS_SUCCESS)
1539 {
1540 fprintf(stderr, "ERROR: problem writing output segment: status = %d, T_REC = %s, input recnum = %lld, histrecnum = %lld\n", status, trecstr, inrec->recnum, histrecnum);
1541 return 0;
1542 }
1543
|
1544 tplarson 1.15 // drms_copykey(outrec, inrec, "T_REC");
|
1545 tplarson 1.1 drms_setkey_int(outrec, "QUALITY", quality);
1546 drms_setkey_int(outrec, "MAPMMAX", mapped_lmax);
1547 drms_setkey_int(outrec, "SINBDIVS", sinb_divisions);
1548 drms_setkey_float(outrec, "CRPIX1", mapcols/2.0 + 0.5);
1549 drms_setkey_float(outrec, "CRVAL1", mapl0);
1550 drms_setkey_float(outrec, "CROTA1", 0.0);
1551 drms_setkey_float(outrec, "CDELT1", longinterval);
1552 drms_setkey_float(outrec, "CRPIX2", maprows/2.0 + 0.5);
1553 drms_setkey_float(outrec, "CRVAL2", 0.0);
1554 drms_setkey_float(outrec, "CROTA2", 0.0);
1555 drms_setkey_float(outrec, "CDELT2", sinBdelta);
1556 drms_setkey_string(outrec, "CTYPE1", "CRLN_CEA");
1557 drms_setkey_string(outrec, "CTYPE2", "CRLT_CEA");
1558 drms_setkey_string(outrec, "CUNIT1", "deg");
1559 drms_setkey_string(outrec, "CUNIT2", "sinlat");
|
1560 tplarson 1.17
|
1561 tplarson 1.22 // set keywords for magnetic pipeline
1562 drms_setkey_float(outrec, "MAPRMAX", rmax);
1563 drms_setkey_float(outrec, "MAPBMAX", bmax);
1564 drms_setkey_float(outrec, "MAPLGMAX", longmax_adjusted);
1565 drms_setkey_float(outrec, "MAPLGMIN", longmin_adjusted);
1566 drms_setkey_int(outrec, "INTERPO", interpolation);
1567 drms_setkey_int(outrec, "LGSHIFT", longitude_shift);
1568 drms_setkey_int(outrec, "MCORLEV", mag_correction);
1569 drms_setkey_int(outrec, "MOFFSET", mag_offset);
1570 drms_setkey_int(outrec, "CARSTRCH", carrStretch);
1571 drms_setkey_float(outrec, "DIFROT_A", diffrotA);
1572 drms_setkey_float(outrec, "DIFROT_B", diffrotB);
1573 drms_setkey_float(outrec, "DIFROT_C", diffrotC);
1574
|
1575 tplarson 1.1 dsignout=vsign*drms_getkey_double(inrec, "DATASIGN", &status);
1576 if (status != DRMS_SUCCESS)
1577 dsignout=vsign;
1578 dsignout/=fabs(dsignout);
1579 drms_setkey_int(outrec, "DATASIGN", (int)dsignout);
1580
1581 tnow = (double)time(NULL);
1582 tnow += UNIX_epoch;
1583 drms_setkey_time(outrec, "DATE", tnow);
1584
|
1585 tplarson 1.14 // drms_close_record(outrec, DRMS_INSERT_RECORD);
|
1586 tplarson 1.1 }
1587
1588
1589 if (verbflag > 1)
1590 {
1591 wt=getwalltime();
1592 ct=getcputime(&ut, &st);
1593 fprintf(stdout,
1594 " remap done, %.2f ms wall time, %.2f ms cpu time\n",
1595 wt-wt2,
1596 ct-ct2);
1597 }
1598
|
1599 tplarson 1.9 if (!h2mflag && !tsflag)
1600 goto skip;
1601
|
1602 tplarson 1.1 inp=v2hptr;
1603 outp=h2mptr;
1604
1605 mean = 0.0;
1606 if (subtract_mean) /* get mean of entire remapped image */
1607 {
1608 nmean = 0;
1609 ip = inp;
1610 for (row=0; row<maprows; row++)
1611 {
1612 for (col=0; col<mapcols; col++)
1613 {
1614 if (!isnan(*ip))
1615 {
1616 nmean++;
1617 mean += *ip;
1618 }
1619 ip++;
1620 }
1621 }
1622 mean /= (nmean ? nmean : 1);
1623 tplarson 1.1 }
1624
1625 for (row=0; row<maprows; row++)
1626 {
1627
1628 if (cent_long == 0)
1629 {
1630
1631 /* Old code with 0 at beginning of array */
1632
1633 nok = 0;
1634 bp = buf;
1635 wp = weight;
1636 ip = inp + mapcols * row;
1637 for (col=0; col<mapcols; col++)
1638 {
1639 if (!isnan(*ip))
1640 {
1641 *bp++ = *wp * (*ip - mean);
1642 nok += 1;
1643 }
1644 tplarson 1.1 else
1645 *bp++ = 0.0;
1646 ip++;
1647 wp++;
1648 }
1649
1650 /* zero fill */
1651 for (i=col; i<nfft; i++)
1652 *bp++ = 0.0;
1653 }
1654 else
1655 {
1656
1657 /* New code with 0 at center meridian */
1658 /* Assumes that input array is symmetric around meridian */
1659
1660 nok = 0;
1661
1662 /* First copy right side of meridian */
1663
1664 bp = buf;
1665 tplarson 1.1 ip = inp + mapcols * row+mapcols2;
1666 wp = weight + mapcols2;
1667 if (subtract_mean)
1668 {
1669 for (col=0; col<=mapcols2; col++)
1670 {
1671 if (!isnan(*ip))
1672 {
1673 *bp++ = *wp * (*ip - mean);
1674 nok += 1;
1675 }
1676 else
1677 *bp++ = 0.0;
1678 ip++;
1679 wp++;
1680 }
1681 }
1682 else
1683 {
1684 for (col=0; col<=mapcols2; col++)
1685 {
1686 tplarson 1.1 if (!isnan(*ip))
1687 {
1688 *bp++ = *wp * *ip;
1689 nok += 1;
1690 }
1691 else
1692 *bp++ = 0.0;
1693 ip++;
1694 wp++;
1695 }
1696 }
1697
1698 /* Then do left side of meridian */
1699 bp = buf+lfft-mapcols2;
1700 ip = inp + mapcols * row;
1701 wp = weight;
1702 if (subtract_mean)
1703 {
1704 for (col=0; col<mapcols2; col++)
1705 {
1706 if (!isnan(*ip))
1707 tplarson 1.1 {
1708 *bp++ = *wp * (*ip - mean);
1709 nok += 1;
1710 }
1711 else
1712 *bp++ = 0.0;
1713 ip++;
1714 wp++;
1715 }
1716 }
1717 else
1718 {
1719 for (col=0; col<mapcols2; col++)
1720 {
1721 if (!isnan(*ip))
1722 {
1723 *bp++ = *wp * *ip;
1724 nok += 1;
1725 }
1726 else
1727 *bp++ = 0.0;
1728 tplarson 1.1 ip++;
1729 wp++;
1730 }
1731 }
1732
1733 /* Finally zero fill */
1734 bp = buf+mapcols2+1;
1735 for (i=0; i<lfft-mapcols; i++)
1736 *bp++ = 0.0;
1737
1738 } /* End of copying if statement */
1739
1740 if ((zero_miss == 0) && (nok != mapcols))
1741 {
1742
1743 /* Stuff with missing */
1744 for (i=0; i<nout; i++)
1745 {
1746 op = outp + i * maprows + row;
1747 *op = DRMS_MISSING_FLOAT;
1748 }
1749 tplarson 1.1 }
1750 else
1751 {
1752 if (normalize)
1753 norm = 1./sqrt(2*(double)nfft/nok)/lfft;
1754 else
1755 norm = 1./lfft;
|
1756 tplarson 1.8
|
1757 tplarson 1.1 /* Fourier transform */
1758 fftwf_execute(fftwp);
1759
1760 /* transpose, normalize, this is where most time is spent */
1761 /* First do real part */
1762 for (i=0; i<nout/2; i++)
1763 {
1764 op = outp + 2*i * maprows + row;
1765 *op = wbuf[i]*norm;
1766 }
1767 /* Do imaginary part */
1768 /* Imaginary part of m=0 is 0*/
1769 *(outp+row+maprows)=0.0;
1770 /* Use normx to get the complex conjugate */
1771 normx=-norm;
1772 for (i=1; i<nout/2; i++)
1773 {
1774 op = outp + 2*i * maprows + row + maprows;
1775 *op = wbuf[lfft-i]*normx;
1776 }
1777
1778 tplarson 1.1 }
1779 } /* Next row */
1780
1781 if (h2mflag)
1782 {
|
1783 tplarson 1.14 // outrec = drms_create_record(drms_env, h2mout, lifetime, &status);
1784 outrec=h2mrecset->records[ih2mrec++];
|
1785 tplarson 1.20 drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
|
1786 tplarson 1.1 DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
1787 DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
|
1788 tplarson 1.14 if (histlink != NULL)
|
1789 tplarson 1.1 drms_setlink_static(outrec, histlinkname, histrecnum);
|
1790 tplarson 1.14 if (srclink != NULL)
|
1791 tplarson 1.1 drms_setlink_static(outrec, srclinkname, inrec->recnum);
1792 if (segoutflag)
1793 segout = drms_segment_lookup(outrec, segnameout);
1794 else
1795 segout = drms_segment_lookupnum(outrec, 0);
1796 length[0]=maprows;
1797 length[1]=nout;
|
1798 tplarson 1.8 outarr = drms_array_create(usetype, 2, length, h2mptr, &status);
|
1799 tplarson 1.22 drms_setkey_int(outrec, "TOTVALS", maprows*nout);
|
1800 tplarson 1.9 set_statistics(segout, outarr, 1);
|
1801 tplarson 1.16 outarr->bzero=segout->bzero;
1802 outarr->bscale=segout->bscale;
|
1803 tplarson 1.8 status=drms_segment_write(segout, outarr, 0);
1804 free(outarr);
1805
1806 if (status != DRMS_SUCCESS)
1807 {
1808 fprintf(stderr, "ERROR: problem writing output segment: status = %d, T_REC = %s, input recnum = %lld, histrecnum = %lld\n", status, trecstr, inrec->recnum, histrecnum);
1809 return 0;
1810 }
1811
|
1812 tplarson 1.15 // drms_copykey(outrec, inrec, "T_REC");
|
1813 tplarson 1.1 drms_setkey_int(outrec, "QUALITY", quality);
1814 drms_setkey_int(outrec, "MAPMMAX", mapped_lmax);
1815 drms_setkey_int(outrec, "SINBDIVS", sinb_divisions);
1816 drms_setkey_int(outrec, "LMAX", lmax);
1817 drms_setkey_double(outrec, "CRPIX1", maprows/2.0 + 0.5);
1818 drms_setkey_double(outrec, "CRVAL1", 0.0);
1819 drms_setkey_double(outrec, "CROTA1", 0.0);
1820 drms_setkey_double(outrec, "CDELT1", sinBdelta);
1821 drms_setkey_double(outrec, "CRPIX2", 1.0);
1822 drms_setkey_double(outrec, "CRVAL2", 0.0);
1823 drms_setkey_double(outrec, "CROTA2", 0.0);
1824 drms_setkey_double(outrec, "CDELT2", 1.0);
1825 drms_setkey_string(outrec, "CTYPE1", "CRLT_CEA");
1826 drms_setkey_string(outrec, "CTYPE2", "CRLN_FFT");
1827 drms_setkey_string(outrec, "CUNIT1", "rad");
1828 drms_setkey_string(outrec, "CUNIT2", "m");
|
1829 tplarson 1.17
|
1830 tplarson 1.1 tnow = (double)time(NULL);
1831 tnow += UNIX_epoch;
1832 drms_setkey_time(outrec, "DATE", tnow);
|
1833 tplarson 1.14 // drms_close_record(outrec, DRMS_INSERT_RECORD);
|
1834 tplarson 1.1 }
1835
1836 if (verbflag > 1)
1837 {
1838 wt2=getwalltime();
1839 ct2=getcputime(&ut2, &st2);
1840 fprintf(stdout,
1841 " fft and transpose done, %.2f ms wall time, %.2f ms cpu time\n",
1842 wt2-wt,
1843 ct2-ct);
1844 }
1845
|
1846 tplarson 1.9 if (!tsflag)
1847 goto skip;
1848
|
1849 tplarson 1.1 inptr=h2mptr;
1850 imageoffset = imagesize * irec;
1851 for (mx = 0; mx < 2*msize; mx++) /* for each m, re and im */
1852 {
1853 moffset = mx * maprows;
1854 mptr = inptr + moffset;
1855 for (latx = 0; latx < nlat; latx++)
1856 {
1857 poslatx = nlat + latx; neglatx = nlat - 1 - latx;
1858 evenpart[latx] = mptr[poslatx] + mptr[neglatx];
1859 oddpart[latx] = mptr[poslatx] - mptr[neglatx];
1860 }
1861 workptr = folded + imageoffset + moffset;
1862 scopy_ (&nlat, evenpart, &increment1, workptr, &increment1);
1863 workptr += nlat;
1864 scopy_ (&nlat, oddpart, &increment1, workptr, &increment1);
1865 }
1866
1867 skip:
|
1868 tplarson 1.8 inrec = drms_recordset_fetchnext(drms_env, inrecset, &fetchstat, &chunkstat, NULL);
|
1869 tplarson 1.1
|
1870 tplarson 1.4 if (inrec != NULL)
|
1871 tplarson 1.1 {
1872 trec = drms_getkey_time(inrec, "T_REC", &status);
1873 PARAMETER_ERROR("T_REC")
|
1874 tplarson 1.8 trecind=(trec-tstart+cadence/2)/cadence;
|
1875 tplarson 1.1 sprint_time(trecstr, trec, "TAI", 0);
1876 }
|
1877 tplarson 1.4
|
1878 tplarson 1.1
1879 } /* end loop on each input image for this timechunk */
1880
1881 if (verbflag)
1882 printf(" number of bad images = %d\n", bc);
|
1883 tplarson 1.4
|
1884 tplarson 1.9 if (!tsflag)
1885 goto continue_outer_loop;
1886
|
1887 tplarson 1.8 //needed if recordset does not extend to end of a timechunk
1888 while (irec < nrecs)
1889 {
1890 bad[bc++]=irec;
1891 irec++;
1892 }
1893
|
1894 tplarson 1.4 if (bc == nrecs)
|
1895 tplarson 1.1 {
|
1896 tplarson 1.4 nodata=1;
1897 }
1898 else
1899 {
1900 while (bc > 0)
1901 {
1902 imageoffset=imagesize*bad[--bc];
1903 for (i=0;i<imagesize;i++) folded[imageoffset+i]=0.0;
1904 }
|
1905 tplarson 1.1 }
1906
1907 if (verbflag)
1908 {
1909 wt2=getwalltime();
1910 ct2=getcputime(&ut2, &st2);
1911 fprintf(stdout,
1912 " images processed, %.2f ms wall time, %.2f ms cpu time\n",
1913 wt2-wt1,
1914 ct2-ct1);
1915 }
1916
|
1917 tplarson 1.8 // skip to here if no input records in a given timechunk
1918 skip_norecs:
|
1919 tplarson 1.1
1920 /* we now have folded data for a chunk of sn's */
1921 /* now do Jesper's tricks */
1922
1923 /* ldone is the last l for which plm's have been set up */
1924 ldone=-1;
1925
1926 /* loop on each chunk of l's */
1927 for (lchunk = lchunkfirst; lchunk <= lchunklast; lchunk++)
1928 {
1929 lfirst = lchunk * lchunksize;
1930 llast = lfirst + lchunksize - 1;
1931 lfirst = maxval(lfirst,lmin);
1932 llast = minval(llast,lmax);
1933 /* get the first and last indexes into the l-m array */
1934 ifirst = lfirst*(lfirst+1)/2;
1935 ilast = llast*(llast+1)/2+llast;
1936
1937 if (verbflag > 1)
1938 {
1939 wt3=getwalltime();
1940 tplarson 1.1 ct3=getcputime(&ut3, &st3);
1941 printf(" processing lchunk %d, lmin = %d, lmax = %d\n", lchunk, lfirst, llast);
1942 }
1943
|
1944 tplarson 1.14 // outrec = drms_create_record(drms_env, outseries, lifetime, &status);
1945 outrec=outrecset->records[irecout++];
1946 /*
|
1947 tplarson 1.1 if (status != DRMS_SUCCESS || outrec == NULL)
1948 {
|
1949 tplarson 1.14 fprintf(stderr,"ERROR: unable to open record in output dataseries %s, status = %d, histrecnum = %lld\n", outseries, status, histrecnum);
|
1950 tplarson 1.1 return 0;
1951 }
|
1952 tplarson 1.14 */
1953 if (histlink != NULL)
|
1954 tplarson 1.1 drms_setlink_static(outrec, histlinkname, histrecnum);
1955
|
1956 tplarson 1.4 if (nodata)
|
1957 tplarson 1.8 goto skip_nodata;
|
1958 tplarson 1.4
|
1959 tplarson 1.1 /* now the size of the output array is known */
1960 length[0]=2*nrecs; /* accomodate re & im parts for each sn */
1961 length[1]=(ilast-ifirst+1); /* for each l & m, lfirst <= l <= llast */
1962
1963 outarr = drms_array_create(usetype, 2, length, NULL, &status);
1964
1965 if (segoutflag)
1966 segout = drms_segment_lookup(outrec, segnameout);
1967 else
1968 segout = drms_segment_lookupnum(outrec, 0);
1969
|
1970 tplarson 1.8 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
|
1971 tplarson 1.1 {
|
1972 tplarson 1.14 fprintf(stderr,"ERROR: problem with output segment or data array: lfirst = %d, llast = %d, length = [%d, %d], status = %d, iset = %d, T_START= %s, histrecnum = %lld",
1973 lfirst, llast, length[0], length[1], status, iset, tstartstr, histrecnum);
|
1974 tplarson 1.1 return 0;
1975 }
1976
1977 outptr = (float *)(outarr->data);
1978
1979 /* loop on each m */
1980 for (m = 0; m <= llast; m++)
1981 {
1982
1983 modd = is_odd(m);
1984 meven = !modd;
1985 lstart = maxval(lfirst,m); /* no l can be smaller than this m */
1986
1987 /* set up masks (plm's) for this m and chunk in l */
1988
1989 if ((lstart-1) == ldone)
1990 {
1991 /* get saved plms if any */
1992 if ((lstart - 2) >= m)
1993 for (latx = 0; latx < nlat; latx++)
1994 plm[(lstart-2)*nlat + latx] = saveplm[m*nlat + latx];
1995 tplarson 1.1 if ((lstart - 1) >= m)
1996 for (latx = 0; latx < nlat; latx++)
1997 plm[(lstart-1)*nlat + latx] = saveplm[msize*nlat + m*nlat + latx];
1998
1999 /* then set up the current chunk */
|
2000 tplarson 1.17 // setplm_ (&lstart, &llast, &m, &nlat, indx, latgrid, &nlat, plm);
2001 setplm2(lstart, llast, m, nlat, indx, latgrid, nlat, plm, NULL);
|
2002 tplarson 1.1 }
2003 else
2004 {
2005 /* This fixes the lmin != 0 problem */
|
2006 tplarson 1.17 // setplm_ (&m, &llast, &m, &nlat, indx, latgrid, &nlat, plm);
2007 setplm2(m, llast, m, nlat, indx, latgrid, nlat, plm, NULL);
|
2008 tplarson 1.1 }
2009
2010 /* save plm's for next chunk */
2011 if ((llast-1) >= m)
2012 for (latx = 0; latx < nlat; latx++)
2013 saveplm[m*nlat + latx] = plm[(llast - 1)*nlat + latx];
2014 for (latx = 0; latx < nlat; latx++)
2015 saveplm[msize*nlat + m*nlat + latx] = plm[llast*nlat + latx];
2016 ldone=llast;
2017
2018 /* copy plm's into masks */
2019 /* note that this converts from double to single precision */
2020 /* the test prevents underflows which gobble CPU time */
2021 /* Hmmm... looks like if statement is not needed. Weird...
2022 for (l = lstart; l <= llast; l++) {
2023 moffset = l * nlat;
2024 for (latx = 0; latx < nlat; latx++) {
2025 plmptr = plm + moffset + latx;
2026 if (is_very_small (*plmptr))
2027 masks [(l-lstart)*nlat + latx] = 0.0;
2028 else
2029 tplarson 1.1 masks [(l-lstart)*nlat + latx] = *plmptr;
2030 }
2031 plmptr = plm + moffset;
2032 maskptr = masks+(l-lstart)*nlat;
2033 dscopy_(&nlat,plmptr,maskptr);
2034 }
2035 */
2036 plmptr=plm+nlat*lstart;
2037 maskptr=masks;
2038 latx=nlat*(llast-lstart+1);
2039 // dscopy_(&latx,plmptr,maskptr);
2040 int ilatx;
2041 for (ilatx=0;ilatx<latx;ilatx++)
2042 maskptr[ilatx]=plmptr[ilatx];
2043
2044 /* for each sn in snchunk */
2045 // for (sn = infsn; sn <= inlsn; sn++) {
2046 for (snx=0; snx<nrecs; snx++)
2047 {
2048 /* select folded data for real/imag l's and this m
2049 into temporay arrays for matrix multiply */
2050 tplarson 1.1 // snx = sn - infsn;
2051 /* TO DO - pull offset calculations out of loop */
2052 /* New code with odd leading dimension */
2053 scopy_ (&nlat,
2054 folded + nlat*(4*m+modd) + snx*imagesize,
2055 &increment1,
2056 real4evenl + snx*nlatx,
2057 &increment1);
2058 scopy_ (&nlat,
2059 folded + nlat*(4*m+meven) + snx*imagesize,
2060 &increment1,
2061 real4oddl + snx*nlatx,
2062 &increment1);
2063 scopy_ (&nlat,
2064 folded + nlat*(4*m+2+modd) + snx*imagesize,
2065 &increment1,
2066 imag4evenl + snx*nlatx,
2067 &increment1);
2068 scopy_ (&nlat,
2069 folded + nlat*(4*m+2+meven) + snx*imagesize,
2070 &increment1,
2071 tplarson 1.1 imag4oddl + snx*nlatx,
2072 &increment1);
2073 } /* end loop through snchunk */
2074
2075
2076 /* do even l's */
2077 lfirsteven = is_even(lstart) ? lstart : lstart+1;
2078 llasteven = is_even(llast) ? llast : llast-1;
2079 nevenl = (llasteven-lfirsteven)/2 + 1; /* number of even l's */
2080 /* do real part */
2081 /* All parts used to have alpha=&one, now have alpha=&cnorm */
2082 sgemm_ (transpose, /* form of op(A) */
2083 normal, /* form of op(B) */
2084 &nsn, /* number of sn's */
2085 &nevenl, /* number of even l's for this m */
2086 &nlat, /* number of latitudes */
2087 &cnorm, /* scalar multiplier of op(A) */
2088 real4evenl, /* matrix A */
2089 &nlatx, /* use every nlat-long row of A */
2090 masks + nlat*(lfirsteven-lstart), /* matrix B */
2091 &maprows, /* 2*nlat, use every other row (nlat long) of B */
2092 tplarson 1.1 &zero, /* scalar multiplier of C */
2093 outx + nsn*2*(lfirsteven-lstart), /* matrix C (output) */
2094 &fournsn, /* use every fourth nsn-long row of C */
2095 1, /* length of transpose character string */
2096 1); /* length of normal character string */
2097 /* do imag part */
2098 sgemm_ (transpose, /* form of op(A) */
2099 normal, /* form of op(B) */
2100 &nsn, /* number of sn's */
2101 &nevenl, /* number of even l's for this m */
2102 &nlat, /* number of latitudes */
2103 &cnorm, /* scalar multiplier of op(A) */
2104 imag4evenl, /* matrix A */
2105 &nlatx, /* use every nlat-long row of A */
2106 masks + nlat*(lfirsteven-lstart), /* matrix B */
2107 &maprows, /* 2*nlat, use every other nlat-long row of B */
2108 &zero, /* scalar multiplier of C */
2109 outx + nsn*(2*(lfirsteven-lstart)+1), /* matrix C (output) */
2110 &fournsn, /* use every fourth nsn-long row of C */
2111 1, /* length of transpose character string */
2112 1); /* length of normal character string */
2113 tplarson 1.1
2114 /* do odd l's */
2115 lfirstodd = is_odd(lstart) ? lstart : lstart+1;
2116 llastodd = is_odd(llast) ? llast : llast-1;
2117 noddl = (llastodd-lfirstodd)/2 + 1; /* number of odd l's */
2118 /* do real part */
2119 sgemm_ (transpose, /* form of op(A) */
2120 normal, /* form of op(B) */
2121 &nsn, /* number of sn's */
2122 &noddl, /* number of odd l's for this m */
2123 &nlat, /* number of latitudes */
2124 &cnorm , /* scalar multiplier of op(A) */
2125 real4oddl, /* matrix A */
2126 &nlatx, /* use every nlat-long row of A */
2127 masks + nlat*(lfirstodd-lstart), /* matrix B */
2128 &maprows, /* 2*nlat, use every other nlat-long row of B */
2129 &zero, /* scalar multiplier of C */
2130 outx + nsn*2*(lfirstodd-lstart), /* matrix C (output) */
2131 &fournsn, /* use every fourth nsn-long row of C */
2132 1, /* length of transpose character string */
2133 1); /* length of normal character string */
2134 tplarson 1.1 /* do imag part */
2135 sgemm_ (transpose, /* form of op(A) */
2136 normal, /* form of op(B) */
2137 &nsn, /* number of sn's */
2138 &noddl, /* number of odd l's for this m */
2139 &nlat, /* number of latitudes */
2140 &cnorm, /* scalar multiplier of op(A) */
2141 imag4oddl, /* matrix A */
2142 &nlatx, /* use every nlat-long row of A */
2143 masks + nlat*(lfirstodd-lstart), /* matrix B */
2144 &maprows, /* 2*nlat, use every other nlat-long row of B */
2145 &zero, /* scalar multiplier of C */
2146 outx + nsn*(2*(lfirstodd-lstart)+1), /* matrix C (output) */
2147 &fournsn, /* use every fourth nsn-long row of C */
2148 1, /* length of transpose character string */
2149 1); /* length of normal character string */
2150
2151 /* copy outx into out sds */
2152 /* alternate real and imaginary values in out - as in pipeLNU */
2153 for (l = lstart; l <= llast; l++)
2154 {
2155 tplarson 1.1 fromoffset = 2*nsn*(l-lstart);
2156 tooffset = 2*nsn*(l*(l+1)/2 + m -ifirst);
2157 scopy_ (&nsn,
2158 outx+fromoffset,
2159 &increment1,
2160 outptr+tooffset,
2161 &increment2);
2162 scopy_ (&nsn,
2163 outx+fromoffset+nsn,
2164 &increment1,
2165 outptr+tooffset+1,
2166 &increment2);
2167 } /* end loop through l's for this m */
2168
2169
2170 } /* end loop on m */
2171
|
2172 tplarson 1.17
|
2173 tplarson 1.8 outarr->bzero=segout->bzero;
2174 outarr->bscale=segout->bscale;
2175 status=drms_segment_write(segout, outarr, 0);
|
2176 tplarson 1.4 drms_free_array(outarr);
2177 nsegments++;
2178
|
2179 tplarson 1.8 if (status != DRMS_SUCCESS)
2180 {
2181 fprintf(stderr, "ERROR: problem writing output segment: status = %d, T_START = %s, LMIN = %d, LMAX = %d, histrecnum = %lld\n", status, tstartstr, lfirst, llast, histrecnum);
2182 return 0;
2183 }
2184
2185 skip_nodata:
|
2186 tplarson 1.4
|
2187 tplarson 1.1 drms_setkey_int(outrec, "LMIN", lfirst);
2188 drms_setkey_int(outrec, "LMAX", llast);
2189 drms_setkey_time(outrec, "T_START", tstart);
2190 drms_setkey_time(outrec, "T_STOP", tstart+chunksecs);
2191 drms_setkey_time(outrec, "T_OBS", tstart+chunksecs/2);
2192 drms_setkey_time(outrec, "DATE__OBS", tstart);
|
2193 tplarson 1.8 drms_setkey_string(outrec, "TAG", tag);
|
2194 tplarson 1.23 if (verflag)
2195 drms_setkey_string(outrec, "VERSION", version);
|
2196 tplarson 1.20
2197 for (i=0;i<16;i++)
2198 setbits(calversout,4*i+3,4,nybblearrout[i]);
|
2199 tplarson 1.16 drms_setkey_longlong(outrec, "CALVER64", calversout);
|
2200 tplarson 1.1
|
2201 tplarson 1.4 if (nodata)
2202 drms_setkey_int(outrec, "QUALITY", QUAL_NODATA);
|
2203 tplarson 1.20 else if (mixflag)
2204 drms_setkey_int(outrec, "QUALITY", QUAL_MIXEDCALVER);
|
2205 tplarson 1.4 else
2206 drms_setkey_int(outrec, "QUALITY", 0);
|
2207 tplarson 1.1
2208 // these could be constant, but set them just in case
|
2209 tplarson 1.6 drms_setkey_int(outrec, "MAPMMAX", mapped_lmax);
2210 drms_setkey_int(outrec, "SINBDIVS", sinb_divisions);
|
2211 tplarson 1.1 drms_setkey_float(outrec, "T_STEP", cadence);
|
2212 tplarson 1.8
|
2213 tplarson 1.6 ndt=chunksecs/cadence;
2214 drms_setkey_int(outrec, "NDT", ndt);
2215
|
2216 tplarson 1.1 tnow = (double)time(NULL);
2217 tnow += UNIX_epoch;
2218 drms_setkey_time(outrec, "DATE", tnow);
2219
|
2220 tplarson 1.14 // drms_close_record(outrec, DRMS_INSERT_RECORD);
|
2221 tplarson 1.4 nrecords++;
|
2222 tplarson 1.1
2223 if (verbflag > 1)
2224 {
2225 wt=getwalltime();
2226 ct=getcputime(&ut, &st);
2227 fprintf(stdout,
2228 " %.2f ms wall time, %.2f ms cpu time\n",
2229 wt-wt3,
2230 ct-ct3);
2231 }
2232
2233
2234 } /* end loop on each chunk of l's */
2235
2236 if (verbflag)
2237 {
2238 wt1=getwalltime();
2239 ct1=getcputime(&ut1, &st1);
2240 fprintf(stdout, "SHT of timechunk %d complete: %.2f ms wall time, %.2f ms cpu time\n", iset,
2241 wt1-wt2, ct1-ct2);
2242 }
2243 tplarson 1.1
|
2244 tplarson 1.8 continue_outer_loop:
|
2245 tplarson 1.1 tstart+=chunksecs;
2246 } /* end loop on each time chunk */
2247
2248
2249 if (chunkstat != kRecChunking_LastInRS && chunkstat != kRecChunking_NoMoreRecs)
2250 fprintf(stderr, "WARNING: input records remain after last output record: chunkstat = %d\n", (int)chunkstat);
2251
|
2252 tplarson 1.8 drms_close_records(inrecset, DRMS_FREE_RECORD);
|
2253 tplarson 1.14 if (tsflag)
2254 drms_close_records(outrecset, DRMS_INSERT_RECORD);
2255 if (v2hflag)
2256 drms_close_records(v2hrecset, DRMS_INSERT_RECORD);
2257 if (h2mflag)
2258 drms_close_records(h2mrecset, DRMS_INSERT_RECORD);
2259
|
2260 tplarson 1.1
2261 wt=getwalltime();
2262 ct=getcputime(&ut, &st);
|
2263 tplarson 1.13 if (verbflag && tsflag)
|
2264 tplarson 1.1 {
|
2265 tplarson 1.4 printf("number of records created = %d\n", nrecords);
2266 printf("number of segments created = %d\n", nsegments);
|
2267 tplarson 1.13 }
2268 if (verbflag)
2269 {
|
2270 tplarson 1.1 fprintf(stdout, "total time spent: %.2f ms wall time, %.2f ms cpu time\n",
2271 wt-wt0, ct-ct0);
2272 }
2273
|
2274 tplarson 1.8 if (!error)
2275 printf("module %s successful completion\n", cmdparams.argv[0]);
2276 else
2277 printf("module %s failed to produce %d timechunks: histrecnum = %lld\n", cmdparams.argv[0], error, histrecnum);
2278
|
2279 tplarson 1.1 return 0;
2280 }
|
Return to
jv2ts.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps