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