OpenCPN Partial API docs
Loading...
Searching...
No Matches
TCDS_Binary_Harmonic.cpp
1/***************************************************************************
2 *
3 * Project: OpenCPN
4 *
5 ***************************************************************************
6 * Copyright (C) 2013 by David S. Register *
7 * *
8 * This program is free software; you can redistribute it and/or modify *
9 * it under the terms of the GNU General Public License as published by *
10 * the Free Software Foundation; either version 2 of the License, or *
11 * (at your option) any later version. *
12 * *
13 * This program is distributed in the hope that it will be useful, *
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16 * GNU General Public License for more details. *
17 * *
18 * You should have received a copy of the GNU General Public License *
19 * along with this program; if not, write to the *
20 * Free Software Foundation, Inc., *
21 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. *
22 **************************************************************************/
23
24#include "TCDS_Binary_Harmonic.h"
25#include "tcmgr.h"
26#include <string>
27/* Declarations for zoneinfo compatibility */
28
29/* Most of these entries are loaded from the tzdata.h include file. That
30 * file was generated from tzdata200c. */
31
32static const char *tz_names[][2] = {
33#include "tzdata.h"
34
35 /* Terminator */
36 {NULL, NULL}};
37
38/* Timelib Time services.
39 * Original XTide source code date: 1997-08-15
40 * Last modified 1998-09-07 by Mike Hopper for WXTide32
41 *
42 * Copyright (C) 1997 David Flater.
43 * Also starring: Geoff Kuenning; Rob Miracle; Dean Pentcheff;
44 * Eric Rosen.
45 *
46 * This program is free software; you can redistribute it and/or modify
47 * it under the terms of the GNU General Public License as published by
48 * the Free Software Foundation; either version 2 of the License, or
49 * (at your option) any later version.
50 *
51 * This program is distributed in the hope that it will be useful,
52 * but WITHOUT ANY WARRANTY; without even the implied warranty of
53 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
54 * GNU General Public License for more details.
55 *
56 * You should have received a copy of the GNU General Public License
57 * along with this program; if not, write to the Free Software
58 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
59 *
60 * Changes by Mike Hopper for WXTide32:
61 * Changed do_timestamp to shorten NT's LONG time zone names just the CAPS
62 * Changed _hpux selector to WIN32 to select HP timezone values.
63 * Added a whole set of remote TZ handler routines, all starting with "tz".
64 */
65
66#ifndef __WXMSW__
67typedef unsigned short WORD;
68typedef long LONG;
69typedef char WCHAR;
70
71typedef struct _SYSTEMTIME {
72 WORD wYear;
73 WORD wMonth;
74 WORD wDayOfWeek;
75 WORD wDay;
76 WORD wHour;
77 WORD wMinute;
78 WORD wSecond;
79 WORD wMilliseconds;
81
82typedef struct _TIME_ZONE_INFORMATION {
83 LONG Bias;
84 WCHAR StandardName[32];
85 SYSTEMTIME StandardDate;
86 LONG StandardBias;
87 WCHAR DaylightName[32];
88 SYSTEMTIME DaylightDate;
89 LONG DaylightBias;
91#else
92#include <windows.h>
93#endif
94
95/*-----------------9/24/2002 4:30PM-----------------
96 * An attempt to get Windoz to work with non-US timezones...
97 * --------------------------------------------------*/
98typedef struct {
100 time_t year_beg;
101 time_t year_end;
102 time_t enter_std;
103 time_t enter_dst;
104 int isdst;
106
107tz_info_entry tz_info_local, tz_info_remote, *tz_info = &tz_info_local;
108
109/*-----------------9/24/2002 8:12AM-----------------
110 * Parse time string in the form [-][hh][:mm][:ss] into seconds.
111 * Returns updated string pointer and signed seconds.
112 * --------------------------------------------------*/
113char *tz_time2sec(char *psrc, long *timesec) {
114 int neg;
115 long temp, mpy;
116 *timesec = 0;
117 mpy = 3600;
118 while (*psrc == ' ') psrc++; /* Skip leading blanks */
119 if (*psrc == '+') psrc++; /* Gobble leading + */
120 if (*psrc == '-') {
121 neg = TRUE;
122 psrc++;
123 } else
124 neg = FALSE;
125
126 do {
127 temp = 0;
128 while (isdigit(*psrc)) temp = temp * 10 + (*(psrc++) - '0');
129
130 *timesec = *timesec + temp * mpy;
131
132 if (*psrc == ':') {
133 mpy /= 60;
134 psrc++;
135 }
136 } while (isdigit(*psrc));
137
138 if (neg) *timesec = 0 - *timesec;
139
140 return (psrc);
141}
142
143/*-----------------9/24/2002 8:16AM-----------------
144 * Parse timezone name string.
145 * Returns string at psrc, updated psrc, and chars copied.
146 * --------------------------------------------------*/
147static char *tz_parse_name(char *psrc, char *pdst, int maxlen) {
148 int nReturn;
149
150 nReturn = 0;
151 while (*psrc == ' ') psrc++; /* Skip leading blanks */
152
153 while (isalpha(*psrc) && nReturn < maxlen) {
154 *(pdst++) = *(psrc++);
155 nReturn++;
156 }
157
158 *pdst = 0;
159 return (psrc);
160}
161
162/*-----------------9/24/2002 8:38AM-----------------
163 * Parse tz rule string into SYSTEMTIME structure.
164 * --------------------------------------------------*/
165static char *tz_parse_rule(char *psrc, SYSTEMTIME *st) {
166 int mol[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
167 long temp, mo;
168 if (*psrc == ',') psrc++; /* Gobble leading comma */
169
170 while (*psrc == ' ') psrc++; /* Skip leading blanks */
171
172 st->wYear = 0;
173 st->wMonth = 0;
174 st->wDay = 0;
175 st->wDayOfWeek = 0;
176 st->wHour = 0;
177 st->wMinute = 0;
178 st->wSecond = 0;
179 st->wMilliseconds = 0;
180
181 if (*psrc == 'J') { /* Julian day (1 <= n <= 365) no leap */
182 psrc++; /* Gobble 'J' */
183 temp = 0;
184 while (isdigit(*psrc)) temp = temp * 10 + (*(psrc++) - '0');
185
186 if (temp < 1 || temp > 365) return (0);
187 temp--;
188 for (mo = 0; temp >= mol[mo]; mo++) temp -= mol[mo];
189 st->wMonth = mo + 1;
190 st->wDay = temp + 1;
191 st->wYear = 1;
192 }
193
194 else if (*psrc == 'M') {
195 psrc++; /* Gobble 'M' */
196
197 temp = 0;
198 while (isdigit(*psrc)) temp = temp * 10 + (*(psrc++) - '0'); /* Get month */
199 if (temp < 1 || temp > 12 || *psrc != '.') return (0);
200 st->wMonth = (unsigned short)temp;
201
202 psrc++; /* Gobble '.' */
203 temp = 0;
204 while (isdigit(*psrc))
205 temp = temp * 10 + (*(psrc++) - '0'); /* Get week number */
206 if (temp < 1 || temp > 5 || *psrc != '.') return (0);
207 st->wDay = (unsigned short)temp;
208
209 psrc++; /* Gobble '.' */
210 temp = 0;
211 while (isdigit(*psrc))
212 temp = temp * 10 + (*(psrc++) - '0'); /* Get day of week number */
213 if (temp < 0 || temp > 6) return (0);
214 st->wDayOfWeek = (unsigned short)temp;
215 }
216
217 if (*psrc == '/') { /* time is specified */
218 psrc++; /* Gobble '/' */
219 psrc = tz_time2sec(psrc, &temp);
220 if (temp < 0 || temp >= 86400) return (0);
221 st->wHour = temp / 3600;
222 temp = temp % 3600;
223 st->wMinute = temp / 60;
224 st->wSecond = temp % 60;
225 }
226 return (psrc);
227}
228
229/*-----------------9/24/2002 3:38PM-----------------
230 * Load tz rule into timezone info data block.
231 * --------------------------------------------------*/
232static void tz_load_rule(char *prule, tz_info_entry *tz_info_remote) {
233 prule = tz_parse_name(prule, (char *)tz_info_remote->tzi.StandardName, 30);
234 prule = tz_time2sec(prule, &tz_info_remote->tzi.Bias);
235 tz_info_remote->tzi.Bias /= 60;
236 tz_info_remote->tzi.StandardBias = 0;
237
238 prule = tz_parse_name(prule, (char *)tz_info_remote->tzi.DaylightName, 30);
239 if (*(char *)tz_info_remote->tzi.DaylightName != '\0') {
240 prule = tz_time2sec(prule, &tz_info_remote->tzi.DaylightBias);
241 tz_info_remote->tzi.DaylightBias /= 60;
242 if (tz_info_remote->tzi.DaylightBias == 0)
243 tz_info_remote->tzi.DaylightBias = -60;
244 else
245 tz_info_remote->tzi.DaylightBias -= tz_info_remote->tzi.Bias;
246
247 if (*prule == ',') {
248 prule = tz_parse_rule(prule, &tz_info_remote->tzi.DaylightDate);
249 if (prule && *prule == ',')
250 tz_parse_rule(prule, &tz_info_remote->tzi.StandardDate);
251 else
252 tz_parse_rule((char *)"M10.5.0/02:00:00",
253 &tz_info_remote->tzi.StandardDate);
254 } else { /* Default is US style tz change */
255 tz_parse_rule((char *)"M4.1.0/02:00:00",
256 &tz_info_remote->tzi.DaylightDate);
257 tz_parse_rule((char *)"M10.5.0/02:00:00",
258 &tz_info_remote->tzi.StandardDate);
259 }
260 } else { /* No DST */
261 tz_info_remote->tzi.DaylightDate.wMonth = 0;
262 tz_info_remote->isdst = 0;
263 }
264}
265
266/* Attempt to load up the local time zone of the location. Moof! */
267static void change_time_zone(const char *tz) {
268 // static char env_string[MAXARGLEN+1];
269 int index;
270
271 if (*tz == ':') tz++; /* Gobble lead-in char */
272 /* Find the translation for the timezone string */
273 index = 0;
274 while (1) {
275 if (tz_names[index][0] == NULL) {
276 tz_info = &tz_info_local;
277 /* Not found. */
278 break;
279 }
280 if (!strcmp(tz_names[index][0], (tz))) {
281 char tz[40];
282 strncpy(tz, tz_names[index][1], 39);
283 tz_load_rule(tz, &tz_info_remote);
284 tz_info = &tz_info_remote;
285 /* Force compute next time this data is used */
286 tz_info->year_beg = 0; // Begin date/time is Jan 1, 1970
287 tz_info->year_end = 0; // End date/time is Jan 1, 1970
288 // sprintf (env_string, "TZ=%s", tz_names[index][1]);
289 break;
290 }
291 index++;
292 }
293}
294
295TCDS_Binary_Harmonic::TCDS_Binary_Harmonic() {
296 m_cst_speeds = NULL;
297 m_cst_nodes = NULL;
298 m_cst_epochs = NULL;
299
300 num_IDX = 0;
301 num_nodes = 0;
302 num_csts = 0;
303 num_epochs = 0;
304
305 // Build the units array
306}
307
308TCDS_Binary_Harmonic::~TCDS_Binary_Harmonic() {}
309
310TC_Error_Code TCDS_Binary_Harmonic::LoadData(const wxString &data_file_path) {
311 if (!open_tide_db(data_file_path.mb_str())) return TC_TCD_FILE_CORRUPT;
312
313 // Build the tables of constituent data
314
315 DB_HEADER_PUBLIC hdr = get_tide_db_header();
316
317 source_ident = wxString(hdr.version, wxConvUTF8);
318
319 num_csts = hdr.constituents;
320 if (0 == num_csts) return TC_GENERIC_ERROR;
321
322 num_nodes = hdr.number_of_years;
323 if (0 == num_nodes) return TC_GENERIC_ERROR;
324
325 // Allocate a working buffer
326 m_work_buffer = (double *)malloc(num_csts * sizeof(double));
327
328 // Constituent speeds
329 m_cst_speeds = (double *)malloc(num_csts * sizeof(double));
330
331 for (int a = 0; a < num_csts; a++) {
332 m_cst_speeds[a] = get_speed(a);
333 m_cst_speeds[a] *= M_PI / 648000; /* Convert to radians per second */
334 }
335
336 // Equilibrium tables by year
337 m_first_year = hdr.start_year;
338 num_epochs = hdr.number_of_years;
339
340 m_cst_epochs = (double **)malloc(num_csts * sizeof(double *));
341 for (int i = 0; i < num_csts; i++)
342 m_cst_epochs[i] = (double *)malloc(num_epochs * sizeof(double));
343
344 for (int i = 0; i < num_csts; i++) {
345 for (int year = 0; year < num_epochs; year++) {
346 m_cst_epochs[i][year] = get_equilibrium(i, year);
347 m_cst_epochs[i][year] *= M_PI / 180.0;
348 }
349 }
350
351 // Node factors
352
353 m_cst_nodes = (double **)malloc(num_csts * sizeof(double *));
354 for (int a = 0; a < num_csts; a++)
355 m_cst_nodes[a] = (double *)malloc(num_nodes * sizeof(double));
356
357 for (int a = 0; a < num_csts; a++) {
358 for (int year = 0; year < num_nodes; year++)
359 m_cst_nodes[a][year] = get_node_factor(a, year);
360 }
361
362 // now load and create the index
363
364 TIDE_RECORD *ptiderec = (TIDE_RECORD *)calloc(sizeof(TIDE_RECORD), 1);
365 for (unsigned int i = 0; i < hdr.number_of_records; i++) {
366 read_tide_record(i, ptiderec);
367
368 num_IDX++; // Keep counting entries for harmonic file stuff
369 IDX_entry *pIDX = new IDX_entry;
370 pIDX->source_data_type = SOURCE_TYPE_BINARY_HARMONIC;
371 pIDX->pDataSource = NULL;
372
373 pIDX->Valid15 = 0;
374
375 pIDX->pref_sta_data = NULL; // no reference data yet
376 pIDX->IDX_Useable = 1; // but assume data is OK
377 pIDX->IDX_tzname = NULL;
378
379 pIDX->IDX_lon = ptiderec->header.longitude;
380 pIDX->IDX_lat = ptiderec->header.latitude;
381
382 const char *tz = get_tzfile(ptiderec->header.tzfile);
383 change_time_zone((char *)tz);
384 if (tz_info) pIDX->IDX_time_zone = -tz_info->tzi.Bias;
385
386 strncpy(pIDX->IDX_station_name, ptiderec->header.name, MAXNAMELEN);
387
388 // Extract a "depth" value from name string, if present.
389 // Name string will contain "(depth xx ft)"
390 std::string name(ptiderec->header.name);
391 size_t n = name.find("depth");
392 if (n != std::string::npos) {
393 std::string d = name.substr(n);
394 std::string dp = d.substr(6);
395 size_t nd = dp.find_first_of(' ');
396 std::string sval = dp.substr(0, nd);
397 int depth = std::stoi(sval);
398 pIDX->current_depth = depth;
399 }
400
401 pIDX->IDX_flood_dir = ptiderec->max_direction;
402 pIDX->IDX_ebb_dir = ptiderec->min_direction;
403
404 if (REFERENCE_STATION == ptiderec->header.record_type) {
405 // Establish Station Type
406 wxString caplin(pIDX->IDX_station_name, wxConvUTF8);
407 caplin.MakeUpper();
408 if (caplin.Contains(_T("CURRENT")))
409 pIDX->IDX_type = 'C';
410 else
411 pIDX->IDX_type = 'T';
412
413 int t1 = ptiderec->zone_offset;
414 double zone_offset = (double)(t1 / 100) + ((double)(t1 % 100)) / 60.;
415 // pIDX->IDX_time_zone = t1a;
416
417 pIDX->IDX_ht_time_off = pIDX->IDX_lt_time_off = 0;
418 pIDX->IDX_ht_mpy = pIDX->IDX_lt_mpy = 1.0;
419 pIDX->IDX_ht_off = pIDX->IDX_lt_off = 0.0;
420 pIDX->IDX_ref_dbIndex = ptiderec->header.reference_station; // will be -1
421
422 // build a Station_Data class, and add to member array
423
424 Station_Data *psd = new Station_Data;
425
426 psd->amplitude = (double *)malloc(num_csts * sizeof(double));
427 psd->epoch = (double *)malloc(num_csts * sizeof(double));
428 psd->station_name = (char *)malloc(ONELINER_LENGTH);
429
430 strncpy(psd->station_name, ptiderec->header.name, MAXNAMELEN);
431 psd->station_type = pIDX->IDX_type;
432
433 // Get meridian, which is seconds difference from UTC, not figuring DST,
434 // so that New York is always (-300 * 60)
435 psd->meridian = -(tz_info->tzi.Bias * 60);
436 psd->zone_offset = zone_offset;
437
438 // Get units
439 strncpy(psd->unit, get_level_units(ptiderec->level_units), 40 - 1);
440 psd->unit[40 - 1] = '\0';
441
442 psd->have_BOGUS = (findunit(psd->unit) != -1) &&
443 (known_units[findunit(psd->unit)].type == BOGUS);
444
445 int unit_c;
446 if (psd->have_BOGUS)
447 unit_c = findunit("knots");
448 else
449 unit_c = findunit(psd->unit);
450
451 if (unit_c != -1) {
452 strncpy(psd->units_conv, known_units[unit_c].name,
453 sizeof(psd->units_conv) - 1);
454 strncpy(psd->units_abbrv, known_units[unit_c].abbrv,
455 sizeof(psd->units_abbrv) - 1);
456 } else {
457 strncpy(psd->units_conv, psd->unit, 40 - 1);
458 psd->units_conv[40 - 1] = '\0';
459 strncpy(psd->units_abbrv, psd->unit, 20 - 1);
460 psd->units_abbrv[20 - 1] = '\0';
461 }
462
463 // Get constituents
464 for (int a = 0; a < num_csts; a++) {
465 psd->amplitude[a] = ptiderec->amplitude[a];
466 psd->epoch[a] = ptiderec->epoch[a] * M_PI / 180.;
467 }
468
469 psd->DATUM = ptiderec->datum_offset;
470
471 m_msd_array.Add(psd); // add it to the member array
472 pIDX->pref_sta_data = psd;
473 pIDX->IDX_ref_dbIndex = i;
474 pIDX->have_offsets = 0;
475 } else if (SUBORDINATE_STATION == ptiderec->header.record_type) {
476 // Establish Station Type
477 wxString caplin(pIDX->IDX_station_name, wxConvUTF8);
478 caplin.MakeUpper();
479 if (caplin.Contains(_T("CURRENT")))
480 pIDX->IDX_type = 'c';
481 else
482 pIDX->IDX_type = 't';
483
484 int t1 = ptiderec->max_time_add;
485 double t1a = (double)(t1 / 100) + ((double)(t1 % 100)) / 60.;
486 t1a *= 60; // Minutes
487 pIDX->IDX_ht_time_off = t1a;
488 pIDX->IDX_ht_mpy = ptiderec->max_level_multiply;
489 if (0. == pIDX->IDX_ht_mpy) pIDX->IDX_ht_mpy = 1.0;
490 pIDX->IDX_ht_off = ptiderec->max_level_add;
491
492 t1 = ptiderec->min_time_add;
493 t1a = (double)(t1 / 100) + ((double)(t1 % 100)) / 60.;
494 t1a *= 60; // Minutes
495 pIDX->IDX_lt_time_off = t1a;
496 pIDX->IDX_lt_mpy = ptiderec->min_level_multiply;
497 if (0. == pIDX->IDX_lt_mpy) pIDX->IDX_lt_mpy = 1.0;
498 pIDX->IDX_lt_off = ptiderec->min_level_add;
499
500 pIDX->IDX_ref_dbIndex = ptiderec->header.reference_station;
501 // strncpy(pIDX->IDX_reference_name, ptiderec->header.name,
502 // MAXNAMELEN);
503
504 if (pIDX->IDX_ht_time_off || pIDX->IDX_ht_off != 0.0 ||
505 pIDX->IDX_lt_off != 0.0 || pIDX->IDX_ht_mpy != 1.0 ||
506 pIDX->IDX_lt_mpy != 1.0)
507 pIDX->have_offsets = 1;
508 }
509
510 m_IDX_array.Add(pIDX);
511 }
512
513 // Mark the index entries individually with invariant harmonic constants
514 unsigned int max_index = GetMaxIndex();
515 for (unsigned int i = 0; i < max_index; i++) {
516 IDX_entry *pIDX = GetIndexEntry(i);
517 if (pIDX) {
518 pIDX->num_nodes = num_nodes;
519 pIDX->num_csts = num_csts;
520 pIDX->num_epochs = num_epochs;
521 pIDX->m_cst_speeds = m_cst_speeds;
522 pIDX->m_cst_nodes = m_cst_nodes;
523 pIDX->m_cst_epochs = m_cst_epochs;
524 pIDX->first_year = m_first_year;
525 pIDX->m_work_buffer = m_work_buffer;
526 }
527 }
528 free(ptiderec);
529
530 return TC_NO_ERROR;
531}
532
533IDX_entry *TCDS_Binary_Harmonic::GetIndexEntry(int n_index) {
534 return &m_IDX_array[n_index];
535}
536
537TC_Error_Code TCDS_Binary_Harmonic::LoadHarmonicData(IDX_entry *pIDX) {
538 // Find the indicated Master station
539 if (!strlen(pIDX->IDX_reference_name)) {
540 strncpy(pIDX->IDX_reference_name, get_station(pIDX->IDX_ref_dbIndex),
541 MAXNAMELEN - 1);
542 pIDX->IDX_reference_name[MAXNAMELEN - 1] = '\0';
543
544 // TIDE_RECORD *ptiderec = (TIDE_RECORD *)calloc(sizeof(TIDE_RECORD),
545 // 1); read_tide_record (pIDX->IDX_ref_dbIndex, ptiderec); free(
546 // ptiderec );
547
548 IDX_entry *pIDX_Ref = &m_IDX_array.Item(pIDX->IDX_ref_dbIndex);
549 Station_Data *pRefSta = pIDX_Ref->pref_sta_data;
550 pIDX->pref_sta_data = pRefSta;
551 pIDX->station_tz_offset =
552 -pRefSta->meridian + (pRefSta->zone_offset * 3600);
553 }
554 return TC_NO_ERROR;
555}
Definition: IDX_entry.h:41