forked from dankelley/oce
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sontek_adp.cpp
145 lines (142 loc) · 6 KB
/
sontek_adp.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
/* vim: set expandtab shiftwidth=2 softtabstop=2 tw=70: */
#include <Rcpp.h>
using namespace Rcpp;
// Cross-reference work:
// 1. update ../src/registerDynamicSymbol.c with an item for this
// 2. main code should use the autogenerated wrapper in ../R/RcppExports.R
//
// [[Rcpp::export]]
IntegerVector do_ldc_sontek_adp(RawVector buf, IntegerVector have_ctd, IntegerVector have_gps, IntegerVector have_bottom_track, IntegerVector pcadp, IntegerVector max)
{
/*
#define DEBUG
*/
// ldc = locate data chunk; _sontek_adp = for a SonTek ADV.
// Arguments:
// buf = buffer with data
// Shave_ctd = 1 if have CTD data (must be 0, for this version)
// Shave_gps = 1 if have GPS data (must be 0, for this version)
// Shave_bottom_track = 1 if have bottom-track data (must be 0, for this version)
// Spcadp = 1 if device is a PCADP (which has longer headers)
// Smax = number of profiles to get (set to <0 to get all)
//
// Method:
// The code checks for bytes as follows, and does a checksum on
// results. It has to determine ncell and nbeam first.
// BYTE Contents
// 1 0xA5 (flag 1)
// 2 0x10 (flag 2)
// 3 0x50 (decimal 80, number of bytes in header)
// 4+ See ADPManual_710.pdf, logical page 84 et seq.
//
// REFERENCES
// 1. see ADPManual_710.pdf, logical pages 82-86.
// 2. hydratools20apr06/adp2cdf.m line 1360 re PCADP's extra header.
unsigned char byte1 = 0xA5;
unsigned char byte2 = 0x10;
unsigned char byte3 = 0x50; /* bytes in header (=80) */
#ifdef DEBUG
Rprintf("have_ctd=%d, have_bottom_track=%d, have_gps=%d, max=%d\n",have_ctd[0],have_bottom_track[0],have_gps[0],max[0]);
#endif
if (have_ctd[0] != 0)
::Rf_error("cannot read SonTek ADP files with CTD data");
if (have_bottom_track[0] != 0)
::Rf_error("cannot read SonTek ADP files with bottom-track data");
if (have_gps[0] != 0)
::Rf_error("cannot read SonTek ADP files with GPS data");
int nbuf = buf.size();
#ifdef DEBUG
Rprintf("nbuf=%d\n", nbuf);
#endif
/* Count matches, so we can allocate the right length */
unsigned int matches = 0;
unsigned short int check_sum_start = ((unsigned short)0xa5<<8) | ((unsigned short)0x96); /* manual p96 says 0xA596; assume little-endian */
if (max[0] < 0)
max[0] = 0;
/* scan first profile to determine ncell and nbeam */
int first_look = 1000;
if (first_look > nbuf)
::Rf_error("cannot read Sontek ADP from a buffer with fewer than 1000 bytes");
int i;
int ncell = -1, nbeam = -1;
for (i = 0; i < first_look - 3; i++) { /* note that we don't look to the very end */
//Rprintf(" %d: %x %x %x (%x %x %x)\n", i, buf[i], buf[i+1], buf[i+2], byte1, byte2, byte3);
if (buf[i] == byte1 && buf[i+1] == byte2 && buf[i+2] == byte3) {
nbeam = (int)buf[i + 26];
ncell = ((unsigned short)buf[i+30]) | ((unsigned short)buf[i+31] << 8);
#ifdef DEBUG
Rprintf("tentative first-profile at buf[%d], yielding nbeam=%d and ncell=%d\n",
i, nbeam, ncell);
#endif
if (nbeam < 2 || nbeam > 3)
::Rf_error("number of beams must be 2 or 3, but it is %d", nbeam);
if (ncell < 1)
::Rf_error("number of cells cannot be less than 1, but it is %d", ncell);
break;
}
}
if (nbeam < 0 || ncell < 0)
::Rf_error("cannot determine #beams or #cells, based on first 1000 bytes in buffer");
// The next line envisions more data streams, e.g. ctd.
int chunk_length = 80 + (have_ctd[0]?16:0) + (have_gps[0]?40:0) + (have_bottom_track[0]?18:0) + 4 * ncell * nbeam;
// Next 2 lines acount for extra header in each PCADP profile; see ref 2.
int max_beams = 4;
int pcadp_extra_header_length = 2*(8+max_beams) + 2*max_beams + max_beams;
if (pcadp[0])
chunk_length += pcadp_extra_header_length;
int bad = 0, maxbad = 100;
#ifdef DEBUG
Rprintf("pcadp=%d pcadp_extra_header_length=%d max_beams=%d\n", pcadp[0], pcadp_extra_header_length, max_beams);
Rprintf("bytes: 0x%x 0x%x 0x%x\n", byte1, byte2, byte3);
Rprintf("chunk_length: %d\n", chunk_length);
#endif
for (int i = 0; i < nbuf - 3 - chunk_length; i++) { // FIXME is 3 right, or needed?
if (buf[i] == byte1 && buf[i+1] == byte2 && buf[i+2] == byte3) {
unsigned short int check_sum = check_sum_start; // RHS is fixed
unsigned short int desired_check_sum = ((unsigned short)buf[i+chunk_length]) | ((unsigned short)buf[i+chunk_length+1] << 8);
for (int c = 0; c < chunk_length; c++)
check_sum += (unsigned short int)buf[i + c];
if (check_sum == desired_check_sum) {
matches++;
#ifdef DEBUG
Rprintf("OK at buf[%d]: check_sum=%d (should be %d); check_sum_start=%d\n",
i, check_sum, desired_check_sum, check_sum_start);
#endif
if (max[0] != 0 && matches >= (unsigned int)max[0])
break;
} else {
#ifdef DEBUG
Rprintf("BAD at buf[%d]: check_sum=%d (should be %d); check_sum_start=%d\n",
i, check_sum, desired_check_sum, check_sum_start);
#endif
if (bad++ > maxbad)
::Rf_error("bad=%d exceeds maxbad=%d\n", bad, maxbad);
}
}
}
/* allocate space, then run through whole buffer again, noting the matches */
unsigned int nres = matches;
IntegerVector res(nres>0?nres:1, 1);
if (nres > 0) {
#ifdef DEBUG
Rprintf("getting space for %d matches\n", nres);
#endif
unsigned int ires = 0;
for (int i=0; i<(nbuf-3-chunk_length); i++) { // FIXME is 3 right, or needed?
if (buf[i] == byte1 && buf[i+1] == byte2 && buf[i+2] == byte3) {
unsigned short int check_sum = check_sum_start; // RHS is fixed
unsigned short int desired_check_sum = ((unsigned short)buf[i+chunk_length]) | ((unsigned short)buf[i+chunk_length+1] << 8);
for (int c = 0; c < chunk_length; c++)
check_sum += (unsigned short int)buf[i + c];
if (check_sum == desired_check_sum)
res[ires++] = i + 1; /* the +1 is to get R pointers */
if (ires > nres) /* FIXME: or +1? */
break;
}
}
return(res);
} else {
res[0] = NA_INTEGER;
}
return(res);
}