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