| 1 |
/***************************************************************** |
| 2 |
* rdvinv.c ~rick/src/rdvinv |
| 3 |
* |
| 4 |
* Responsible: Charles Baldner, Rick Bogart, Sarbani Basu |
| 5 |
* RBogart@spd.aas.org |
| 6 |
* |
| 7 |
* Velocity inversion code for ring diagram power spectra fits |
| 8 |
* |
| 9 |
* Parameters: (type default description) |
| 10 |
* in string - Input data set descriptor or filename |
| 11 |
* out string fort.10.hmixy Output filename |
| 12 |
* seg string fit.out Input data segment name |
| 13 |
* cr int NS Carrington rotation common to dataset |
| 14 |
* clon float NS CM Longitude common to dataset |
| 15 |
* lon float NS Carrington longitude common to dataset |
| 16 |
* lat float NS Heliographic latitude common to dataset |
| 17 |
* amu double 0.005 Error trade-off parameter |
| 18 |
* ob double 1.0 Lower frequency limit (mHz) |
| 19 |
* oe double 5.2 Upper frequency limit (mHz) |
| 20 |
* rb double 0.97 Lower radius limit |
| 21 |
* re double 1.00 Upper radius limit |
| 22 |
* num int 40 Number of target inversion points |
| 23 |
* kernel string - Data record or pathname of file |
| 24 |
* containing mode kernels for inversions |
| 25 |
* ave string Not Specified Output file for averaging kernels |
| 26 |
* (if not set, kernels are not written) |
| 27 |
* coef string Not Specified Output file for ?? coefficients |
| 28 |
* (if not set, coefficients are not written) |
| 29 |
* |
| 30 |
* Flags |
| 31 |
* -v run verbose |
| 32 |
* |
| 33 |
* Notes: |
| 34 |
* Output is specified by argment 'out'. If 'out' is a drms data series |
| 35 |
* that can be opened, the inversions will be saved to the appropriate |
| 36 |
* drms records; otherwise, 'out' will be used as the base filename for |
| 37 |
* the inversion outputs. |
| 38 |
* |
| 39 |
* Bugs: |
| 40 |
* If the input data is not in a drms series, an output series can't be |
| 41 |
* created - rdvinv will exit with an error message. |
| 42 |
* Argument 'out', when a file name, cannot contain any path information. |
| 43 |
* When writing output to a file, rather than a drms record set, only one |
| 44 |
* set of files can be specified; if the input contains more than one |
| 45 |
* data record, each inversion will overwrite the last. |
| 46 |
* There is no verification that the requested output segments have the |
| 47 |
* correct protocol; probably doesn't matter anyway |
| 48 |
* The constant value of amu is written needlessly for each target location |
| 49 |
* for backward compatability |
| 50 |
* The reading of the kernel file is still in the FORTRAN subroutine |
| 51 |
* |
| 52 |
* Revision history is at end of file |
| 53 |
******************************************************************/ |
| 54 |
|
| 55 |
|
| 56 |
#include <stdio.h> |
| 57 |
#include <string.h> |
| 58 |
#include <math.h> |
| 59 |
#include <jsoc_main.h> |
| 60 |
#include "keystuff.c" |
| 61 |
#include "old_rdutil.c" |
| 62 |
#include "ola_xy.c" |
| 63 |
|
| 64 |
/* prototypes */ |
| 65 |
|
| 66 |
extern void ola_(double *, int *, double *, double *, double *, |
| 67 |
double *, double *, int *, char *, double *, double *, double *, |
| 68 |
double *, double *, double *, double *, double *, double *, |
| 69 |
char *, char *, int *, int *, int *, double *, double *, |
| 70 |
int *, double *, double *, double *, int *, int, int, int); |
| 71 |
|
| 72 |
char *module_name = "rdvinv"; |
| 73 |
char *version_id = "0.91"; |
| 74 |
|
| 75 |
ModuleArgs_t module_args[] = { |
| 76 |
{ARG_STRING, "in", "", "Input data series or recordset"}, |
| 77 |
{ARG_STRING, "seg", "fit.out", "Input data segment name"}, |
| 78 |
{ARG_STRING, "uxseg", "Ux", "Output data segment name for Ux inversions"}, |
| 79 |
{ARG_STRING, "uyseg", "Uy", "Output data segment name for Uy inversions"}, |
| 80 |
{ARG_INT, "cr", "Not Specified", "Carrington rotation for all regions"}, |
| 81 |
{ARG_FLOAT, "clon", "Not Specified", |
| 82 |
"Carrington longitude of CM for all regions"}, |
| 83 |
{ARG_FLOAT, "lon", "Not Specified", "Carrington longitude of all regions"}, |
| 84 |
{ARG_FLOAT, "lat", "Not Specified", "Carrington latitude of all regions"}, |
| 85 |
{ARG_DOUBLE, "amu", "0.005", "Error trade-off parameter"}, |
| 86 |
{ARG_DOUBLE, "ob", "1.0", "Lower frequency limit (mHz)"}, |
| 87 |
{ARG_DOUBLE, "oe", "5.2", "Upper frequency limit (mHz)"}, |
| 88 |
{ARG_DOUBLE, "rb", "0.97", "Lower radius limit"}, |
| 89 |
{ARG_DOUBLE, "re", "1.00", "Upper radius limit"}, |
| 90 |
{ARG_INT, "num", "40", "Number of target inversion points"}, |
| 91 |
{ARG_STRING, "out", "fort.10.hmixy", "Output filename"}, |
| 92 |
{ARG_STRING, "kernel", "", ""}, |
| 93 |
{ARG_STRING, "ave", "Not Specified", |
| 94 |
"output file for averaging kernels (if not set, kernels are not written)"}, |
| 95 |
{ARG_STRING, "coef", "Not Specified", ""}, |
| 96 |
{ARG_FLAG, "v", "", "run in verbose mode"}, |
| 97 |
{ARG_END} |
| 98 |
}; |
| 99 |
|
| 100 |
/* list of keywords to propagate (if possible) from input to output */ |
| 101 |
char *propagate[] = {"CarrTime", "CarrRot", "CMLon", "LonHG", "LatHG", "LonCM", |
| 102 |
"MidTime", "Duration", "MapProj", "MapScale", "Map_PA", "Width", "Height", |
| 103 |
"Size", "Cadence", "ZonalTrk", "ZonalVel", "MeridTrk", "MeridVel", "MAI"}; |
| 104 |
|
| 105 |
int process_record (const char *filename, int drms_output, char *outfilex, |
| 106 |
char *outfiley, char *kernel, char *kername, char *ave, char *coef, |
| 107 |
int qave, int qcoef, int verbose, double ob, double oe, int num, |
| 108 |
double rb, double re, double amu, char *sourcerec, char *codever) { |
| 109 |
FILE *infile = fopen (filename, "r"); |
| 110 |
FILE *filex = fopen (outfilex, "w"); |
| 111 |
FILE *filey = fopen (outfiley, "w"); |
| 112 |
double *rcgx, *rcgy, *quartx, *quarty, *solnx, *solny, *errx, *erry; |
| 113 |
double *l, *f, *ef, *ux, *eux, *uy, *euy, *rtrg; |
| 114 |
int *n, *mask; |
| 115 |
int npts, i, j, status; |
| 116 |
int lenkern, lenave, lencoef; |
| 117 |
status = read_fit_v (infile, &npts, &n, &l, &f, &ef, &ux, &eux, &uy, &euy); |
| 118 |
fclose (infile); |
| 119 |
if(status) { |
| 120 |
fprintf(stderr, "File %s could not be read\n", filename); |
| 121 |
return 1; |
| 122 |
} |
| 123 |
|
| 124 |
rtrg = (double *)malloc (num * sizeof (double)); |
| 125 |
rcgx = (double *)malloc (num * sizeof (double)); |
| 126 |
rcgy = (double *)malloc (num * sizeof (double)); |
| 127 |
quartx = (double *)malloc (3 * num * sizeof (double)); |
| 128 |
quarty = (double *)malloc (3 * num * sizeof (double)); |
| 129 |
solnx = (double *)malloc (num * sizeof (double)); |
| 130 |
solny = (double *)malloc (num * sizeof (double)); |
| 131 |
errx = (double *)malloc (num * sizeof (double)); |
| 132 |
erry = (double *)malloc (num * sizeof (double)); |
| 133 |
mask = (int *)malloc (npts * sizeof(int)); |
| 134 |
interp_vel (n, l, f, ef, ux, eux, uy, euy, npts); |
| 135 |
autoweed_vel (n, l, ux, uy, mask, npts); |
| 136 |
j = 0; |
| 137 |
for (i=0; i<npts; i++) { |
| 138 |
if (mask[i]) { |
| 139 |
n[j] = n[i]; |
| 140 |
l[j] = l[i]; |
| 141 |
f[j] = f[i]; |
| 142 |
ef[j] = ef[i]; |
| 143 |
ux[j] = ux[i]; |
| 144 |
eux[j] = eux[i]; |
| 145 |
uy[j] = uy[i]; |
| 146 |
euy[j] = euy[i]; |
| 147 |
j++; |
| 148 |
} |
| 149 |
} |
| 150 |
npts = j; |
| 151 |
|
| 152 |
fprintf (filex, "# Flow inversion of %s\n", sourcerec); |
| 153 |
fprintf (filey, "# Flow inversion of %s\n", sourcerec); |
| 154 |
fprintf (filex, "# against kernel %s\n", kername); |
| 155 |
fprintf (filey, "# against kernel %s\n", kername); |
| 156 |
fprintf (filex, "# %s\n", codever); |
| 157 |
fprintf (filey, "# %s\n", codever); |
| 158 |
fprintf (filex, "# amu = %.5f ob = %.3f mHz, oe = %.3f mHz\n", amu, ob, oe); |
| 159 |
fprintf (filey, "# amu = %.5f ob = %.3f mHz, oe = %.3f mHz\n", amu, ob, oe); |
| 160 |
fprintf (filex, "# rb = %.3f, re= %.3f, num = %d\n",rb, re, num); |
| 161 |
fprintf (filey, "# rb = %.3f, re= %.3f, num = %d\n",rb, re, num); |
| 162 |
fprintf (filex, |
| 163 |
"# Rtrg amu CG of av. kernel, quartiles soln err\n"); |
| 164 |
fprintf (filex, "# 1 2 3 3-1\n"); |
| 165 |
fprintf (filey, |
| 166 |
"# Rtrg amu CG of av. kernel, quartiles soln err\n"); |
| 167 |
fprintf (filey, "# 1 2 3 3-1\n"); |
| 168 |
|
| 169 |
lenkern = strlen (kernel); |
| 170 |
lenave = (qave) ? strlen (ave) : 0; |
| 171 |
lencoef = (qcoef) ? strlen (coef) : 0; |
| 172 |
|
| 173 |
ola_ (l, n, f, ux, eux, uy, euy, &npts, kernel, rtrg, rcgx, rcgy, |
| 174 |
quartx, quarty, solnx, solny, errx, erry, ave, coef, &qave, &qcoef, |
| 175 |
&verbose, &ob, &oe, &num, &rb, &re, &amu, &status, |
| 176 |
lenkern, lenave, lencoef); |
| 177 |
if (status) return status; |
| 178 |
|
| 179 |
for (i = 0; i < num; i++) |
| 180 |
fprintf (filex, "%7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %12.5e %12.5e\n", |
| 181 |
rtrg[i], amu, rcgx[i], quartx[2 + 3*i], quartx[1 + 3*i], quartx[3*i], |
| 182 |
fabs (quartx[3*i] - quartx[2+ 3*i]), solnx[i], errx[i]); |
| 183 |
for (i = 0; i < num; i++) |
| 184 |
fprintf (filey, "%7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %12.5e %12.5e\n", |
| 185 |
rtrg[i], amu, rcgy[i], quarty[2 + 3*i], quarty[1 + 3*i], quarty[3*i], |
| 186 |
fabs (quarty[3*i] - quarty[2+ 3*i]), solny[i], erry[i]); |
| 187 |
|
| 188 |
free (n); |
| 189 |
free (mask); |
| 190 |
free (l); |
| 191 |
free (f); |
| 192 |
free (ef); |
| 193 |
free (ux); |
| 194 |
free (eux); |
| 195 |
free (uy); |
| 196 |
free (euy); |
| 197 |
|
| 198 |
free (rtrg); |
| 199 |
free (rcgx); |
| 200 |
free (rcgy); |
| 201 |
free (quartx); |
| 202 |
free (quarty); |
| 203 |
free (solnx); |
| 204 |
free (solny); |
| 205 |
free (errx); |
| 206 |
free (erry); |
| 207 |
fclose (filex); |
| 208 |
fclose (filey); |
| 209 |
|
| 210 |
return 0; |
| 211 |
} |
| 212 |
/* these functions should probably be added to keystuff.c */ |
| 213 |
/* they are also in rdsinv_v07 */ |
| 214 |
int key_is_prime (DRMS_Record_t *rec, const char *key) { |
| 215 |
DRMS_Keyword_t *pkey, *keywd; |
| 216 |
int n, npct = rec->seriesinfo->pidx_num; |
| 217 |
|
| 218 |
if (!(keywd = drms_keyword_lookup (rec, key, 1))) return 0; |
| 219 |
if (npct) { |
| 220 |
for (n = 0; n < npct; n++) { |
| 221 |
pkey = rec->seriesinfo->pidx_keywords[n]; |
| 222 |
if (pkey->info->recscope > 1) pkey = drms_keyword_slotfromindex (pkey); |
| 223 |
if (strcmp (key, pkey->info->name)) continue; |
| 224 |
return 1; |
| 225 |
} |
| 226 |
} |
| 227 |
return 0; |
| 228 |
} |
| 229 |
|
| 230 |
int unique_key_values (const char *dsqry, const char *key, DRMS_Record_t *rec) { |
| 231 |
int status, values; |
| 232 |
DRMS_Array_t *value; |
| 233 |
DRMS_Keyword_t *pkey, *keywd; |
| 234 |
|
| 235 |
if (!(keywd = drms_keyword_lookup (rec, key, 1))) return 0; |
| 236 |
value = drms_record_getvector (drms_env, dsqry, key, |
| 237 |
DRMS_TYPE_FLOAT, 1, &status); |
| 238 |
if (!value) return 0; |
| 239 |
values = value->axis[1]; |
| 240 |
drms_free_array (value); |
| 241 |
return values; |
| 242 |
} |
| 243 |
|
| 244 |
char *printcrcl_loc (int cr, double cl, double lnhg, double lncm, double lat) { |
| 245 |
static char fnam[32]; |
| 246 |
double dlon = lncm; |
| 247 |
|
| 248 |
if (cr > 0) { |
| 249 |
if (isnan (cl)) { |
| 250 |
if (isnan (lncm)) { |
| 251 |
/* construct a filename CR_lonHG+lat */ |
| 252 |
sprintf (fnam, "%04d_%04.1f%+05.1f", cr, lnhg, lat); |
| 253 |
} else { |
| 254 |
/* construct a filename CR_lonEWlatNS */ |
| 255 |
if (dlon >= 0) { |
| 256 |
if (lat >= 0) sprintf (fnam, "%04d_%04.1fW%04.1fN", cr, dlon, lat); |
| 257 |
else if (lat > -0.05) sprintf (fnam, "%04d_%04.1fW%04.1fN", |
| 258 |
cr, dlon, -lat); |
| 259 |
else sprintf (fnam, "%04d_%04.1fW%04.1fS", cr, dlon, -lat); |
| 260 |
} else if (dlon > -0.05) { |
| 261 |
if (lat >= 0) sprintf (fnam, "%04d_%04.1fW%04.1fN", cr, -dlon, lat); |
| 262 |
else if (lat > -0.05) sprintf (fnam, "%04d_%04.1fW%04.1fN", |
| 263 |
cr, -dlon, -lat); |
| 264 |
else sprintf (fnam, "%04d_%04.1fW%04.1fS", cr, -dlon, -lat); |
| 265 |
} else { |
| 266 |
if (lat >= 0) sprintf (fnam, "%04d_%04.1fE%04.1fN", cr, -dlon, lat); |
| 267 |
else if (lat > -0.05) sprintf (fnam, "%04d_%04.1fE%04.1fN", |
| 268 |
cr, -dlon, -lat); |
| 269 |
else sprintf (fnam, "%04d_%04.1fE%04.1fS", cr, -dlon, -lat); |
| 270 |
} |
| 271 |
} |
| 272 |
return fnam; |
| 273 |
} |
| 274 |
/* construct a filename CR:CL_lonEWlatNS */ |
| 275 |
while (cl < 0) { |
| 276 |
cl += 360; |
| 277 |
cr++; |
| 278 |
} while (cl > 360) { |
| 279 |
cl -= 360; |
| 280 |
cr--; |
| 281 |
} |
| 282 |
if (isnan (lncm)) { |
| 283 |
dlon = lnhg - cl; |
| 284 |
while (dlon > 180) dlon -= 360; |
| 285 |
while (dlon < -180) dlon += 360; |
| 286 |
} |
| 287 |
if (dlon >= 0) { |
| 288 |
if (lat >= 0) |
| 289 |
sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, dlon, lat); |
| 290 |
else if (lat > -0.05) |
| 291 |
sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, dlon, -lat); |
| 292 |
else sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fS", cr, cl, dlon, -lat); |
| 293 |
} else if (dlon > -0.05) { |
| 294 |
if (lat >= 0) |
| 295 |
sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, -dlon, lat); |
| 296 |
else if (lat > -0.05) |
| 297 |
sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, -dlon, -lat); |
| 298 |
else sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fS", cr, cl, -dlon, -lat); |
| 299 |
} else { |
| 300 |
if (lat >= 0) |
| 301 |
sprintf (fnam, "%04d:%05.1f_%04.1fE%04.1fN", cr, cl, -dlon, lat); |
| 302 |
else if (lat > -0.05) |
| 303 |
sprintf (fnam, "%04d:%05.1f_%04.1fE%04.1fN", cr, cl, -dlon, -lat); |
| 304 |
else sprintf (fnam, "%04d:%05.1f_%04.1fE%04.1fS", cr, cl, -dlon, -lat); |
| 305 |
} |
| 306 |
} else { |
| 307 |
if (isnan (dlon)) { |
| 308 |
/* construct a filename lonHG+lat */ |
| 309 |
sprintf (fnam, "%04.1f%+05.1f", lnhg, lat); |
| 310 |
return fnam; |
| 311 |
} |
| 312 |
/* construct a filename lonEWlatNS */ |
| 313 |
if (dlon >= 0) { |
| 314 |
if (lat >= 0) sprintf (fnam, "%04.1fW%04.1fN", dlon, lat); |
| 315 |
else if (lat > -0.05) sprintf (fnam, "%04.1fW%04.1fN", dlon, -lat); |
| 316 |
else sprintf (fnam, "%04.1fW%04.1fS", dlon, -lat); |
| 317 |
} else if (dlon > -0.05) { |
| 318 |
if (lat >= 0) sprintf (fnam, "%04.1fW%04.1fN", -dlon, lat); |
| 319 |
else if (lat > -0.05) sprintf (fnam, "%04.1fW%04.1fN", -dlon, -lat); |
| 320 |
else sprintf (fnam, "%04.1fW%04.1fS", -dlon, -lat); |
| 321 |
} else { |
| 322 |
if (lat >= 0) sprintf (fnam, "%04.1fE%04.1fN", -dlon, lat); |
| 323 |
else if (lat > -0.05) sprintf (fnam, "%04.1fE%04.1fN", -dlon, -lat); |
| 324 |
else sprintf (fnam, "%04.1fE%04.1fS", -dlon, -lat); |
| 325 |
} |
| 326 |
} |
| 327 |
return fnam; |
| 328 |
} |
| 329 |
/* module body */ |
| 330 |
int DoIt(void) { |
| 331 |
CmdParams_t *params = &cmdparams; |
| 332 |
int status = 0; |
| 333 |
int drms_input, drms_output; |
| 334 |
int nrec, rec_i; |
| 335 |
DRMS_RecordSet_t *recordSet = NULL, *kern_set = NULL; |
| 336 |
DRMS_Record_t *irec, *orec; |
| 337 |
DRMS_Segment_t *segment, *oseg; |
| 338 |
FILE *fpt; |
| 339 |
double latc, loncm, lonhg, loncCar; |
| 340 |
float fval; |
| 341 |
int ival; |
| 342 |
int n; |
| 343 |
char odir[DRMS_MAXPATHLEN], filename[DRMS_MAXPATHLEN+5]; |
| 344 |
char buffer[1024], outfilex[DRMS_MAXPATHLEN], outfiley[DRMS_MAXPATHLEN]; |
| 345 |
char outdir[DRMS_MAXPATHLEN], outdirx[DRMS_MAXPATHLEN], |
| 346 |
outdiry[DRMS_MAXPATHLEN], kernfile[DRMS_MAXPATHLEN]; |
| 347 |
char *carstr; |
| 348 |
char rec_query[DRMS_MAXQUERYLEN], source[DRMS_MAXQUERYLEN]; |
| 349 |
char module_ident[64]; |
| 350 |
|
| 351 |
int twosegs = 0; |
| 352 |
|
| 353 |
char *in = strdup (params_get_str (params, "in")); |
| 354 |
int cr = params_get_int (params, "cr"); |
| 355 |
float cl = params_get_float (params, "clon"); |
| 356 |
float lat = params_get_float (params, "lat"); |
| 357 |
float lon = params_get_float (params, "lon"); |
| 358 |
char *seg = strdup (params_get_str (params, "seg")); |
| 359 |
char *uxseg = strdup (params_get_str (params, "uxseg")); |
| 360 |
char *uyseg = strdup (params_get_str (params, "uyseg")); |
| 361 |
char *kernel = strdup (params_get_str (params, "kernel")); |
| 362 |
char *out = strdup (params_get_str (params, "out")); |
| 363 |
char *ave = strdup (params_get_str (params, "ave")); |
| 364 |
char *coef = strdup (params_get_str (params, "coef")); |
| 365 |
double ob = params_get_double (params, "ob"); |
| 366 |
double oe = params_get_double (params, "oe"); |
| 367 |
double rb = params_get_double (params, "rb"); |
| 368 |
double re = params_get_double (params, "re"); |
| 369 |
double amu = params_get_double (params, "amu"); |
| 370 |
int num = params_get_int (params, "num"); |
| 371 |
int verbose = params_isflagset (params, "v"); |
| 372 |
|
| 373 |
int qave = strcmp (ave, "Not Specified"); |
| 374 |
int qcoef = strcmp (coef, "Not Specified"); |
| 375 |
int anycr = (cr <= 0); |
| 376 |
int anycl = isnan (cl); |
| 377 |
int anylon = isnan (lon); |
| 378 |
int anylat = isnan (lat); |
| 379 |
|
| 380 |
snprintf (module_ident, 64, "%s v %s", module_name, version_id); |
| 381 |
if (verbose) printf ("%s:\n", module_ident); |
| 382 |
|
| 383 |
drms_input = 1; |
| 384 |
/* attempt to open the kernel as a DRMS record set */ |
| 385 |
kern_set = drms_open_records (drms_env, kernel, &status); |
| 386 |
if (status) { |
| 387 |
/* attempt to open the kernel as a file name */ |
| 388 |
fpt = fopen (kernel, "r"); |
| 389 |
if (!fpt) { |
| 390 |
fprintf (stderr, "Error: kernel specification \"%s\"\n", kernel); |
| 391 |
fprintf (stderr, |
| 392 |
" does not specify a readable data record nor file\n"); |
| 393 |
return 1; |
| 394 |
} |
| 395 |
fclose (fpt); |
| 396 |
strcpy (kernfile, kernel); |
| 397 |
} else { |
| 398 |
/* check that the kernel set specifies unique record */ |
| 399 |
if (kern_set->n != 1) { |
| 400 |
fprintf (stderr, "Error: no data records in set %s\n", kernel); |
| 401 |
return 1; |
| 402 |
} |
| 403 |
irec = kern_set->records[0]; |
| 404 |
drms_record_directory (irec, kernfile, 1); |
| 405 |
segment = drms_segment_lookup (irec, "kernel"); |
| 406 |
if (!segment) { |
| 407 |
/* and that it contains the segment "kernel" */ |
| 408 |
fprintf (stderr, "Error: kernel record \"%s\"\n", kernel); |
| 409 |
fprintf (stderr, " does not contain a segment \"kernel\"\n"); |
| 410 |
drms_close_records (kern_set, DRMS_FREE_RECORD); |
| 411 |
return 1; |
| 412 |
} |
| 413 |
sprintf (kernfile, "%s/S%05i/kernel", segment->record->su->sudir, |
| 414 |
segment->record->slotnum); |
| 415 |
drms_close_records (kern_set, DRMS_FREE_RECORD); |
| 416 |
} |
| 417 |
|
| 418 |
if (anycr && anycl && anylon && anylat) { |
| 419 |
/* no target range specified, assume implicit in input set */ |
| 420 |
/* but it might be implicit in output series spec! */ |
| 421 |
recordSet = drms_open_records (drms_env, in, &status); |
| 422 |
if (status) { |
| 423 |
fpt = fopen (in, "r"); |
| 424 |
if (!fpt) { |
| 425 |
fprintf (stderr, "Error: in specification \"%s\"\n", in); |
| 426 |
fprintf (stderr, |
| 427 |
" does not specify a readable data set nor file\n"); |
| 428 |
return 1; |
| 429 |
} |
| 430 |
fclose (fpt); |
| 431 |
drms_input = 0; |
| 432 |
nrec = 1; |
| 433 |
} else nrec = recordSet->n; |
| 434 |
if (nrec < 1) { |
| 435 |
fprintf (stderr, "Error: no data records in set %s\n", in); |
| 436 |
return 1; |
| 437 |
} |
| 438 |
} else { |
| 439 |
strncpy (rec_query, in, DRMS_MAXQUERYLEN); |
| 440 |
if (!anycr) { |
| 441 |
snprintf (source, DRMS_MAXQUERYLEN, "[?CarrRot=%d?]", cr); |
| 442 |
strncat (rec_query, source, DRMS_MAXQUERYLEN); |
| 443 |
} |
| 444 |
if (!anycl) { |
| 445 |
float clmin = cl - 0.01; |
| 446 |
float clmax = cl + 0.01; |
| 447 |
/* need to deal with 0/360 */ |
| 448 |
if (clmin < 0.0) { |
| 449 |
clmin += 360.0; |
| 450 |
snprintf (source, DRMS_MAXQUERYLEN, "[?CMLon>%g or CMLon<%g?]", |
| 451 |
clmin, clmax); |
| 452 |
} else if (clmax > 360.0) { |
| 453 |
clmax -= 360.0; |
| 454 |
snprintf (source, DRMS_MAXQUERYLEN, "[?CMLon>%g or CMLon<%g?]", |
| 455 |
clmin, clmax); |
| 456 |
} else snprintf (source, DRMS_MAXQUERYLEN, "[?CMLon>%g and CMLon<%g?]", |
| 457 |
clmin, clmax); |
| 458 |
strncat (rec_query, source, DRMS_MAXQUERYLEN); |
| 459 |
} |
| 460 |
if (!anylon) { |
| 461 |
float lonmin = lon - 0.01; |
| 462 |
float lonmax = lon + 0.01; |
| 463 |
/* need to deal with 0/360 */ |
| 464 |
if (lonmin < 0.0) { |
| 465 |
lonmin += 360.0; |
| 466 |
snprintf (source, DRMS_MAXQUERYLEN, "[?LonHG>%g or LonHG<%g?]", |
| 467 |
lonmin, lonmax); |
| 468 |
} else if (lonmax > 360.0) { |
| 469 |
lonmax -= 360.0; |
| 470 |
snprintf (source, DRMS_MAXQUERYLEN, "[?LonHG>%g or LonHG<%g?]", |
| 471 |
lonmin, lonmax); |
| 472 |
} else snprintf (source, DRMS_MAXQUERYLEN, "[?LonHG>%g and LonHG<%g?]", |
| 473 |
lonmin, lonmax); |
| 474 |
strncat (rec_query, source, DRMS_MAXQUERYLEN); |
| 475 |
} |
| 476 |
if (!anylat) { |
| 477 |
float latmin = lat - 0.01; |
| 478 |
float latmax = lat + 0.01; |
| 479 |
snprintf (source, DRMS_MAXQUERYLEN, "[?LatHG>%g and LatHG<%g?]", |
| 480 |
latmin, latmax); |
| 481 |
strncat (rec_query, source, DRMS_MAXQUERYLEN); |
| 482 |
} |
| 483 |
if (!(recordSet = drms_open_records (drms_env, rec_query, &status))) { |
| 484 |
fprintf (stderr, "Error: unable to open input data set %s\n", rec_query); |
| 485 |
fprintf (stderr, " status = %d\n", status); |
| 486 |
return 1; |
| 487 |
} |
| 488 |
nrec = recordSet->n; |
| 489 |
} |
| 490 |
if (nrec < 1) { |
| 491 |
fprintf (stderr, "Error: no records found in input data set\n"); |
| 492 |
fprintf (stderr, " %s\n", rec_query); |
| 493 |
} |
| 494 |
if (drms_input) { |
| 495 |
/* check that range is restricted by dataset specification */ |
| 496 |
irec = recordSet->records[0]; |
| 497 |
if (anycr) { |
| 498 |
anycr = unique_key_values (in, "CarrRot", irec) - 1; |
| 499 |
ival = drms_getkey_int (irec, "CarrRot", &status); |
| 500 |
if (!status) cr = ival; |
| 501 |
} |
| 502 |
if (anycl) { |
| 503 |
anycl = unique_key_values (in, "CMLon", irec) - 1; |
| 504 |
fval = drms_getkey_float (irec, "CMLon", &status); |
| 505 |
if (!status) cl = fval; |
| 506 |
} |
| 507 |
if (anylon) { |
| 508 |
anylon = unique_key_values (in, "LonHG", irec) - 1; |
| 509 |
fval = drms_getkey_float (irec, "LonHG", &status); |
| 510 |
if (!status) lon = fval; |
| 511 |
} |
| 512 |
if (anylat) { |
| 513 |
anylat = unique_key_values (in, "LatHG", irec) - 1; |
| 514 |
fval = drms_getkey_float (irec, "LatHG", &status); |
| 515 |
if (!status) lat = fval; |
| 516 |
} |
| 517 |
if (verbose) printf ("Processing %d records\n", nrec); |
| 518 |
} |
| 519 |
/* Attempt to create output drms record */ |
| 520 |
drms_output = 1; |
| 521 |
orec = drms_create_record (drms_env, out, DRMS_PERMANENT, &status); |
| 522 |
if (status) { |
| 523 |
drms_output = 0; |
| 524 |
fprintf (stderr, |
| 525 |
"Warning: drms_create_record() returned %d for data series:\n", status); |
| 526 |
fprintf (stderr, |
| 527 |
" %s\n will be interpreted as directory name\n", out); |
| 528 |
strncpy (odir, out, DRMS_MAXPATHLEN); |
| 529 |
} else if (drms_record_directory (orec, odir, 1)) { |
| 530 |
fprintf (stderr, |
| 531 |
"Error: drms_record_directory() failure: SUMS may not be available\n"); |
| 532 |
fprintf (stderr, |
| 533 |
" or series %s may contain no segments\n", out); |
| 534 |
return 1; |
| 535 |
} |
| 536 |
if (drms_output) { |
| 537 |
int kstat = 0; |
| 538 |
int keyct = sizeof (propagate) / sizeof (char *); |
| 539 |
|
| 540 |
kstat += check_and_set_key_str (orec, "Module", module_ident); |
| 541 |
kstat += check_and_set_key_time (orec, "Created", CURRENT_SYSTEM_TIME); |
| 542 |
/* module specific keywords */ |
| 543 |
kstat += check_and_set_key_float (orec, "amu", amu); |
| 544 |
kstat += check_and_set_key_float (orec, "freqmin", ob); |
| 545 |
kstat += check_and_set_key_float (orec, "freqmax", oe); |
| 546 |
kstat += check_and_set_key_float (orec, "radmin", rb); |
| 547 |
kstat += check_and_set_key_float (orec, "radmax", re); |
| 548 |
kstat += check_and_set_key_str (orec, "BLD_VERS", jsoc_version); |
| 549 |
/* set possible prime keys as appropriate */ |
| 550 |
if (anycr) { |
| 551 |
if (key_is_prime (orec, "CarrRot")) { |
| 552 |
fprintf (stderr, |
| 553 |
"Error: CarrRot is prime key for output series %s\n", out); |
| 554 |
fprintf (stderr, |
| 555 |
" but input data set contains multiple values\n"); |
| 556 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 557 |
/* this pretty much has to be true... */ |
| 558 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 559 |
return 1; |
| 560 |
} |
| 561 |
} else kstat += check_and_set_key_int (orec, "CarrRot", cr); |
| 562 |
if (anycl) { |
| 563 |
if (key_is_prime (orec, "CMLon")) { |
| 564 |
fprintf (stderr, |
| 565 |
"Error: CMLon is prime key for output series %s\n", out); |
| 566 |
fprintf (stderr, |
| 567 |
" but input data set contains multiple values\n"); |
| 568 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 569 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 570 |
return 1; |
| 571 |
} |
| 572 |
} else kstat += check_and_set_key_double (orec, "CMLon", cl); |
| 573 |
if (anycr || anycl) { |
| 574 |
if (key_is_prime (orec, "CarrTime")) { |
| 575 |
fprintf (stderr, |
| 576 |
"Error: CarrTime is prime key for output series %s\n", out); |
| 577 |
fprintf (stderr, |
| 578 |
" but input data set contains multiple values\n"); |
| 579 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 580 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 581 |
return 1; |
| 582 |
} |
| 583 |
} else { |
| 584 |
char crcl[32]; |
| 585 |
sprintf (crcl, "%d:%03.0f", cr, cl); |
| 586 |
kstat += check_and_set_key_str (orec, "CarrTime", crcl); |
| 587 |
} |
| 588 |
if (anylon) { |
| 589 |
if (key_is_prime (orec, "LonHG")) { |
| 590 |
fprintf (stderr, |
| 591 |
"Error: LonHG is prime key for output series %s\n", out); |
| 592 |
fprintf (stderr, |
| 593 |
" but input data set contains multiple values\n"); |
| 594 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 595 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 596 |
return 1; |
| 597 |
} |
| 598 |
} else kstat += check_and_set_key_double (orec, "LonHG", lon); |
| 599 |
if (anylat) { |
| 600 |
if (key_is_prime (orec, "LatHG")) { |
| 601 |
fprintf (stderr, |
| 602 |
"Error: LatHG is prime key for output series %s\n", out); |
| 603 |
fprintf (stderr, |
| 604 |
" but input data set contains multiple values\n"); |
| 605 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 606 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 607 |
return 1; |
| 608 |
} |
| 609 |
} else kstat += check_and_set_key_double (orec, "LatHG", lat); |
| 610 |
|
| 611 |
if (kstat) { |
| 612 |
fprintf (stderr, "Error writing key value(s) to record in series %s:\n", |
| 613 |
out); |
| 614 |
fprintf (stderr, " series may not have appropriate structure;\n"); |
| 615 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 616 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 617 |
return 1; |
| 618 |
} |
| 619 |
/* search for separate output segments ux and uy */ |
| 620 |
if (drms_segment_lookup (orec, uxseg) && drms_segment_lookup (orec, uyseg)) |
| 621 |
twosegs = 1; |
| 622 |
else { |
| 623 |
/* search for 1st generic segment and use it */ |
| 624 |
int segct = drms_record_numsegments (orec); |
| 625 |
if (segct < 1) { |
| 626 |
fprintf (stderr, |
| 627 |
"Error: no data segment in output series %s\n", out); |
| 628 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 629 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 630 |
return 1; |
| 631 |
} |
| 632 |
for (n = 0; n < segct; n++) { |
| 633 |
oseg = drms_segment_lookupnum (orec, n); |
| 634 |
if (oseg->info->protocol != DRMS_GENERIC) continue; |
| 635 |
break; |
| 636 |
} |
| 637 |
if (n == segct) { |
| 638 |
fprintf (stderr, |
| 639 |
"Error: no generic data segment in output series %s\n", out); |
| 640 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 641 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 642 |
return 1; |
| 643 |
} |
| 644 |
fprintf (stderr, " writing to segment %s\n", oseg->info->name); |
| 645 |
} |
| 646 |
} |
| 647 |
if (twosegs) { |
| 648 |
sprintf (outdirx, "%s/%s", odir, uxseg); |
| 649 |
sprintf (outdiry, "%s/%s", odir, uyseg); |
| 650 |
mkdir (outdirx, 01755); |
| 651 |
mkdir (outdiry, 01755); |
| 652 |
} else { |
| 653 |
if (drms_output) { |
| 654 |
sprintf (outdir, "%s/%s", odir, oseg->info->name); |
| 655 |
mkdir (outdir, 01755); |
| 656 |
} else strcpy (outdir, out); |
| 657 |
} |
| 658 |
/* main processing loop */ |
| 659 |
for (rec_i=0; rec_i<nrec; rec_i++) { |
| 660 |
if (verbose) printf (" Processing record %i\n", rec_i); |
| 661 |
if (drms_input) { |
| 662 |
irec = recordSet->records[rec_i]; |
| 663 |
/* should use call to drms_record_directory() instead */ |
| 664 |
if (irec->sunum != -1LL && irec->su == NULL) { |
| 665 |
irec->su = drms_getunit (irec->env, irec->seriesinfo->seriesname, |
| 666 |
irec->sunum, 1, &status); |
| 667 |
} |
| 668 |
segment = drms_segment_lookup(irec, seg); |
| 669 |
sprintf (filename, "%s/S%05i/%s", segment->record->su->sudir, |
| 670 |
segment->record->slotnum, seg); |
| 671 |
drms_sprint_rec_query (source, irec); |
| 672 |
cr = drms_getkey_int (irec, "CarrRot", &status); |
| 673 |
latc = drms_getkey_double (irec, "LatHG", &status); |
| 674 |
lonhg = drms_getkey_double (irec, "LonHG", &status); |
| 675 |
carstr = drms_getkey_string (irec, "CarrTime", &status); |
| 676 |
if (status) loncCar = drms_getkey_double (irec, "CMLon", &status); |
| 677 |
else sscanf (carstr, "%i:%lf", &cr, &loncCar); |
| 678 |
loncm = drms_getkey_double (irec, "LonCM", &status); |
| 679 |
} else { |
| 680 |
strcpy (filename, in); |
| 681 |
strcpy (source, in); |
| 682 |
latc = loncm = lonhg = 0.0; |
| 683 |
} |
| 684 |
if (twosegs) { |
| 685 |
sprintf (outfilex, "%s/%s/%s.Ux", odir, uxseg, |
| 686 |
printcrcl_loc (cr, loncCar, lonhg, loncm, latc)); |
| 687 |
sprintf (outfiley, "%s/%s/%s.Uy", odir, uyseg, |
| 688 |
printcrcl_loc (cr, loncCar, lonhg, loncm, latc)); |
| 689 |
} else { |
| 690 |
sprintf (outfilex, "%s/%s.Ux", outdir, |
| 691 |
printcrcl_loc (cr, loncCar, lonhg, loncm, latc)); |
| 692 |
sprintf (outfiley, "%s/%s.Uy", outdir, |
| 693 |
printcrcl_loc (cr, loncCar, lonhg, loncm, latc)); |
| 694 |
} |
| 695 |
status = process_record (filename, drms_output, outfilex, outfiley, |
| 696 |
kernfile, kernel, ave, coef, qave, qcoef, verbose, ob, oe, num, rb, re, amu, |
| 697 |
source, module_ident); |
| 698 |
if (status) { |
| 699 |
fprintf (stderr, "Error processing record %d (%s); aborted\n", rec_i, |
| 700 |
printcrcl_loc (cr, loncCar, lonhg, loncm, latc)); |
| 701 |
if (drms_output) { |
| 702 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 703 |
drms_close_record (orec, DRMS_FREE_RECORD); |
| 704 |
return 1; |
| 705 |
} |
| 706 |
} |
| 707 |
} |
| 708 |
|
| 709 |
if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD); |
| 710 |
if (drms_output) drms_close_record (orec, DRMS_INSERT_RECORD); |
| 711 |
|
| 712 |
return 0; |
| 713 |
} |
| 714 |
|
| 715 |
/* |
| 716 |
* Revision History |
| 717 |
* |
| 718 |
* June 2009 - v 0.1 module created, no real drms functionality |
| 719 |
* 09.10.2009 - v 0.2 drms functionality added - input and output from drms |
| 720 |
* records now possible, but input/output to rooted filenames also |
| 721 |
* possible |
| 722 |
* 0.3 2009.10.22 fixed format of name of output files |
| 723 |
* 0.4 2010.02.11 unified processing calls for named file and drms input |
| 724 |
* branches |
| 725 |
* 0.5 2010.02.12 Combined branches into one, added some very obvious |
| 726 |
* error handling (no given input argument, input cannot be opened) |
| 727 |
* removed rdutil.h inclusion |
| 728 |
* v 0.5 frozen 2010.03.08 |
| 729 |
* 0.6 Added logic for specifying input grouping of data set records |
| 730 |
* by Carr Rot, CM Longitude, Longitude, and/or Latitude |
| 731 |
* Added keyword propagation list (but not propagation!) |
| 732 |
* Output to multiple files named by lat and lon in segment |
| 733 |
* directory |
| 734 |
* Fixed processing of longitude queries near 0/360 |
| 735 |
* 2010.03.10 Fixed setting of keys for run parameters; put in trap |
| 736 |
* for empty input data sets and widened clon match range |
| 737 |
* v 0.6 frozen 2010.04.23 |
| 738 |
* 0.7 Output data files are now written into per segment |
| 739 |
* subdirectories; added option for naming individual segments |
| 740 |
* (if there are two) |
| 741 |
* v 0.7 frozen 2010.08.19 |
| 742 |
* 0.8 Created 2011.04.22; Fixed typos in ola_xy.c (not yet used) |
| 743 |
* Fixed bugs in printing of filenames (particularly the |
| 744 |
* reversal of east and west), and dropped requirement for CMLon; |
| 745 |
* heliographic option for file naming |
| 746 |
* Removed default for kernel |
| 747 |
* 0.9 Created 2011.12.04; Fixed typos in ola_xy.c (not yet used) |
| 748 |
* v 0.8 frozen 2012.02.01 (but included in JSOC release 6.1 2011.12.05 in |
| 749 |
* identical form except for comments) |
| 750 |
* 0.9 added (preferred) option for specifying kernel as DRMS data |
| 751 |
* record rather than filename |
| 752 |
* v 0.9 frozen 2012.04.24 |
| 753 |
* 0.91 similar to 0.9, but includes fixes in FORTRAN code ola_xy_v12 |
| 754 |
* and provision for return status |
| 755 |
* v 0.91 frozen 2012.04.24 |
| 756 |
*/ |
