1 // $Id: HarmonicsFile.cc 7306 2020-06-28 00:24:42Z flaterco $
2
3 /* HarmonicsFile Hide details of interaction with libtcd.
4
5 Although a few method signatures elsewhere in XTide have been
6 optimized for compatibility with libtcd, at least 95% of the work
7 of integrating XTide with a different database should be
8 concentrated in re-implementing this class.
9
10 Copyright (C) 1998 David Flater w/ contributions by Jan Depner.
11
12 This program is free software: you can redistribute it and/or modify
13 it under the terms of the GNU General Public License as published by
14 the Free Software Foundation, either version 3 of the License, or
15 (at your option) any later version.
16
17 This program is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 GNU General Public License for more details.
21
22 You should have received a copy of the GNU General Public License
23 along with this program. If not, see <http://www.gnu.org/licenses/>.
24 */
25
26 #include "libxtide.hh"
27 #include "HarmonicsFile.hh"
28 #include "SubordinateStation.hh"
29 #include <tcd.h>
30
31 namespace libxtide {
32
33
34 static bool haveInstance = false;
35
36
HarmonicsFile(const Dstr & filename)37 HarmonicsFile::HarmonicsFile (const Dstr &filename):
38 _filename(filename) {
39
40 // libtcd is stateful and cannot handle multiple harmonics files
41 // simultaneously. XTide currently has no need to open multiple
42 // harmonics files simultaneously, so for now, the constructor will
43 // barf if the attempt is made to have more than one instance at a
44 // time. If this class is modified to deal with different
45 // databases, that trap can be removed.
46 assert (!haveInstance);
47 haveInstance = true;
48
49 // Do some sanity checks before invoking libtcd. (open_tide_db just
50 // returns false regardless of what the problem was.)
51 {
52 FILE *fp = fopen (filename.aschar(), "rb");
53 if (fp) {
54 char c = fgetc (fp);
55 if (c != '[') {
56 Dstr details (filename);
57 details += " is apparently not a TCD file.\n\
58 We do not use harmonics.txt or offsets.xml anymore. Please see\n\
59 https://flaterco.com/xtide/files.html for a link to the current data.";
60 Global::barf (Error::CORRUPT_HARMONICS_FILE, details);
61 }
62 fclose (fp);
63 } else {
64 // This should never happen. We statted this file in Global.cc.
65 Global::cantOpenFile (filename, Error::fatal);
66 }
67 }
68
69 if (!open_tide_db (_filename.aschar())) {
70 Dstr details (_filename);
71 details += ": libtcd returned generic failure.";
72 Global::barf (Error::CORRUPT_HARMONICS_FILE, details);
73 }
74 DB_HEADER_PUBLIC db = get_tide_db_header();
75 _versionString = _filename;
76 {
77 int i;
78 if ((i = _versionString.strrchr ('/')) != -1)
79 _versionString /= (i+1);
80 }
81 _versionString += ' ';
82 _versionString += (char *) db.last_modified;
83 _versionString += " / ";
84 _versionString += (char *) db.version;
85 }
86
87
versionString()88 const Dstr &HarmonicsFile::versionString() {
89 return _versionString;
90 }
91
92
93 // libtcd's open_tide_db is called each time a HarmonicsFile is
94 // created. Intentionally, there is no matching call to
95 // close_tide_db. See https://flaterco.com/xtide/tcd_notes.html
96 // and XTide changelog for version 2.6.3.
97
~HarmonicsFile()98 HarmonicsFile::~HarmonicsFile() {
99 assert (haveInstance);
100 haveInstance = false;
101 }
102
103
getNextStationRef()104 StationRef * const HarmonicsFile::getNextStationRef () {
105 TIDE_STATION_HEADER rec;
106 long i;
107 if ((i = get_next_partial_tide_record (&rec)) == -1) return NULL;
108 assert (i >= 0);
109
110 // FIXME: currents kludge
111 // We don't get units in the libtcd header (partial tide record), so this
112 // is a kludge that depends on the naming convention.
113 char *name = (char *)rec.name;
114 bool isCurrent (false);
115 size_t l = strlen (name);
116 if ((l >= 8) &&
117 ((strstr (name, " Current") == name + l - 8) ||
118 (strstr (name, " Current "))))
119 isCurrent = true;
120
121 StationRef *sr = new StationRef (_filename,
122 i,
123 name,
124 ((rec.latitude != 0.0 || rec.longitude != 0.0) ?
125 Coordinates(rec.latitude, rec.longitude) : Coordinates()),
126 (char *)get_tzfile(rec.tzfile),
127 (rec.record_type == REFERENCE_STATION),
128 isCurrent);
129 return sr;
130 }
131
132
parse_xfields(MetaFieldVector & metadata,constCharPointer xfields)133 static void parse_xfields (MetaFieldVector &metadata,
134 constCharPointer xfields) {
135 assert (xfields);
136 Dstr x (xfields);
137 Dstr linebuf, name, value;
138 x.getline (linebuf);
139 while (!(linebuf.isNull())) {
140 if (linebuf[0] == ' ') {
141 if (!(name.isNull())) {
142 linebuf /= 1;
143 value += '\n';
144 value += linebuf;
145 }
146 } else {
147 if (!(name.isNull())) {
148 metadata.push_back (MetaField (name, value));
149 name = (char *)NULL;
150 value = (char *)NULL;
151 }
152 int i = linebuf.strchr (':');
153 if (i > 0) {
154 name = linebuf;
155 name -= (unsigned)i;
156 value = linebuf.ascharfrom((unsigned)i+1);
157 }
158 }
159 x.getline (linebuf);
160 }
161 if (!(name.isNull()))
162 metadata.push_back (MetaField (name, value));
163 }
164
165
166 // Analogous to datum, levelAdds for hydraulic currents are always in
167 // knots not knots^2. The database can specify either knots or
168 // knots^2 in the sub station record, but there is only one sensible
169 // interpretation.
levelAddUnits(const TIDE_RECORD & rec)170 static const Units::PredictionUnits levelAddUnits (const TIDE_RECORD &rec) {
171 return Units::flatten (Units::parse (get_level_units(rec.level_units)));
172 }
173
174
appendOffsetsMetadata(MetaFieldVector & metadata,const TIDE_RECORD & rec)175 static void appendOffsetsMetadata (MetaFieldVector &metadata,
176 const TIDE_RECORD &rec) {
177 char tmp[80];
178 Units::PredictionUnits lu = levelAddUnits(rec);
179
180 metadata.push_back (MetaField ("Min time add",
181 rec.min_time_add ? ret_time_neat (rec.min_time_add) : "NULL"));
182 sprintf (tmp, "%+2.2f %s", rec.min_level_add, Units::shortName(lu));
183 metadata.push_back (MetaField ("Min level add",
184 rec.min_level_add ? tmp : "NULL"));
185 sprintf (tmp, "%0.3f", rec.min_level_multiply);
186 metadata.push_back (MetaField ("Min level mult",
187 rec.min_level_multiply > 0.0 ? tmp : "NULL"));
188
189 metadata.push_back (MetaField ("Max time add",
190 rec.max_time_add ? ret_time_neat (rec.max_time_add) : "NULL"));
191 sprintf (tmp, "%+2.2f %s", rec.max_level_add, Units::shortName(lu));
192 metadata.push_back (MetaField ("Max level add",
193 rec.max_level_add ? tmp : "NULL"));
194 sprintf (tmp, "%0.3f", rec.max_level_multiply);
195 metadata.push_back (MetaField ("Max level mult",
196 rec.max_level_multiply > 0.0 ? tmp : "NULL"));
197
198 if (Units::isCurrent(lu)) {
199 metadata.push_back (MetaField ("Flood begins",
200 rec.flood_begins == NULLSLACKOFFSET ? "NULL"
201 : ret_time_neat (rec.flood_begins)));
202 metadata.push_back (MetaField ("Ebb begins",
203 rec.ebb_begins == NULLSLACKOFFSET ? "NULL"
204 : ret_time_neat (rec.ebb_begins)));
205 }
206 }
207
208
209 // refStationNativeUnits: determination of hydraulic vs. regular
210 // current is made based on units at reference station. Units are
211 // flattened for levelAdd offsets.
buildMetadata(const StationRef & sr,MetaFieldVector & metadata,const TIDE_RECORD & rec,Units::PredictionUnits refStationNativeUnits,CurrentBearing minCurrentBearing,CurrentBearing maxCurrentBearing)212 static void buildMetadata (const StationRef &sr,
213 MetaFieldVector &metadata,
214 const TIDE_RECORD &rec,
215 Units::PredictionUnits refStationNativeUnits,
216 CurrentBearing minCurrentBearing,
217 CurrentBearing maxCurrentBearing) {
218 Dstr tmpbuf;
219
220 metadata.push_back (MetaField ("Name", sr.name));
221 metadata.push_back (MetaField ("In file", sr.harmonicsFileName));
222 if (rec.legalese)
223 metadata.push_back (MetaField ("Legalese", get_legalese(rec.legalese)));
224 if (rec.station_id_context[0])
225 metadata.push_back (MetaField ("Station ID context",
226 rec.station_id_context));
227 if (rec.station_id[0])
228 metadata.push_back (MetaField ("Station ID", rec.station_id));
229 if (rec.date_imported)
230 metadata.push_back (MetaField ("Date imported",
231 ret_date(rec.date_imported)));
232 sr.coordinates.print (tmpbuf);
233 metadata.push_back (MetaField ("Coordinates", tmpbuf));
234 metadata.push_back (MetaField ("Country", get_country(rec.country)));
235 metadata.push_back (MetaField ("Time zone", sr.timezone));
236 metadata.push_back (MetaField ("Native units",
237 Units::longName(refStationNativeUnits)));
238 if (!(maxCurrentBearing.isNull())) {
239 maxCurrentBearing.print (tmpbuf);
240 metadata.push_back (MetaField ("Flood direction", tmpbuf));
241 }
242 if (!(minCurrentBearing.isNull())) {
243 minCurrentBearing.print (tmpbuf);
244 metadata.push_back (MetaField ("Ebb direction", tmpbuf));
245 }
246 if (rec.source[0])
247 metadata.push_back (MetaField ("Source", rec.source));
248 metadata.push_back (MetaField ("Restriction",
249 get_restriction(rec.restriction)));
250 if (rec.comments[0])
251 metadata.push_back (MetaField ("Comments", rec.comments));
252 if (rec.notes[0])
253 metadata.push_back (MetaField ("Notes", rec.notes));
254 parse_xfields (metadata, rec.xfields);
255
256 switch (rec.header.record_type) {
257 case REFERENCE_STATION:
258 metadata.push_back (MetaField ("Type",
259 (Units::isCurrent(refStationNativeUnits) ?
260 (Units::isHydraulicCurrent(refStationNativeUnits) ?
261 "Reference station, hydraulic current" :
262 "Reference station, current")
263 : "Reference station, tide")));
264 metadata.push_back (MetaField ("Meridian",
265 ret_time_neat(rec.zone_offset)));
266 if (!Units::isCurrent(refStationNativeUnits))
267 metadata.push_back (MetaField ("Datum", get_datum(rec.datum)));
268 if (rec.months_on_station)
269 metadata.push_back (MetaField ("Months on station",
270 rec.months_on_station));
271 if (rec.last_date_on_station)
272 metadata.push_back (MetaField ("Last date on station",
273 ret_date(rec.last_date_on_station)));
274 if (rec.expiration_date)
275 metadata.push_back (MetaField ("Expiration",
276 ret_date(rec.expiration_date)));
277 metadata.push_back (MetaField ("Confidence", rec.confidence));
278 break;
279
280 case SUBORDINATE_STATION:
281 metadata.push_back (MetaField ("Type",
282 (Units::isCurrent(refStationNativeUnits) ?
283 (Units::isHydraulicCurrent(refStationNativeUnits) ?
284 "Subordinate station, hydraulic current" :
285 "Subordinate station, current")
286 : "Subordinate station, tide")));
287 metadata.push_back (MetaField ("Reference",
288 get_station(rec.header.reference_station)));
289 appendOffsetsMetadata (metadata, rec);
290 break;
291
292 default:
293 assert (false);
294 }
295 }
296
297
298 // Get constituents from a TIDE_RECORD, adjusting if needed.
getConstituents(const TIDE_RECORD & rec,SimpleOffsets adjustments)299 static const ConstituentSet getConstituents (const TIDE_RECORD &rec,
300 SimpleOffsets adjustments) {
301 assert (rec.header.record_type == REFERENCE_STATION);
302
303 DB_HEADER_PUBLIC db = get_tide_db_header();
304 Units::PredictionUnits amp_units
305 (Units::parse (get_level_units (rec.level_units)));
306 SafeVector<Constituent> constituents;
307
308 // Null constituents are eliminated here.
309 for (unsigned i=0; i<db.constituents; ++i)
310 if (rec.amplitude[i] > 0.0)
311 constituents.push_back (
312 Constituent (get_speed(i),
313 db.start_year,
314 db.number_of_years,
315 get_equilibriums(i),
316 get_node_factors(i),
317 Amplitude (amp_units, rec.amplitude[i]),
318 rec.epoch[i]));
319
320 assert (!constituents.empty());
321
322 PredictionValue datum (Units::flatten(amp_units), rec.datum_offset);
323
324 // Normalize the meridian to UTC.
325 // To compensate for a negative meridian requires a positive offset.
326 // (This adjustment is combined with any that were passed in.)
327
328 // This is the only place where mutable offsets would make even a
329 // little bit of sense.
330 adjustments = SimpleOffsets (
331 adjustments.timeAdd() - Interval(ret_time_neat(rec.zone_offset)),
332 adjustments.levelAdd(),
333 adjustments.levelMultiply());
334
335 ConstituentSet cs (constituents, datum, adjustments);
336
337 Dstr u (Global::settings["u"].s);
338 if (u != "x")
339 cs.setUnits (Units::parse (u));
340
341 return cs;
342 }
343
344
getTideRecord(uint32_t recordNumber,TIDE_RECORD & rec)345 static void getTideRecord (uint32_t recordNumber, TIDE_RECORD &rec) {
346 require (read_tide_record ((NV_INT32)recordNumber, &rec)
347 == (NV_INT32)recordNumber);
348 if (Global::settings["in"].c == 'y' &&
349 rec.header.record_type == REFERENCE_STATION)
350 infer_constituents (&rec);
351 }
352
353
getStation(const StationRef & stationRef)354 Station * const HarmonicsFile::getStation (const StationRef &stationRef) {
355 Station *s = NULL;
356 TIDE_RECORD rec;
357 getTideRecord (stationRef.recordNumber, rec);
358
359 Dstr note;
360 if (rec.notes[0])
361 note = rec.notes;
362 CurrentBearing minCurrentBearing, maxCurrentBearing;
363 bool isDegreesTrue = (!(strcmp((char *)get_dir_units(rec.direction_units),
364 "degrees true")));
365 if (rec.min_direction >= 0 && rec.min_direction < 360)
366 minCurrentBearing = CurrentBearing (rec.min_direction, isDegreesTrue);
367 if (rec.max_direction >= 0 && rec.max_direction < 360)
368 maxCurrentBearing = CurrentBearing (rec.max_direction, isDegreesTrue);
369 Dstr name ((char *)rec.header.name);
370 if (rec.legalese) {
371 name += " - ";
372 name += get_legalese(rec.legalese);
373 }
374
375 switch (rec.header.record_type) {
376 case REFERENCE_STATION:
377 {
378 MetaFieldVector metadata;
379 buildMetadata (stationRef,
380 metadata,
381 rec,
382 Units::parse(get_level_units(rec.level_units)),
383 minCurrentBearing,
384 maxCurrentBearing);
385 s = new Station (name,
386 stationRef,
387 getConstituents (rec, SimpleOffsets()),
388 note,
389 minCurrentBearing,
390 maxCurrentBearing,
391 metadata);
392 }
393 break;
394
395 case SUBORDINATE_STATION:
396 {
397 TIDE_RECORD referenceStationRec;
398 assert (rec.header.reference_station >= 0);
399 getTideRecord (rec.header.reference_station, referenceStationRec);
400 Units::PredictionUnits refStationNativeUnits
401 (Units::parse(get_level_units(referenceStationRec.level_units)));
402
403 MetaFieldVector metadata;
404 buildMetadata (stationRef,
405 metadata,
406 rec,
407 refStationNativeUnits,
408 minCurrentBearing,
409 maxCurrentBearing);
410
411 NullableInterval tempFloodBegins, tempEbbBegins;
412
413 // For these, zero is not the same as null.
414 if (rec.flood_begins != NULLSLACKOFFSET)
415 tempFloodBegins = Interval (ret_time_neat (rec.flood_begins));
416 if (rec.ebb_begins != NULLSLACKOFFSET)
417 tempEbbBegins = Interval (ret_time_neat (rec.ebb_begins));
418
419 Units::PredictionUnits lu = levelAddUnits(rec);
420 HairyOffsets ho (
421 SimpleOffsets (Interval(ret_time_neat(rec.max_time_add)),
422 PredictionValue(lu, rec.max_level_add),
423 rec.max_level_multiply),
424 SimpleOffsets (Interval(ret_time_neat(rec.min_time_add)),
425 PredictionValue(lu, rec.min_level_add),
426 rec.min_level_multiply),
427 tempFloodBegins,
428 tempEbbBegins);
429
430 // If the offsets can be reduced to simple, then we can adjust
431 // the constituents and be done with it. However, in the case
432 // of hydraulic currents, level multipliers cannot be treated as
433 // simple and applied to the constants because of the square
434 // root operation--i.e., it's nonlinear.
435
436 SimpleOffsets so;
437 if ((!Units::isHydraulicCurrent(refStationNativeUnits) ||
438 ho.maxLevelMultiply() == 1.0) && ho.trySimplify (so))
439
440 // Can simplify.
441 s = new Station (name,
442 stationRef,
443 getConstituents (referenceStationRec, so),
444 note,
445 minCurrentBearing,
446 maxCurrentBearing,
447 metadata);
448
449 else
450
451 // Cannot simplify.
452 s = new SubordinateStation (name,
453 stationRef,
454 getConstituents (referenceStationRec,
455 SimpleOffsets()),
456 note,
457 minCurrentBearing,
458 maxCurrentBearing,
459 metadata,
460 ho);
461 }
462 break;
463
464 default:
465 assert (false);
466 }
467
468 // If this sanity check is deferred until Graph::draw it has
469 // unhealthy consequences (clock windows can be killed while only
470 // partially constructed; graph windows can double-barf due to
471 // events queued on the X server).
472 if (s->minLevelHeuristic() >= s->maxLevelHeuristic())
473 Global::barf (Error::ABSURD_OFFSETS, s->name);
474
475 return s;
476 }
477
478 }
479