1 // $Id: Station.cc 5748 2014-10-11 19:38:53Z flaterco $
2
3 /* Station A tide station.
4
5 Station has a subclass SubordinateStation. The superclass is used
6 for reference stations and that rare subordinate station where the
7 offsets can be reduced to simple corrections to the constituents
8 and datum. After such corrections are made, there is no
9 operational difference between that and a reference station.
10
11 Copyright (C) 1998 David Flater.
12
13 This program is free software: you can redistribute it and/or modify
14 it under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 This program is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU General Public License for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with this program. If not, see <http://www.gnu.org/licenses/>.
25 */
26
27 #include "libxtide.hh"
28 #include "Calendar.hh"
29 #include "Graph.hh"
30 #include "PixelatedGraph.hh"
31 #include "TTYGraph.hh"
32 #include "Banner.hh"
33 #include "RGBGraph.hh"
34 #include "Skycal.hh"
35 #include "SVGGraph.hh"
36
37 #ifdef HAVE_SYS_RESOURCE_H
38 #include <sys/resource.h>
39 #endif
40
41 namespace libxtide {
42
43
44 #ifndef HAVE_LLROUND
45 // Round to nearest integer, away from zero.
46 // (Interval_rep_t is what we want; long long int is possibly longer.)
llround(double x)47 static interval_rep_t llround (double x) {
48 interval_rep_t ret;
49 if (x < 0)
50 ret = x - .5;
51 else
52 ret = x + .5;
53 return ret;
54 }
55 #endif
56
57
Station(const Dstr & name_,const StationRef & stationRef,const ConstituentSet & constituents,const Dstr & note_,CurrentBearing minCurrentBearing_,CurrentBearing maxCurrentBearing_,const MetaFieldVector & metadata)58 Station::Station (const Dstr &name_,
59 const StationRef &stationRef,
60 const ConstituentSet &constituents,
61 const Dstr ¬e_,
62 CurrentBearing minCurrentBearing_,
63 CurrentBearing maxCurrentBearing_,
64 const MetaFieldVector &metadata):
65 name(name_),
66 coordinates(stationRef.coordinates),
67 timezone(stationRef.timezone),
68 minCurrentBearing(minCurrentBearing_),
69 maxCurrentBearing(maxCurrentBearing_),
70 note(note_),
71 isCurrent(Units::isCurrent(constituents.predictUnits())),
72 aspect(Global::settings["ga"].d),
73 step(Global::hour),
74 _stationRef(stationRef),
75 _constituents(constituents),
76 _metadata(metadata) {}
77
78
~Station()79 Station::~Station() {}
80
81
clone() const82 Station * const Station::clone() const {
83 return new Station (*this);
84 }
85
86
reload() const87 Station * const Station::reload() const {
88 Station *s = _stationRef.load();
89 s->markLevel = markLevel;
90 if (!markLevel.isNull())
91 if (markLevel.Units() != s->predictUnits())
92 s->markLevel.Units (s->predictUnits());
93 s->step = step;
94 return s;
95 }
96
97
minLevelHeuristic() const98 const PredictionValue Station::minLevelHeuristic() const {
99 return _constituents.datum() - _constituents.maxAmplitudeHeuristic();
100 }
101
102
maxLevelHeuristic() const103 const PredictionValue Station::maxLevelHeuristic() const {
104 return _constituents.datum() + _constituents.maxAmplitudeHeuristic();
105 }
106
107
isSubordinateStation()108 const bool Station::isSubordinateStation() {
109 return false;
110 }
111
112
haveFloodBegins()113 const bool Station::haveFloodBegins() {
114 return true;
115 }
116
117
haveEbbBegins()118 const bool Station::haveEbbBegins() {
119 return true;
120 }
121
122
123 // The following block of methods is slightly revised from the code
124 // delivered by Geoffrey T. Dairiki for XTide 1. Jeff's original
125 // comments (modulo a few global replacements) are shown in C-style,
126 // while mine are in C++ style. As usual, see also Station.hh.
127
128 /*************************************************************************
129 *
130 * Geoffrey T. Dairiki Fri Jul 19 15:44:21 PDT 1996
131 *
132 ************************************************************************/
133
134
135
136 /*
137 * We are guaranteed to find all high and low tides as long as their
138 * spacing is greater than Global::eventPrecision.
139 */
140
141
maxMinZeroFn(Timestamp t,unsigned deriv,PredictionValue marklev unusedParameter)142 const PredictionValue Station::maxMinZeroFn (Timestamp t,
143 unsigned deriv,
144 PredictionValue marklev unusedParameter) {
145 return _constituents.tideDerivative (t, deriv+1);
146 }
147
148
markZeroFn(Timestamp t,unsigned deriv,PredictionValue marklev)149 const PredictionValue Station::markZeroFn (Timestamp t,
150 unsigned deriv,
151 PredictionValue marklev) {
152 PredictionValue pv_out = _constituents.tideDerivative (t, deriv);
153 if (deriv == 0)
154 pv_out -= marklev;
155 return pv_out;
156 }
157
158
159 /* findZero (time_t t1, time_t t2, double (*f)(time_t t, int deriv))
160 * Find a zero of the function f, which is bracketed by t1 and t2.
161 * Returns a value which is either an exact zero of f, or slightly
162 * past the zero of f.
163 */
164
165 /*
166 * Here's a root finder based upon a modified Newton-Raphson method.
167 */
168
findZero(Timestamp tl,Timestamp tr,const PredictionValue (Station::* f)(Timestamp t,unsigned deriv,PredictionValue marklev),PredictionValue marklev)169 const Timestamp Station::findZero (Timestamp tl,
170 Timestamp tr,
171 const PredictionValue (Station::*f) (Timestamp t,
172 unsigned deriv,
173 PredictionValue marklev),
174 PredictionValue marklev) {
175 PredictionValue fl = (this->*f)(tl, 0, marklev);
176 PredictionValue fr = (this->*f)(tr, 0, marklev);
177 double scale = 1.0;
178 Interval dt;
179 Timestamp t;
180 PredictionValue fp, ft, f_thresh;
181
182 assert (fl.val() != 0.0 && fr.val() != 0.0);
183 assert (tl < tr);
184 if (fl.val() > 0) {
185 scale = -1.0;
186 fl = -fl;
187 fr = -fr;
188 }
189 assert (fl.val() < 0.0 && fr.val() > 0.0);
190
191 while (tr - tl > Global::eventPrecision) {
192 if (t.isNull())
193 dt = Global::zeroInterval; // Force bisection on first step
194 else if (abs(ft) > f_thresh /* not decreasing fast enough */
195 || (ft.val() > 0.0 ? /* newton step would go outside bracket */
196 (fp <= ft / (t - tl).s()) :
197 (fp <= -ft / (tr - t).s())))
198 dt = Global::zeroInterval; /* Force bisection */
199 else {
200 /* Attempt a newton step */
201 assert (fp.val() != 0.0);
202 // Here I actually do want to round away from zero.
203 dt = Interval(llround(-ft/fp));
204
205 /* Since our goal specifically is to reduce our bracket size
206 as quickly as possible (rather than getting as close to
207 the zero as possible) we should ensure that we don't take
208 steps which are too small. (We'd much rather step over
209 the root than take a series of steps that approach the
210 root rapidly but from only one side.) */
211 if (abs(dt) < Global::eventPrecision)
212 dt = (ft.val() < 0.0 ? Global::eventPrecision
213 : -Global::eventPrecision);
214
215 t += dt;
216 if (t >= tr || t <= tl)
217 dt = Global::zeroInterval; /* Force bisection if outside bracket */
218 f_thresh = abs(ft) / 2.0;
219 }
220
221 if (dt == Global::zeroInterval) {
222 /* Newton step failed, do bisection */
223 t = tl + (tr - tl) / 2;
224 f_thresh = fr > -fl ? fr : -fl;
225 }
226
227 if ((ft = scale * (this->*f)(t,0,marklev)).val() == 0.0)
228 return t; /* Exact zero */
229 else if (ft.val() > 0.0)
230 tr = t, fr = ft;
231 else
232 tl = t, fl = ft;
233
234 fp = scale * (this->*f)(t,1,marklev);
235 }
236
237 return tr;
238 }
239
240
241 /* next_zero(time_t t, double (*f)(), double max_fp, double max_fpp)
242 * Find the next zero of the function f which occurs after time t.
243 * The arguments max_fp and max_fpp give the maximum possible magnitudes
244 * that the first and second derivative of f can achieve.
245 *
246 * Algorithm: Our goal here is to bracket the next zero of f ---
247 * then we can use findZero() to quickly refine the root.
248 * So, we will step forward in time until the sign of f changes,
249 * at which point we know we have bracketed a root.
250 * The trick is to use large steps in our search, making
251 * sure the steps are not so large that we inadvertently
252 * step over more than one root.
253 *
254 * The big trick, is that since the tides (and derivatives of
255 * the tides) are all just harmonic series's, it is easy to place
256 * absolute bounds on their values.
257 */
258
259 // This method is only used in one place and is only used for finding
260 // maxima and minima, so I renamed it to nextMaxMin, got rid of the
261 // max_fp and max_fpp parameters, and installed a more convenient out
262 // parameter.
263
264 // Since by definition the tide cannot change direction between maxima
265 // and minima, there is at most one crossing of a given mark level
266 // between min/max points. Therefore, we already have a bracket of
267 // the mark level to give to findZero, and there is no need for a
268 // function like this to find the next mark crossing.
269
nextMaxMin(Timestamp t,TideEvent & tideEvent_out)270 void Station::nextMaxMin (Timestamp t, TideEvent &tideEvent_out) {
271
272 const Amplitude max_fp (_constituents.tideDerivativeMax(2));
273 const Amplitude max_fpp (_constituents.tideDerivativeMax(3));
274
275 Timestamp t_left, t_right;
276 Interval step, step1, step2;
277 PredictionValue f_left, df_left, f_right, df_right, junk;
278 double scale = 1.0;
279
280 t_left = t;
281
282 /* If we start at a zero, step forward until we're past it. */
283 while ((f_left = maxMinZeroFn(t_left,0,junk)).val() == 0.0)
284 t_left += Global::eventPrecision;
285
286 if (f_left.val() < 0.0) {
287 tideEvent_out.eventType = TideEvent::min;
288 } else {
289 tideEvent_out.eventType = TideEvent::max;
290 scale = -1.0;
291 f_left = -f_left;
292 }
293
294 while (true) {
295
296 /* Minimum time to next zero: */
297 step1 = Interval((interval_rep_t)(abs(f_left) / max_fp));
298
299 /* Minimum time to next turning point: */
300 df_left = scale * maxMinZeroFn(t_left,1,junk);
301 step2 = Interval((interval_rep_t)(abs(df_left) / max_fpp));
302
303 if (df_left.val() < 0.0)
304 /* Derivative is in the wrong direction. */
305 step = step1 + step2;
306 else
307 step = step1 > step2 ? step1 : step2;
308
309 if (step < Global::eventPrecision)
310 step = Global::eventPrecision; /* No ridiculously small steps */
311
312 t_right = t_left + step;
313
314 /*
315 * If we hit upon an exact zero, step right until we're off the
316 * zero. If the sign has changed, we are bracketing a desired
317 * root. If the sign hasn't changed, then the zero was at an
318 * inflection point (i.e. a double-zero to within
319 * Global::eventPrecision) and we want to ignore it.
320 */
321 while ((f_right = scale * maxMinZeroFn(t_right,0,junk)).val() == 0.0)
322 t_right += Global::eventPrecision;
323
324 if (f_right.val() > 0.0) { /* Found a bracket */
325 tideEvent_out.eventTime = findZero (t_left,
326 t_right,
327 &Station::maxMinZeroFn,
328 junk);
329 return;
330 }
331
332 t_left = t_right, f_left = f_right;
333 }
334 }
335
336
findMarkCrossing_Dairiki(Timestamp t1,Timestamp t2,PredictionValue marklev,bool & isRising_out)337 const Timestamp Station::findMarkCrossing_Dairiki (Timestamp t1,
338 Timestamp t2,
339 PredictionValue marklev,
340 bool &isRising_out) {
341 if (t1 > t2)
342 std::swap (t1, t2);
343
344 PredictionValue f1 (markZeroFn(t1,0,marklev));
345 PredictionValue f2 (markZeroFn(t2,0,marklev));
346
347 // Fail gently on rotten brackets. (This used to be an assertion.)
348 if (f1 == f2)
349 return Timestamp(); // return null timestamp
350
351 // We need || instead of && to set isRising_out correctly in the
352 // case where there's a zero exactly at t1 or t2.
353 if (!(isRising_out = (f1.val() < 0.0 || f2.val() > 0.0))) {
354 f1 = -f1;
355 f2 = -f2;
356 }
357
358 // Since f1 != f2, we can't get two zeros, so it doesn't matter which
359 // one we check first.
360 if (f1.val() == 0.0)
361 return t1;
362 else if (f2.val() == 0.0)
363 return t2;
364
365 if (f1.val() < 0.0 && f2.val() > 0.0)
366 return findZero (t1, t2, &Station::markZeroFn, marklev);
367
368 return Timestamp(); // Don't have a bracket, return null timestamp.
369 }
370
371
372 /*************************************************************************/
373
374
findSimpleMarkCrossing(Timestamp t1,Timestamp t2,PredictionValue marklev,bool & isRising_out)375 const Timestamp Station::findSimpleMarkCrossing (Timestamp t1,
376 Timestamp t2,
377 PredictionValue marklev,
378 bool &isRising_out) {
379
380 // marklev must compensate for datum and KnotsSquared. See markZeroFn.
381 // Units should already be comparable to datum.
382 marklev -= _constituents.datum();
383 // Correct knots / knots squared
384 if (_constituents.predictUnits() != marklev.Units())
385 marklev.Units (_constituents.predictUnits());
386
387 return findMarkCrossing_Dairiki (t1, t2, marklev, isRising_out);
388 }
389
390
predictTideEvents(Timestamp startTime,Timestamp endTime,TideEventsOrganizer & organizer,TideEventsFilter filter)391 void Station::predictTideEvents (Timestamp startTime,
392 Timestamp endTime,
393 TideEventsOrganizer &organizer,
394 TideEventsFilter filter) {
395 assert (Global::eventPrecision > Global::zeroInterval);
396 if (startTime >= endTime)
397 return;
398
399 addSimpleTideEvents (startTime, endTime, organizer, filter);
400
401 if (filter == noFilter)
402 addSunMoonEvents (startTime, endTime, organizer);
403 }
404
405
addSimpleTideEvents(Timestamp startTime,Timestamp endTime,TideEventsOrganizer & organizer,TideEventsFilter filter)406 void Station::addSimpleTideEvents (Timestamp startTime,
407 Timestamp endTime,
408 TideEventsOrganizer &organizer,
409 TideEventsFilter filter) {
410 bool isRising;
411 TideEvent te;
412
413 // loopTime is the "internal" timestamp for scanning the reference
414 // station. The timestamps of each event get mangled for sub
415 // stations.
416 Timestamp loopTime (startTime - maximumTimeOffset);
417 Timestamp loopEndTime (endTime - minimumTimeOffset);
418
419 // Patience... range is correctly enforced below.
420 while (loopTime <= loopEndTime) {
421 Timestamp previousLoopTime (loopTime);
422
423 // Get next max or min.
424 nextMaxMin (loopTime, te);
425 loopTime = te.eventTime;
426 finishTideEvent (te);
427 if (te.eventTime >= startTime && te.eventTime < endTime)
428 organizer.add (te);
429
430 // Check for slacks, if applicable. Skip the ones that need
431 // interpolation; those are done in
432 // SubordinateStation::addInterpolatedSubstationMarkCrossingEvents.
433 if (filter != maxMin && isCurrent &&
434 ((te.eventType == TideEvent::max && haveFloodBegins()) ||
435 (te.eventType == TideEvent::min && haveEbbBegins()))) {
436 te.eventTime = findSimpleMarkCrossing (previousLoopTime,
437 loopTime,
438 PredictionValue(predictUnits(), 0.0),
439 isRising);
440 if (!(te.eventTime.isNull())) {
441 te.eventType = (isRising ? TideEvent::slackrise
442 : TideEvent::slackfall);
443 finishTideEvent (te);
444 if (te.eventTime >= startTime && te.eventTime < endTime)
445 organizer.add (te);
446 }
447 }
448
449 // Check for mark, if applicable.
450 if ((!isSubordinateStation()) &&
451 (!markLevel.isNull()) &&
452 (filter == noFilter)) {
453 te.eventTime = findSimpleMarkCrossing (previousLoopTime,
454 loopTime,
455 markLevel,
456 isRising);
457 if (!(te.eventTime.isNull())) {
458 te.eventType = (isRising ? TideEvent::markrise
459 : TideEvent::markfall);
460 finishTideEvent (te);
461 if (te.eventTime >= startTime && te.eventTime < endTime)
462 organizer.add (te);
463 }
464 }
465 }
466 }
467
468
469 // Submethod of predictTideEvents.
addSunMoonEvents(Timestamp startTime,Timestamp endTime,TideEventsOrganizer & organizer)470 void Station::addSunMoonEvents (Timestamp startTime,
471 Timestamp endTime,
472 TideEventsOrganizer &organizer) {
473 TideEvent te;
474
475 const Dstr &em = Global::settings["em"].s;
476
477 if (!(coordinates.isNull())) {
478
479 bool S (em.strchr('S') != -1);
480 bool s (em.strchr('s') != -1);
481 bool M (em.strchr('M') != -1);
482 bool m (em.strchr('m') != -1);
483
484 // Add sunrises and sunsets.
485 if (!(S && s)) {
486 te.eventTime = startTime;
487 Skycal::findNextRiseOrSet (te.eventTime,
488 coordinates,
489 Skycal::solar,
490 te);
491 while (te.eventTime < endTime) {
492 if ((te.eventType == TideEvent::sunrise && !S) ||
493 (te.eventType == TideEvent::sunset && !s)) {
494 finishTideEvent (te);
495 organizer.add (te);
496 }
497 Skycal::findNextRiseOrSet (te.eventTime,
498 coordinates,
499 Skycal::solar,
500 te);
501 }
502 }
503
504 // Add moonrises and moonsets.
505 if (!(M && m)) {
506 te.eventTime = startTime;
507 Skycal::findNextRiseOrSet (te.eventTime,
508 coordinates,
509 Skycal::lunar,
510 te);
511 while (te.eventTime < endTime) {
512 if ((te.eventType == TideEvent::moonrise && !M) ||
513 (te.eventType == TideEvent::moonset && !m)) {
514 finishTideEvent (te);
515 organizer.add (te);
516 }
517 Skycal::findNextRiseOrSet (te.eventTime,
518 coordinates,
519 Skycal::lunar,
520 te);
521 }
522 }
523 }
524
525 // Add moon phases.
526 if (em.strchr('p') == -1) {
527 te.eventTime = startTime;
528 Skycal::findNextMoonPhase (te.eventTime, te);
529 while (te.eventTime < endTime) {
530 finishTideEvent (te);
531 organizer.add (te);
532 Skycal::findNextMoonPhase (te.eventTime, te);
533 }
534 }
535 }
536
537
538 // Analogous to predictTideEvents for raw readings.
predictRawEvents(Timestamp startTime,Timestamp endTime,TideEventsOrganizer & organizer)539 void Station::predictRawEvents (Timestamp startTime,
540 Timestamp endTime,
541 TideEventsOrganizer &organizer) {
542 assert (step > Global::zeroInterval);
543 assert (startTime <= endTime);
544 TideEvent te;
545 te.eventType = TideEvent::rawreading;
546 for (Timestamp t = startTime; t < endTime; t += step) {
547 te.eventTime = t;
548 finishTideEvent (te);
549 organizer.add (te);
550 }
551 }
552
553
extendRange(TideEventsOrganizer & organizer,Direction direction,Interval howMuch,TideEventsFilter filter)554 void Station::extendRange (TideEventsOrganizer &organizer,
555 Direction direction,
556 Interval howMuch,
557 TideEventsFilter filter) {
558 assert (howMuch > Global::zeroInterval);
559 Timestamp startTime, endTime;
560 if (direction == forward) {
561 TideEventsReverseIterator it = organizer.rbegin();
562 assert (it != organizer.rend());
563 startTime = it->second.eventTime;
564 endTime = startTime + howMuch;
565 startTime -= Global::eventSafetyMargin;
566 } else {
567 TideEventsIterator it = organizer.begin();
568 assert (it != organizer.end());
569 endTime = it->second.eventTime;
570 startTime = endTime - howMuch;
571 endTime += Global::eventSafetyMargin;
572 }
573 predictTideEvents (startTime, endTime, organizer, filter);
574 }
575
576
extendRange(TideEventsOrganizer & organizer,Direction direction,unsigned howMany)577 void Station::extendRange (TideEventsOrganizer &organizer,
578 Direction direction,
579 unsigned howMany) {
580 assert (howMany);
581 assert (step > Global::zeroInterval);
582 Timestamp startTime, endTime;
583 if (direction == forward) {
584 TideEventsReverseIterator it = organizer.rbegin();
585 assert (it != organizer.rend());
586 startTime = it->second.eventTime + step;
587 endTime = startTime + step * howMany;
588 } else {
589 TideEventsIterator it = organizer.begin();
590 assert (it != organizer.end());
591 endTime = it->second.eventTime;
592 startTime = endTime - step * howMany;
593 }
594 predictRawEvents (startTime, endTime, organizer);
595 }
596
597
finishPredictionValue(PredictionValue pv)598 const PredictionValue Station::finishPredictionValue (PredictionValue pv) {
599 if (Units::isHydraulicCurrent (pv.Units()))
600 pv.Units (Units::flatten (pv.Units()));
601 pv += _constituents.datum();
602 return pv;
603 }
604
605
predictTideLevel(Timestamp predictTime)606 const PredictionValue Station::predictTideLevel (Timestamp predictTime) {
607 return finishPredictionValue (_constituents.tideDerivative (predictTime, 0));
608 }
609
610
611 #ifdef blendingTest
tideLevelBlendValues(Timestamp predictTime,NullablePredictionValue & firstYear_out,NullablePredictionValue & secondYear_out)612 void Station::tideLevelBlendValues (Timestamp predictTime,
613 NullablePredictionValue &firstYear_out,
614 NullablePredictionValue &secondYear_out) {
615 assert (!isSubordinateStation());
616 _constituents.tideDerivativeBlendValues (predictTime,
617 0,
618 firstYear_out,
619 secondYear_out);
620 if (!firstYear_out.isNull())
621 firstYear_out = finishPredictionValue (firstYear_out);
622 if (!secondYear_out.isNull())
623 secondYear_out = finishPredictionValue (secondYear_out);
624 }
625 #endif
626
627
predictUnits() const628 const Units::PredictionUnits Station::predictUnits () const {
629 return Units::flatten (_constituents.predictUnits());
630 }
631
632
setUnits(Units::PredictionUnits units)633 void Station::setUnits (Units::PredictionUnits units) {
634 if (!isCurrent) {
635 _constituents.setUnits (units);
636 if (!markLevel.isNull())
637 if (markLevel.Units() != units)
638 markLevel.Units (units);
639 }
640 }
641
642
aboutMode(Dstr & text_out,Format::Format form,const Dstr & codeset) const643 void Station::aboutMode (Dstr &text_out,
644 Format::Format form,
645 const Dstr &codeset) const {
646 unsigned maximumNameLength = 0;
647 assert (form == Format::text || form == Format::HTML);
648 if (form == Format::HTML)
649 text_out = "<table>\n";
650 else {
651 text_out = (codeset == "VT100" ? Global::VT100_init : (char*)NULL);
652 MetaFieldVector::const_iterator it = _metadata.begin();
653 while (it != _metadata.end()) {
654 if (it->name.length() > maximumNameLength)
655 maximumNameLength = it->name.length();
656 ++it;
657 }
658 }
659 MetaFieldVector::const_iterator it = _metadata.begin();
660 while (it != _metadata.end()) {
661 if (form == Format::HTML) {
662 text_out += "<tr><td valign=top>";
663 text_out += it->name;
664 text_out += "</td><td valign=top><font face=\"monospace\">";
665 Dstr temp (it->value);
666 temp.repstr ("\n", "<br>\n");
667 text_out += temp;
668 text_out += "</font></td></td>\n";
669 } else {
670 Dstr tmp1 (it->name), tmp2 (it->value), tmp3;
671 // Ignore Global::degreeSign if codeset was overridden.
672 if (codeset == "VT100" && (tmp1 == "Coordinates" ||
673 tmp1 == "Flood direction" ||
674 tmp1 == "Ebb direction"))
675 tmp2.repstr ("�", Global::degreeSign);
676 tmp1.pad (maximumNameLength+2);
677 tmp2.getline (tmp3);
678 tmp1 += tmp3;
679 tmp1 += '\n';
680 while (tmp2.length()) {
681 tmp3 = "";
682 tmp3.pad (maximumNameLength+2);
683 tmp1 += tmp3;
684 tmp2.getline (tmp3);
685 tmp1 += tmp3;
686 tmp1 += '\n';
687 }
688 text_out += tmp1;
689 }
690 ++it;
691 }
692 if (form == Format::HTML)
693 text_out += "</table>\n";
694 Global::finalizeCodeset (text_out, codeset, form);
695 }
696
697
finishTideEvent(TideEvent & te)698 void Station::finishTideEvent (TideEvent &te) {
699 te.isCurrent = isCurrent;
700 te.uncorrectedEventTime.makeNull();
701 te.uncorrectedEventLevel.makeNull();
702 if (te.isSunMoonEvent())
703 te.eventLevel.makeNull();
704 else
705 te.eventLevel = predictTideLevel (te.eventTime);
706 }
707
708
709 // Legal forms are c, h, i, l, or t, but c does nothing.
710 // LaTeX has a catch-22 so that this has to be empty if firstpage.
711 // Possibly the concept of document header should be factored out
712 // (VT100 init codes go there).
713 // FIXME, messy duplication between LaTeX and the others and messy
714 // conditionals.
textBoilerplate(Dstr & text_out,Format::Format form,bool firstpage,double textWidth) const715 void Station::textBoilerplate (Dstr &text_out, Format::Format form,
716 bool firstpage, double textWidth) const {
717 text_out = (char *)NULL;
718 if (form == Format::CSV)
719 return;
720
721 if (form == Format::LaTeX) {
722 if (firstpage)
723 return;
724 Dstr temp (name);
725 temp.LaTeX_mangle();
726 text_out += "{\\Large\\bf \\begin{tabularx}{";
727 text_out += textWidth;
728 text_out += "mm}{Lr}\n";
729 text_out += temp;
730 text_out += " & \\hspace{5mm}";
731 if (coordinates.isNull())
732 text_out += "Coordinates unknown";
733 else {
734 coordinates.print (temp);
735 text_out += temp;
736 }
737 text_out += "\\\\\n\\end{tabularx}}\n\n";
738
739 if (isCurrent) {
740 text_out += "{\\large Flood direction ";
741 if (maxCurrentBearing.isNull())
742 text_out += "unspecified";
743 else {
744 maxCurrentBearing.print (temp);
745 text_out += temp;
746 }
747 text_out += " \\hfill Ebb direction ";
748 if (minCurrentBearing.isNull())
749 text_out += "unspecified";
750 else {
751 minCurrentBearing.print (temp);
752 text_out += temp;
753 }
754 text_out += "}\n\n";
755 }
756
757 if (Global::settings["ou"].c == 'y') {
758 text_out += "Prediction units are ";
759 text_out += Units::longName(predictUnits());
760 MetaFieldVector::const_iterator it = _metadata.begin();
761 while (it != _metadata.end()) {
762 if (it->name == "Datum") {
763 text_out += " relative to ";
764 text_out += it->value;
765 break;
766 }
767 ++it;
768 }
769 text_out += "\n\n";
770 }
771
772 if (!(note.isNull())) {
773 text_out += "Note: ";
774 temp = note;
775 temp.LaTeX_mangle();
776 text_out += temp;
777 text_out += "\n\n";
778 }
779 } else {
780
781 assert (form == Format::HTML ||
782 form == Format::iCalendar ||
783 form == Format::text);
784
785 if (form == Format::iCalendar) {
786
787 // RFC2445 doesn't allow putting very much outside of the VEVENTs.
788 // This makes sense considering that a calendaring tool is only
789 // equipped to display metadata corresponding to specific events.
790
791 // RFC2445 does clearly specify CRLF (CRNL) line discipline.
792
793 text_out += "BEGIN:VCALENDAR\r\n\
794 VERSION:2.0\r\n\
795 PRODID:";
796 // ISO 9070 compliance not mandatory.
797 text_out += XTideVersionString;
798 text_out += "\r\n\
799 CALSCALE:GREGORIAN\r\n\
800 METHOD:PUBLISH\r\n";
801 } else {
802
803 if (form == Format::text && Global::codeset == "VT100" && firstpage)
804 text_out += Global::VT100_init;
805 if (form == Format::HTML) {
806 if (firstpage)
807 text_out += "<h3>";
808 else
809 text_out += "<h3 style=\"page-break-before:always;\">";
810 }
811 text_out += name;
812 if (form == Format::HTML)
813 text_out += "<br>";
814 text_out += '\n';
815 if (coordinates.isNull())
816 text_out += "Coordinates unknown\n";
817 else {
818 Dstr t;
819 coordinates.print (t);
820 if (form == Format::text && Global::needDegrees())
821 t.repstr ("�", Global::degreeSign);
822 text_out += t;
823 text_out += '\n';
824 }
825
826 // When known, append the direction of currents. (The offending
827 // attributes should be null if it's not a current station.)
828 if (!(maxCurrentBearing.isNull())) {
829 if (form == Format::HTML)
830 text_out += "<br>";
831 text_out += "Flood direction ";
832 Dstr tmpbuf;
833 maxCurrentBearing.print (tmpbuf);
834 if (form == Format::text && Global::needDegrees())
835 tmpbuf.repstr ("�", Global::degreeSign);
836 text_out += tmpbuf;
837 text_out += '\n';
838 }
839 if (!(minCurrentBearing.isNull())) {
840 if (form == Format::HTML)
841 text_out += "<br>";
842 text_out += "Ebb direction ";
843 Dstr tmpbuf;
844 minCurrentBearing.print (tmpbuf);
845 if (form == Format::text && Global::needDegrees())
846 tmpbuf.repstr ("�", Global::degreeSign);
847 text_out += tmpbuf;
848 text_out += '\n';
849 }
850
851 if (Global::settings["ou"].c == 'y') {
852 if (form == Format::HTML)
853 text_out += "<br>";
854 text_out += "Prediction units are ";
855 text_out += Units::longName(predictUnits());
856 MetaFieldVector::const_iterator it = _metadata.begin();
857 while (it != _metadata.end()) {
858 if (it->name == "Datum") {
859 text_out += " relative to ";
860 text_out += it->value;
861 break;
862 }
863 ++it;
864 }
865 text_out += '\n';
866 }
867
868 // Similarly for notes
869 if (!(note.isNull())) {
870 if (form == Format::HTML)
871 text_out += "<br>Note: ";
872 else
873 text_out += "Note: ";
874 text_out += note;
875 text_out += '\n';
876 }
877
878 if (form == Format::HTML)
879 text_out += "</h3>";
880 text_out += '\n';
881 Global::finalizeCodeset (text_out, Global::codeset, form);
882 }
883 }
884 }
885
886
887 // iCalendar format output is actually produced by plainMode. From
888 // an engineering perspective this makes perfect sense. But from a
889 // usability perspective, iCalendar output is a calendar and ought
890 // to appear in calendar mode. So calendarMode falls through to
891 // plainMode when i format is chosen.
892
plainMode(Dstr & text_out,Timestamp startTime,Timestamp endTime,Format::Format form)893 void Station::plainMode (Dstr &text_out,
894 Timestamp startTime,
895 Timestamp endTime,
896 Format::Format form) {
897 textBoilerplate (text_out, form, true);
898 TideEventsOrganizer organizer;
899 predictTideEvents (startTime, endTime, organizer);
900 TideEventsIterator it = organizer.begin();
901 while (it != organizer.end()) {
902 Dstr line;
903 it->second.print (line, Mode::plain, form, *this);
904 text_out += line;
905 text_out += '\n';
906 ++it;
907 }
908 if (form == Format::iCalendar)
909 text_out += "END:VCALENDAR\r\n";
910 }
911
912
statsMode(Dstr & text_out,Timestamp startTime,Timestamp endTime)913 void Station::statsMode (Dstr &text_out,
914 Timestamp startTime,
915 Timestamp endTime) {
916
917 textBoilerplate (text_out, Format::text, true);
918 PredictionValue maxl = maxLevelHeuristic();
919 PredictionValue minl = minLevelHeuristic();
920 assert (minl < maxl);
921 PredictionValue meanl = (maxl + minl) / 2.0;
922 Dstr temp;
923 text_out += "Estimated upper bound: ";
924 maxl.print (temp);
925 text_out += temp;
926 text_out += "\nEstimated lower bound: ";
927 minl.print (temp);
928 text_out += temp;
929 text_out += "\nMean, assuming symmetry: ";
930 meanl.print (temp);
931 text_out += temp;
932 text_out += "\n\n";
933
934 bool first (true);
935 TideEventsOrganizer organizer;
936 predictTideEvents (startTime, endTime, organizer, maxMin);
937 Timestamp maxt, mint, lastTidalDay (startTime);
938 PredictionValue sumLevels, sumLLW;
939 NullablePredictionValue LLW;
940 unsigned long numberOfSamples (0), numberOfMLLWSamples (0);
941
942 for (TideEventsIterator it = organizer.begin();
943 it != organizer.end();
944 ++it) {
945 TideEvent &te = it->second;
946 assert (!te.isSunMoonEvent());
947
948 if (!isCurrent) {
949 // MLLW estimation uses the lowest low tide in each tidal day.
950 while (te.eventTime - lastTidalDay >= Global::tidalDay) {
951 if (!LLW.isNull()) {
952 sumLLW += LLW;
953 ++numberOfMLLWSamples;
954 LLW.makeNull();
955 }
956 lastTidalDay += Global::tidalDay;
957 }
958 if (te.eventType == TideEvent::min) {
959 if (LLW.isNull())
960 LLW = te.eventLevel;
961 else if (te.eventLevel < LLW)
962 LLW = te.eventLevel;
963 }
964 }
965
966 sumLevels += te.eventLevel;
967 ++numberOfSamples;
968 if (first || (te.eventLevel < minl)) {
969 mint = te.eventTime;
970 minl = te.eventLevel;
971 }
972 if (first || (te.eventLevel > maxl)) {
973 maxt = te.eventTime;
974 maxl = te.eventLevel;
975 }
976 first = false;
977 }
978 if (!isCurrent)
979 if (endTime - lastTidalDay >= Global::tidalDay && !LLW.isNull()) {
980 sumLLW += LLW;
981 ++numberOfMLLWSamples;
982 }
983
984 text_out += "Searched interval from ";
985 startTime.print (temp, timezone);
986 text_out += temp;
987 text_out += " to ";
988 endTime.print (temp, timezone);
989 text_out += temp;
990 text_out += "\n";
991 if (!first) {
992 text_out += "Maximum was ";
993 maxl.print (temp);
994 text_out += temp;
995 text_out += " at ";
996 maxt.print (temp, timezone);
997 text_out += temp;
998 text_out += '\n';
999
1000 text_out += "Minimum was ";
1001 minl.print (temp);
1002 text_out += temp;
1003 text_out += " at ";
1004 mint.print (temp, timezone);
1005 text_out += temp;
1006 text_out += '\n';
1007
1008 sumLevels /= numberOfSamples;
1009 text_out += "Mean of maxima and minima was ";
1010 sumLevels.print (temp);
1011 text_out += temp;
1012 text_out += '\n';
1013
1014 if (!isCurrent) {
1015 if (numberOfMLLWSamples) {
1016 sumLLW /= numberOfMLLWSamples;
1017 text_out += "Estimated MLLW: ";
1018 sumLLW.print (temp);
1019 text_out += temp;
1020 text_out += '\n';
1021 } else
1022 text_out += "Insufficient data to estimate MLLW.\n";
1023 }
1024 } else {
1025 text_out += "Found no tide events.\n";
1026 }
1027
1028 #ifdef HAVE_SYS_RESOURCE_H
1029 rusage r;
1030 require (getrusage (RUSAGE_SELF, &r) == 0);
1031 text_out += "\nCPU time used: ";
1032 text_out += r.ru_utime.tv_sec + r.ru_utime.tv_usec / 1000000.0;
1033 text_out += " s\n";
1034 #endif
1035 }
1036
1037
calendarMode(Dstr & text_out,Timestamp startTime,Timestamp endTime,Mode::Mode mode,Format::Format form)1038 void Station::calendarMode (Dstr &text_out,
1039 Timestamp startTime,
1040 Timestamp endTime,
1041 Mode::Mode mode,
1042 Format::Format form) {
1043 assert (mode == Mode::calendar || mode == Mode::altCalendar);
1044 assert ((form == Format::CSV && mode == Mode::calendar) ||
1045 form == Format::HTML ||
1046 form == Format::iCalendar ||
1047 form == Format::LaTeX ||
1048 form == Format::text);
1049
1050 if (form == Format::iCalendar)
1051 plainMode (text_out, startTime, endTime, form);
1052 else {
1053 textBoilerplate (text_out, form, true);
1054 std::auto_ptr<Calendar> cal (Calendar::factory (*this,
1055 startTime,
1056 endTime,
1057 mode,
1058 form));
1059 Dstr temp;
1060 cal->print (temp);
1061 text_out += temp;
1062 }
1063 }
1064
1065
rareModes(Dstr & text_out,Timestamp startTime,Timestamp endTime,Mode::Mode mode,Format::Format form)1066 void Station::rareModes (Dstr &text_out,
1067 Timestamp startTime,
1068 Timestamp endTime,
1069 Mode::Mode mode,
1070 Format::Format form) {
1071 assert (form == Format::text || form == Format::CSV);
1072 assert (mode == Mode::raw || mode == Mode::mediumRare);
1073 text_out = (char *)NULL;
1074
1075 TideEventsOrganizer organizer;
1076 predictRawEvents (startTime, endTime, organizer);
1077 TideEventsIterator it = organizer.begin();
1078 while (it != organizer.end()) {
1079 Dstr line;
1080 it->second.print (line, mode, form, *this);
1081 text_out += line;
1082 text_out += '\n';
1083 ++it;
1084 }
1085 }
1086
1087
bannerMode(Dstr & text_out,Timestamp startTime,Timestamp endTime)1088 void Station::bannerMode (Dstr &text_out,
1089 Timestamp startTime,
1090 Timestamp endTime) {
1091 textBoilerplate (text_out, Format::text, true);
1092 std::auto_ptr<Banner> banner (Banner::factory (*this,
1093 Global::settings["tw"].u,
1094 startTime,
1095 endTime));
1096 Dstr temp;
1097 banner->drawTides (this, startTime);
1098 banner->print (temp);
1099 text_out += temp;
1100 }
1101
1102
graphMode(Dstr & text_out,Timestamp startTime,Format::Format form)1103 void Station::graphMode (Dstr &text_out,
1104 Timestamp startTime,
1105 Format::Format form) {
1106 switch (form) {
1107 // print method is not in Graph / don't need yet another class
1108 case Format::text:
1109 {
1110 TTYGraph g (Global::settings["tw"].u, Global::settings["th"].u);
1111 g.drawTides (this, startTime);
1112 g.print (text_out);
1113 }
1114 break;
1115 case Format::SVG:
1116 {
1117 SVGGraph g (Global::settings["gw"].u, Global::settings["gh"].u);
1118 g.drawTides (this, startTime);
1119 g.print (text_out);
1120 }
1121 break;
1122 default:
1123 assert (false);
1124 }
1125 }
1126
1127
clockMode(Dstr & text_out,Format::Format form)1128 void Station::clockMode (Dstr &text_out, Format::Format form) {
1129 switch (form) {
1130 // print method is not in Graph / don't need yet another class
1131 case Format::text:
1132 {
1133 TTYGraph g (Global::settings["tw"].u, Global::settings["th"].u,
1134 Graph::clock);
1135 g.drawTides (this, (time_t)time(NULL));
1136 g.print (text_out);
1137 }
1138 break;
1139 case Format::SVG:
1140 {
1141 SVGGraph g (Global::settings["cw"].u, Global::settings["gh"].u,
1142 Graph::clock);
1143 g.drawTides (this, (time_t)time(NULL));
1144 g.print (text_out);
1145 }
1146 break;
1147 default:
1148 assert (false);
1149 }
1150 }
1151
1152
graphModePNG(FILE * fp,Timestamp startTime)1153 void Station::graphModePNG (FILE *fp, Timestamp startTime) {
1154 RGBGraph g (Global::settings["gw"].u, Global::settings["gh"].u);
1155 g.drawTides (this, startTime);
1156 Global::PNGFile = fp;
1157 g.writeAsPNG (Global::writePNGToFile);
1158 }
1159
1160
clockModePNG(FILE * fp)1161 void Station::clockModePNG (FILE *fp) {
1162 RGBGraph g (Global::settings["cw"].u,
1163 Global::settings["gh"].u,
1164 Graph::clock);
1165 g.drawTides (this, (time_t)time(NULL));
1166 Global::PNGFile = fp;
1167 g.writeAsPNG (Global::writePNGToFile);
1168 }
1169
1170
print(Dstr & text_out,Timestamp startTime,Timestamp endTime,Mode::Mode mode,Format::Format form)1171 void Station::print (Dstr &text_out,
1172 Timestamp startTime,
1173 Timestamp endTime,
1174 Mode::Mode mode,
1175 Format::Format form) {
1176 switch (mode) {
1177 case Mode::about:
1178 if (form != Format::text && form != Format::HTML)
1179 Global::formatBarf (mode, form);
1180 aboutMode (text_out, form, Global::codeset); // Timestamps are redundant
1181 break;
1182
1183 case Mode::plain:
1184 if (form != Format::text && form != Format::CSV)
1185 Global::formatBarf (mode, form);
1186 plainMode (text_out, startTime, endTime, form);
1187 break;
1188
1189 case Mode::stats:
1190 if (form != Format::text)
1191 Global::formatBarf (mode, form);
1192 statsMode (text_out, startTime, endTime);
1193 break;
1194
1195 case Mode::calendar:
1196 switch (form) {
1197 case Format::CSV:
1198 case Format::HTML:
1199 case Format::iCalendar:
1200 case Format::LaTeX:
1201 case Format::text:
1202 calendarMode (text_out, startTime, endTime, mode, form);
1203 break;
1204 default:
1205 Global::formatBarf (mode, form);
1206 }
1207 break;
1208
1209 case Mode::altCalendar:
1210 switch (form) {
1211 case Format::HTML:
1212 case Format::iCalendar:
1213 case Format::LaTeX:
1214 case Format::text:
1215 calendarMode (text_out, startTime, endTime, mode, form);
1216 break;
1217 default:
1218 Global::formatBarf (mode, form);
1219 }
1220 break;
1221
1222 case Mode::raw:
1223 case Mode::mediumRare:
1224 if (form != Format::text && form != Format::CSV)
1225 Global::formatBarf (mode, form);
1226 rareModes (text_out, startTime, endTime, mode, form);
1227 break;
1228
1229 case Mode::banner:
1230 if (form != Format::text)
1231 Global::formatBarf (mode, form);
1232 bannerMode (text_out, startTime, endTime);
1233 break;
1234
1235 case Mode::graph:
1236 switch (form) {
1237 case Format::PNG:
1238 Global::log ("Can't happen: Station::print called for graph mode, PNG form: use graphModePNG instead\n", LOG_ERR);
1239 assert (false);
1240 case Format::text:
1241 case Format::SVG:
1242 graphMode (text_out, startTime, form);
1243 break;
1244 default:
1245 Global::formatBarf (mode, form);
1246 }
1247 break;
1248
1249 case Mode::clock:
1250 switch (form) {
1251 case Format::PNG:
1252 Global::log ("Can't happen: Station::print called for clock mode, PNG form: use graphModePNG instead\n", LOG_ERR);
1253 assert (false);
1254 case Format::text:
1255 case Format::SVG:
1256 clockMode (text_out, form);
1257 break;
1258 default:
1259 Global::formatBarf (mode, form);
1260 }
1261 break;
1262
1263 default:
1264 {
1265 Dstr details ("Unsupported mode: ");
1266 details += (char)mode;
1267 Global::barf (Error::BAD_MODE, details);
1268 }
1269 }
1270 }
1271
1272 }
1273