42 # Jan 5, 2016: - added &velEarthToInstrument(@), &velInstrumentToBeam(@) |
42 # Jan 5, 2016: - added &velEarthToInstrument(@), &velInstrumentToBeam(@) |
43 # Jan 9, 2016: - added &velEarthToBeam(), &velBeamToEarth() |
43 # Jan 9, 2016: - added &velEarthToBeam(), &velBeamToEarth() |
44 # Feb 29, 2016: - debugged & verified velEarthToInstrument(), velInstrumentToBeam() |
44 # Feb 29, 2016: - debugged & verified velEarthToInstrument(), velInstrumentToBeam() |
45 # - added velBeamToEarth() |
45 # - added velBeamToEarth() |
46 # May 19, 2016: - begin implemeting bin interpolation |
46 # May 19, 2016: - begin implemeting bin interpolation |
|
47 # May 25, 2016: - continued |
|
48 # May 26, 2016: - made it work |
47 |
49 |
48 use strict; |
50 use strict; |
49 use POSIX; |
51 use POSIX; |
50 |
52 |
51 my($PI) = 3.14159265358979; |
53 my($PI) = 3.14159265358979; |
52 |
54 |
53 sub rad(@) { return $_[0]/180 * $PI; } |
55 sub rad(@) { return $_[0]/180 * $PI; } |
54 sub deg(@) { return $_[0]/$PI * 180; } |
56 sub deg(@) { return $_[0]/$PI * 180; } |
55 |
57 |
56 $RDI_Coords::minValidVels = 3; # 3-beam solutions ok |
58 #---------------------------------------------------------------------- |
57 |
59 # Tweakables |
58 $RDI_Coords::threeBeam_1 = 0; # stats |
60 #---------------------------------------------------------------------- |
|
61 |
|
62 $RDI_Coords::minValidVels = 3; # 3-beam solutions ok (velBeamToInstrument) |
|
63 $RDI_Coords::binMapping = 'linterp'; # 'linterp' or 'none' (earthVels, BPearthVels) |
|
64 |
|
65 #---------------------------------------------------------------------- |
|
66 # beam to earth transformation (from RDI coord trans manual) |
|
67 #---------------------------------------------------------------------- |
|
68 |
|
69 $RDI_Coords::threeBeam_1 = 0; # stats from velBeamToInstrument |
59 $RDI_Coords::threeBeam_2 = 0; |
70 $RDI_Coords::threeBeam_2 = 0; |
60 $RDI_Coords::threeBeam_3 = 0; |
71 $RDI_Coords::threeBeam_3 = 0; |
61 $RDI_Coords::threeBeam_4 = 0; |
72 $RDI_Coords::threeBeam_4 = 0; |
62 $RDI_Coords::fourBeam = 0; |
73 $RDI_Coords::fourBeam = 0; |
63 |
74 |
66 { # STATIC SCOPE |
77 { # STATIC SCOPE |
67 my(@B2I); |
78 my(@B2I); |
68 |
79 |
69 sub velBeamToInstrument(@) |
80 sub velBeamToInstrument(@) |
70 { |
81 { |
71 my($dta,$ens,$v1,$v2,$v3,$v4) = @_; |
82 my($ADCP,$ens,$v1,$v2,$v3,$v4) = @_; |
72 return undef unless (defined($v1) + defined($v2) + |
83 return undef unless (defined($v1) + defined($v2) + |
73 defined($v3) + defined($v4) |
84 defined($v3) + defined($v4) |
74 >= $RDI_Coords::minValidVels); |
85 >= $RDI_Coords::minValidVels); |
75 |
86 |
76 unless (@B2I) { |
87 unless (@B2I) { |
77 my($a) = 1 / (2 * sin(rad($dta->{BEAM_ANGLE}))); |
88 my($a) = 1 / (2 * sin(rad($ADCP->{BEAM_ANGLE}))); |
78 my($b) = 1 / (4 * cos(rad($dta->{BEAM_ANGLE}))); |
89 my($b) = 1 / (4 * cos(rad($ADCP->{BEAM_ANGLE}))); |
79 my($c) = $dta->{CONVEX_BEAM_PATTERN} ? 1 : -1; |
90 my($c) = $ADCP->{CONVEX_BEAM_PATTERN} ? 1 : -1; |
80 my($d) = $a / sqrt(2); |
91 my($d) = $a / sqrt(2); |
81 @B2I = ([$c*$a, -$c*$a, 0, 0 ], |
92 @B2I = ([$c*$a, -$c*$a, 0, 0 ], |
82 [0, 0, -$c*$a, $c*$a], |
93 [0, 0, -$c*$a, $c*$a], |
83 [$b, $b, $b, $b ], |
94 [$b, $b, $b, $b ], |
84 [$d, $d, -$d, -$d ]); |
95 [$d, $d, -$d, -$d ]); |
115 { # STATIC SCOPE |
126 { # STATIC SCOPE |
116 my($hdg,$pitch,$roll,@I2E); |
127 my($hdg,$pitch,$roll,@I2E); |
117 |
128 |
118 sub velInstrumentToEarth(@) |
129 sub velInstrumentToEarth(@) |
119 { |
130 { |
120 my($dta,$ens,$v1,$v2,$v3,$v4) = @_; |
131 my($ADCP,$ens,$v1,$v2,$v3,$v4) = @_; |
121 return undef unless (defined($v1) && defined($v2) && |
132 return undef unless (defined($v1) && defined($v2) && |
122 defined($v3) && defined($v4) && |
133 defined($v3) && defined($v4) && |
123 defined($dta->{ENSEMBLE}[$ens]->{PITCH}) && |
134 defined($ADCP->{ENSEMBLE}[$ens]->{PITCH}) && |
124 defined($dta->{ENSEMBLE}[$ens]->{ROLL})); |
135 defined($ADCP->{ENSEMBLE}[$ens]->{ROLL})); |
125 |
136 |
126 unless (@I2E && |
137 unless (@I2E && |
127 $pitch == $dta->{ENSEMBLE}[$ens]->{PITCH} && |
138 $pitch == $ADCP->{ENSEMBLE}[$ens]->{PITCH} && |
128 $roll == $dta->{ENSEMBLE}[$ens]->{ROLL}) { |
139 $roll == $ADCP->{ENSEMBLE}[$ens]->{ROLL}) { |
129 printf(STDERR "$0: warning HEADING_ALIGNMENT == %g ignored\n", |
140 printf(STDERR "$0: warning HEADING_ALIGNMENT == %g ignored\n", |
130 $dta->{HEADING_ALIGNMENT}) |
141 $ADCP->{HEADING_ALIGNMENT}) |
131 if ($dta->{HEADING_ALIGNMENT}); |
142 if ($ADCP->{HEADING_ALIGNMENT}); |
132 $hdg = $dta->{ENSEMBLE}[$ens]->{HEADING} - $dta->{HEADING_BIAS} |
143 $hdg = $ADCP->{ENSEMBLE}[$ens]->{HEADING} - $ADCP->{HEADING_BIAS} |
133 if defined($dta->{ENSEMBLE}[$ens]->{HEADING}); |
144 if defined($ADCP->{ENSEMBLE}[$ens]->{HEADING}); |
134 $pitch = $dta->{ENSEMBLE}[$ens]->{PITCH}; |
145 $pitch = $ADCP->{ENSEMBLE}[$ens]->{PITCH}; |
135 $roll = $dta->{ENSEMBLE}[$ens]->{ROLL}; |
146 $roll = $ADCP->{ENSEMBLE}[$ens]->{ROLL}; |
136 my($rad_gimbal_pitch) = atan(tan(rad($pitch)) * cos(rad($roll))); |
147 my($rad_gimbal_pitch) = atan(tan(rad($pitch)) * cos(rad($roll))); |
137 my($sh,$ch) = (sin(rad($hdg)),cos(rad($hdg))) |
148 my($sh,$ch) = (sin(rad($hdg)),cos(rad($hdg))) |
138 if defined($hdg); |
149 if defined($hdg); |
139 my($sp,$cp) = (sin($rad_gimbal_pitch),cos($rad_gimbal_pitch)); |
150 my($sp,$cp) = (sin($rad_gimbal_pitch),cos($rad_gimbal_pitch)); |
140 my($sr,$cr) = (sin(rad($roll)), cos(rad($roll))); |
151 my($sr,$cr) = (sin(rad($roll)), cos(rad($roll))); |
141 @I2E = $dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} |
152 @I2E = $ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} |
142 ? ( |
153 ? ( |
143 [-$ch*$cr-$sh*$sp*$sr, $sh*$cp,-$ch*$sr+$sh*$sp*$cr], |
154 [-$ch*$cr-$sh*$sp*$sr, $sh*$cp,-$ch*$sr+$sh*$sp*$cr], |
144 [-$ch*$sp*$sr+$sh*$cr, $ch*$cp, $sh*$sr+$ch*$sp*$cr], |
155 [-$ch*$sp*$sr+$sh*$cr, $ch*$cp, $sh*$sr+$ch*$sp*$cr], |
145 [+$cp*$sr, $sp, -$cp*$cr, ], |
156 [+$cp*$sr, $sp, -$cp*$cr, ], |
146 ) : ( |
157 ) : ( |
147 [$ch*$cr+$sh*$sp*$sr, $sh*$cp, $ch*$sr-$sh*$sp*$cr], |
158 [$ch*$cr+$sh*$sp*$sr, $sh*$cp, $ch*$sr-$sh*$sp*$cr], |
148 [$ch*$sp*$sr-$sh*$cr, $ch*$cp,-$sh*$sr-$ch*$sp*$cr], |
159 [$ch*$sp*$sr-$sh*$cr, $ch*$cp,-$sh*$sr-$ch*$sp*$cr], |
149 [-$cp*$sr, $sp, $cp*$cr, ], |
160 [-$cp*$sr, $sp, $cp*$cr, ], |
150 ); |
161 ); |
151 } |
162 } |
152 return defined($dta->{ENSEMBLE}[$ens]->{HEADING}) |
163 return defined($ADCP->{ENSEMBLE}[$ens]->{HEADING}) |
153 ? ($v1*$I2E[0][0]+$v2*$I2E[0][1]+$v3*$I2E[0][2], |
164 ? ($v1*$I2E[0][0]+$v2*$I2E[0][1]+$v3*$I2E[0][2], |
154 $v1*$I2E[1][0]+$v2*$I2E[1][1]+$v3*$I2E[1][2], |
165 $v1*$I2E[1][0]+$v2*$I2E[1][1]+$v3*$I2E[1][2], |
155 $v1*$I2E[2][0]+$v2*$I2E[2][1]+$v3*$I2E[2][2], |
166 $v1*$I2E[2][0]+$v2*$I2E[2][1]+$v3*$I2E[2][2], |
156 $v4) |
167 $v4) |
157 : (undef,undef, |
168 : (undef,undef, |
177 { # STATIC SCOPE |
189 { # STATIC SCOPE |
178 my($hdg,$pitch,$roll,@E2I); |
190 my($hdg,$pitch,$roll,@E2I); |
179 |
191 |
180 sub velEarthToInstrument(@) |
192 sub velEarthToInstrument(@) |
181 { |
193 { |
182 my($dta,$ens,$u,$v,$w,$ev) = @_; |
194 my($ADCP,$ens,$u,$v,$w,$ev) = @_; |
183 |
195 |
184 unless (@E2I && |
196 unless (@E2I && |
185 $pitch == $dta->{ENSEMBLE}[$ens]->{PITCH} && |
197 $pitch == $ADCP->{ENSEMBLE}[$ens]->{PITCH} && |
186 $roll == $dta->{ENSEMBLE}[$ens]->{ROLL}) { |
198 $roll == $ADCP->{ENSEMBLE}[$ens]->{ROLL}) { |
187 $hdg = $dta->{ENSEMBLE}[$ens]->{HEADING} - $dta->{HEADING_BIAS} |
199 $hdg = $ADCP->{ENSEMBLE}[$ens]->{HEADING} - $ADCP->{HEADING_BIAS} |
188 if defined($dta->{ENSEMBLE}[$ens]->{HEADING}); |
200 if defined($ADCP->{ENSEMBLE}[$ens]->{HEADING}); |
189 $pitch = $dta->{ENSEMBLE}[$ens]->{PITCH}; |
201 $pitch = $ADCP->{ENSEMBLE}[$ens]->{PITCH}; |
190 $roll = $dta->{ENSEMBLE}[$ens]->{ROLL}; |
202 $roll = $ADCP->{ENSEMBLE}[$ens]->{ROLL}; |
191 my($rad_gimbal_pitch) = atan(tan(rad($pitch)) * cos(rad($roll))); |
203 my($rad_gimbal_pitch) = atan(tan(rad($pitch)) * cos(rad($roll))); |
192 my($useRoll) = ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) ? $roll+180 : $roll; |
204 my($useRoll) = ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) ? $roll+180 : $roll; |
193 my($sh,$ch) = (sin(rad($hdg)),cos(rad($hdg))) |
205 my($sh,$ch) = (sin(rad($hdg)),cos(rad($hdg))) |
194 if defined($hdg); |
206 if defined($hdg); |
195 my($sp,$cp) = (sin($rad_gimbal_pitch),cos($rad_gimbal_pitch)); |
207 my($sp,$cp) = (sin($rad_gimbal_pitch),cos($rad_gimbal_pitch)); |
196 my($sr,$cr) = (sin(rad($useRoll)), cos(rad($useRoll))); |
208 my($sr,$cr) = (sin(rad($useRoll)), cos(rad($useRoll))); |
197 @E2I = ([$ch*$cr+$sh*$sp*$sr, $ch*$sp*$sr-$sh*$cr, -$cp*$sr], # M^-1 = R^-1 * P^-1 * R^-1 |
209 @E2I = ([$ch*$cr+$sh*$sp*$sr, $ch*$sp*$sr-$sh*$cr, -$cp*$sr], # M^-1 = R^-1 * P^-1 * R^-1 |
198 [$sh*$cp, $ch*$cp, $sp ], |
210 [$sh*$cp, $ch*$cp, $sp ], |
199 [$ch*$sr-$sh*$sp*$cr, -$sh*$sr-$ch*$sp*$cr, $cp*$cr]); |
211 [$ch*$sr-$sh*$sp*$cr, -$sh*$sr-$ch*$sp*$cr, $cp*$cr]); |
200 } |
212 } |
201 |
213 |
202 return defined($dta->{ENSEMBLE}[$ens]->{HEADING}) |
214 return defined($ADCP->{ENSEMBLE}[$ens]->{HEADING}) |
203 ? ($u*$E2I[0][0]+$v*$E2I[0][1]+$w*$E2I[0][2], |
215 ? ($u*$E2I[0][0]+$v*$E2I[0][1]+$w*$E2I[0][2], |
204 $u*$E2I[1][0]+$v*$E2I[1][1]+$w*$E2I[1][2], |
216 $u*$E2I[1][0]+$v*$E2I[1][1]+$w*$E2I[1][2], |
205 $u*$E2I[2][0]+$v*$E2I[2][1]+$w*$E2I[2][2], |
217 $u*$E2I[2][0]+$v*$E2I[2][1]+$w*$E2I[2][2], |
206 $ev) |
218 $ev) |
207 : (undef,undef,undef,undef); |
219 : (undef,undef,undef,undef); |
220 { # STATIC SCOPE |
232 { # STATIC SCOPE |
221 my($a,$b,$c,$d); |
233 my($a,$b,$c,$d); |
222 |
234 |
223 sub velInstrumentToBeam(@) |
235 sub velInstrumentToBeam(@) |
224 { |
236 { |
225 my($dta,$x,$y,$z,$ev) = @_; |
237 my($ADCP,$x,$y,$z,$ev) = @_; |
226 return undef unless (defined($x) + defined($y) + |
238 return undef unless (defined($x) + defined($y) + |
227 defined($z) + defined($ev) == 4); |
239 defined($z) + defined($ev) == 4); |
228 |
240 |
229 unless (defined($a)) { |
241 unless (defined($a)) { |
230 $a = 1 / (2 * sin(rad($dta->{BEAM_ANGLE}))); |
242 $a = 1 / (2 * sin(rad($ADCP->{BEAM_ANGLE}))); |
231 $b = 1 / (4 * cos(rad($dta->{BEAM_ANGLE}))); |
243 $b = 1 / (4 * cos(rad($ADCP->{BEAM_ANGLE}))); |
232 $c = $dta->{CONVEX_BEAM_PATTERN} ? 1 : -1; |
244 $c = $ADCP->{CONVEX_BEAM_PATTERN} ? 1 : -1; |
233 $d = $a / sqrt(2); |
245 $d = $a / sqrt(2); |
234 } |
246 } |
235 |
247 |
236 return ( $x/(2*$a*$c) + $z/(4*$b) + $ev/(4*$d), |
248 return ( $x/(2*$a*$c) + $z/(4*$b) + $ev/(4*$d), |
237 -$x/(2*$a*$c) + $z/(4*$b) + $ev/(4*$d), |
249 -$x/(2*$a*$c) + $z/(4*$b) + $ev/(4*$d), |
245 # velEarthToBeam() combines velEarthToInstrument and velInstrumentToBeam |
257 # velEarthToBeam() combines velEarthToInstrument and velInstrumentToBeam |
246 #---------------------------------------------------------------------- |
258 #---------------------------------------------------------------------- |
247 |
259 |
248 sub velEarthToBeam(@) |
260 sub velEarthToBeam(@) |
249 { |
261 { |
250 my($dta,$ens,$u,$v,$w,$ev) = @_; |
262 my($ADCP,$ens,$u,$v,$w,$ev) = @_; |
251 return velInstrumentToBeam($dta, |
263 return velInstrumentToBeam($ADCP, |
252 velEarthToInstrument($dta,$ens,$u,$v,$w,$ev)); |
264 velEarthToInstrument($ADCP,$ens,$u,$v,$w,$ev)); |
253 } |
265 } |
254 |
266 |
255 #====================================================================== |
267 #====================================================================== |
256 # velBeamToBPEarth(@) calculates the vertical- and horizontal vels |
268 # velBeamToBPEarth(@) calculates the vertical- and horizontal vels |
257 # from the two beam pairs separately. Note that (w1+w2)/2 is |
269 # from the two beam pairs separately. Note that (w1+w2)/2 is |
262 { # STATIC SCOPE |
274 { # STATIC SCOPE |
263 my($TwoCosBAngle,$TwoSinBAngle); |
275 my($TwoCosBAngle,$TwoSinBAngle); |
264 |
276 |
265 sub velBeamToBPEarth(@) |
277 sub velBeamToBPEarth(@) |
266 { |
278 { |
267 my($dta,$ens,$b1,$b2,$b3,$b4) = @_; |
279 my($ADCP,$ens,$b1,$b2,$b3,$b4) = @_; |
268 my($v12,$w12,$v34,$w34); |
280 my($v12,$w12,$v34,$w34); |
269 |
281 |
270 return (undef,undef,undef,undef) |
282 return (undef,undef,undef,undef) |
271 unless (defined($dta->{ENSEMBLE}[$ens]->{PITCH}) && |
283 unless (defined($ADCP->{ENSEMBLE}[$ens]->{PITCH}) && |
272 defined($dta->{ENSEMBLE}[$ens]->{ROLL})); |
284 defined($ADCP->{ENSEMBLE}[$ens]->{ROLL})); |
273 |
285 |
274 unless (defined($TwoCosBAngle)) { |
286 unless (defined($TwoCosBAngle)) { |
275 $TwoCosBAngle = 2 * cos(rad($dta->{BEAM_ANGLE})); |
287 $TwoCosBAngle = 2 * cos(rad($ADCP->{BEAM_ANGLE})); |
276 $TwoSinBAngle = 2 * sin(rad($dta->{BEAM_ANGLE})); |
288 $TwoSinBAngle = 2 * sin(rad($ADCP->{BEAM_ANGLE})); |
277 } |
289 } |
278 my($roll) = rad($dta->{ENSEMBLE}[$ens]->{ROLL}); |
290 my($roll) = rad($ADCP->{ENSEMBLE}[$ens]->{ROLL}); |
279 my($sr) = sin($roll); my($cr) = cos($roll); |
291 my($sr) = sin($roll); my($cr) = cos($roll); |
280 my($pitch) = atan(tan(rad($dta->{ENSEMBLE}[$ens]->{PITCH})) * $cr); # gimbal pitch |
292 my($pitch) = atan(tan(rad($ADCP->{ENSEMBLE}[$ens]->{PITCH})) * $cr); # gimbal pitch |
281 my($sp) = sin($pitch); my($cp) = cos($pitch); |
293 my($sp) = sin($pitch); my($cp) = cos($pitch); |
282 |
294 |
283 # Sign convention: |
295 # Sign convention: |
284 # - refer to Coord manual Fig. 3 |
296 # - refer to Coord manual Fig. 3 |
285 # - v12 is horizontal velocity from beam1 to beam2, i.e. westward for upward-looking ADCP |
297 # - v12 is horizontal velocity from beam1 to beam2, i.e. westward for upward-looking ADCP |
286 # with beam 3 pointing north (heading = 0) |
298 # with beam 3 pointing north (heading = 0) |
287 # - w is +ve upward, regardless of instrument orientation |
299 # - w is +ve upward, regardless of instrument orientation |
288 |
300 |
289 my($v12_ic) = ($b1-$b2)/$TwoSinBAngle; # instrument coords with w vertical up |
301 my($v12_ic) = ($b1-$b2)/$TwoSinBAngle; # instrument coords with w vertical up |
290 my($w12_ic) = ($b1+$b2)/$TwoCosBAngle; |
302 my($w12_ic) = ($b1+$b2)/$TwoCosBAngle; |
291 $w12_ic *= -1 if ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
303 $w12_ic *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
292 my($v34_ic) = ($b3-$b4)/$TwoSinBAngle; |
304 my($v34_ic) = ($b3-$b4)/$TwoSinBAngle; |
293 my($w34_ic) = ($b3+$b4)/$TwoCosBAngle; |
305 my($w34_ic) = ($b3+$b4)/$TwoCosBAngle; |
294 $w34_ic *= -1 if ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
306 $w34_ic *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
295 |
307 |
296 if ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) { # beampair Earth coords |
308 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) { # beampair Earth coords |
297 $w12 = $w12_ic*$cr + $v12_ic*$sr - $v34_ic*$sp; |
309 $w12 = $w12_ic*$cr + $v12_ic*$sr - $v34_ic*$sp; |
298 $v12 = $v12_ic*$cr - $w12_ic*$sr + $w34_ic*$sp; |
310 $v12 = $v12_ic*$cr - $w12_ic*$sr + $w34_ic*$sp; |
299 $w34 = $w34_ic*$cp - $v34_ic*$sp + $v12_ic*$sr; |
311 $w34 = $w34_ic*$cp - $v34_ic*$sp + $v12_ic*$sr; |
300 $v34 = $v34_ic*$cp + $w34_ic*$sp - $w12_ic*$sr; |
312 $v34 = $v34_ic*$cp + $w34_ic*$sp - $w12_ic*$sr; |
301 } else { |
313 } else { |
320 { # STATIC SCOPE |
332 { # STATIC SCOPE |
321 my($TwoCosBAngle,$TwoSinBAngle); |
333 my($TwoCosBAngle,$TwoSinBAngle); |
322 |
334 |
323 sub velBeamToBPInstrument(@) |
335 sub velBeamToBPInstrument(@) |
324 { |
336 { |
325 my($dta,$ens,$b1,$b2,$b3,$b4) = @_; |
337 my($ADCP,$ens,$b1,$b2,$b3,$b4) = @_; |
326 my($v12,$w12,$v34,$w34); |
338 my($v12,$w12,$v34,$w34); |
327 |
339 |
328 return (undef,undef,undef,undef) |
340 return (undef,undef,undef,undef) |
329 unless (defined($dta->{ENSEMBLE}[$ens]->{PITCH}) && |
341 unless (defined($ADCP->{ENSEMBLE}[$ens]->{PITCH}) && |
330 defined($dta->{ENSEMBLE}[$ens]->{ROLL})); |
342 defined($ADCP->{ENSEMBLE}[$ens]->{ROLL})); |
331 |
343 |
332 unless (defined($TwoCosBAngle)) { |
344 unless (defined($TwoCosBAngle)) { |
333 $TwoCosBAngle = 2 * cos(rad($dta->{BEAM_ANGLE})); |
345 $TwoCosBAngle = 2 * cos(rad($ADCP->{BEAM_ANGLE})); |
334 $TwoSinBAngle = 2 * sin(rad($dta->{BEAM_ANGLE})); |
346 $TwoSinBAngle = 2 * sin(rad($ADCP->{BEAM_ANGLE})); |
335 } |
347 } |
336 |
348 |
337 # Sign convention: |
349 # Sign convention: |
338 # - refer to Coord manual Fig. 3 |
350 # - refer to Coord manual Fig. 3 |
339 # - v12 is horizontal velocity from beam1 to beam2 |
351 # - v12 is horizontal velocity from beam1 to beam2 |
340 # - w is +ve upward, regardless of instrument orientation |
352 # - w is +ve upward, regardless of instrument orientation |
341 |
353 |
342 if (defined($b1) && defined($b2)) { |
354 if (defined($b1) && defined($b2)) { |
343 $v12 = ($b1-$b2)/$TwoSinBAngle; |
355 $v12 = ($b1-$b2)/$TwoSinBAngle; |
344 $w12 = ($b1+$b2)/$TwoCosBAngle; |
356 $w12 = ($b1+$b2)/$TwoCosBAngle; |
345 $w12 *= -1 if ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
357 $w12 *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
346 } |
358 } |
347 if (defined($b3) && defined($b4)) { |
359 if (defined($b3) && defined($b4)) { |
348 $v34 = ($b3-$b4)/$TwoSinBAngle; |
360 $v34 = ($b3-$b4)/$TwoSinBAngle; |
349 $w34 = ($b3+$b4)/$TwoCosBAngle; |
361 $w34 = ($b3+$b4)/$TwoCosBAngle; |
350 $w34 *= -1 if ($dta->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
362 $w34 *= -1 if ($ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}); |
351 } |
363 } |
352 |
364 |
353 return ($v12,$w12,$v34,$w34); |
365 return ($v12,$w12,$v34,$w34); |
354 } |
366 } |
355 } |
367 } |
360 # Bias correction for beam-coordinate data is done in velInstrumentToEarth() |
372 # Bias correction for beam-coordinate data is done in velInstrumentToEarth() |
361 #====================================================================== |
373 #====================================================================== |
362 |
374 |
363 sub velApplyHdgBias(@) |
375 sub velApplyHdgBias(@) |
364 { |
376 { |
365 my($dta,$ens,$v1,$v2,$v3,$v4) = @_; |
377 my($ADCP,$ens,$v1,$v2,$v3,$v4) = @_; |
366 return (undef,undef,undef,undef) |
378 return (undef,undef,undef,undef) |
367 unless (defined($v1) && defined($v2) && |
379 unless (defined($v1) && defined($v2) && |
368 defined($dta->{ENSEMBLE}[$ens]->{HEADING})); |
380 defined($ADCP->{ENSEMBLE}[$ens]->{HEADING})); |
369 |
381 |
370 my($sh) = sin(rad(-$dta->{HEADING_BIAS})); |
382 my($sh) = sin(rad(-$ADCP->{HEADING_BIAS})); |
371 my($ch) = cos(rad(-$dta->{HEADING_BIAS})); |
383 my($ch) = cos(rad(-$ADCP->{HEADING_BIAS})); |
372 |
384 |
373 return ( $v1*$ch + $v2*$sh, |
385 return ( $v1*$ch + $v2*$sh, |
374 -$v1*$sh + $v2*$ch, |
386 -$v1*$sh + $v2*$ch, |
375 $v3 , |
387 $v3 , |
376 $v4 ); |
388 $v4 ); |
418 return 'nan' unless defined($RDI_pitch) && defined($RDI_roll); |
430 return 'nan' unless defined($RDI_pitch) && defined($RDI_roll); |
419 my($rad_pitch) = atan(tan(rad($RDI_pitch)) * cos(rad($RDI_roll))); |
431 my($rad_pitch) = atan(tan(rad($RDI_pitch)) * cos(rad($RDI_roll))); |
420 return deg(acos(cos($rad_pitch) * cos(rad($RDI_roll)))); |
432 return deg(acos(cos($rad_pitch) * cos(rad($RDI_roll)))); |
421 } |
433 } |
422 |
434 |
|
435 #---------------------------------------------------------------------- |
|
436 # alongBeamDZ(ADCP_dta,ens,beam) => (dz_to_bin1_center,bin_dz) |
|
437 # - calculate vertical distances: |
|
438 # - between transducer and bin1 |
|
439 # - between adjacent bins |
|
440 # - no soundspeed correction |
|
441 # - for UL (Fig. 3 Coord Manual): |
|
442 # b1 = phi + roll b2 = phi - roll |
|
443 # b3 = phi - pitch b4 = phi + pitch |
|
444 # - for DL: |
|
445 # b1 = phi + roll b2 = phi - roll |
|
446 # b3 = phi + pitch b4 = phi - pitch |
|
447 #---------------------------------------------------------------------- |
|
448 |
|
449 sub alongBeamDZ($$$) |
|
450 { |
|
451 my($ADCP,$ens,$beam) = @_; |
|
452 |
|
453 my($tilt); # determine tilt of given beam |
|
454 my($pitch) = $ADCP->{ENSEMBLE}[$ens]->{PITCH}; |
|
455 my($roll) = $ADCP->{ENSEMBLE}[$ens]->{ROLL}; |
|
456 if ($beam == 0) { # beam 1 |
|
457 $tilt = &angle_from_vertical($pitch,$ADCP->{BEAM_ANGLE}+$roll); |
|
458 } elsif ($beam == 1) { # beam 2 |
|
459 $tilt = &angle_from_vertical($pitch,$ADCP->{BEAM_ANGLE}-$roll); |
|
460 } elsif ($beam == 2) { # beam 3 |
|
461 $tilt = $ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} |
|
462 ? &angle_from_vertical($ADCP->{BEAM_ANGLE}-$pitch,$roll) |
|
463 : &angle_from_vertical($ADCP->{BEAM_ANGLE}+$pitch,$roll); |
|
464 } else { # beam 4 |
|
465 $tilt = $ADCP->{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} |
|
466 ? &angle_from_vertical($ADCP->{BEAM_ANGLE}+$pitch,$roll) |
|
467 : &angle_from_vertical($ADCP->{BEAM_ANGLE}-$pitch,$roll); |
|
468 } |
|
469 return ($ADCP->{DISTANCE_TO_BIN1_CENTER}*cos(rad($tilt)), |
|
470 $ADCP->{BIN_LENGTH}*cos(rad($tilt))); |
|
471 } |
|
472 |
|
473 #---------------------------------------------------------------------- |
|
474 # binterp(ADCP_dta,ens,bin,ADCP_field) => @interpolated_vals |
|
475 # - interpolate beam velocities to nominal bin center |
|
476 # - field can be VELOCITY, ECHO_AMPLITUDE, ... |
|
477 # |
|
478 # earthVels(ADCP_dta,ens,bin) => (u,v,w,err_vel) |
|
479 # BPEarthVels(ADCP_dta,ens,bin) => (v12,w12,v34,w34) |
|
480 # - new interface (V1.7) |
|
481 #---------------------------------------------------------------------- |
|
482 |
|
483 sub binterp1($$$$$) # interpolate along a single beam |
|
484 { |
|
485 my($ADCP,$ens,$target_dz,$ADCP_field,$beam) = @_; |
|
486 |
|
487 my($dz2bin1,$bin_dz) = &alongBeamDZ($ADCP,$ens,$beam); |
|
488 my($floor_bin) = int(($target_dz-$dz2bin1) / $bin_dz); |
|
489 $floor_bin-- if ($floor_bin == $ADCP->{N_BINS}-1); |
|
490 |
|
491 my($y1) = $ADCP->{ENSEMBLE}[$ens]->{$ADCP_field}[$floor_bin][$beam]; |
|
492 my($y2) = $ADCP->{ENSEMBLE}[$ens]->{$ADCP_field}[$floor_bin+1][$beam]; |
|
493 $y2 = $y1 unless defined($y2); |
|
494 $y1 = $y2 unless defined($y1); |
|
495 return undef unless defined($y1); |
|
496 |
|
497 my($dz1) = $dz2bin1 + $floor_bin * $bin_dz; |
|
498 my($dz2) = $dz1 + $bin_dz; |
|
499 my($ifac) = ($target_dz - $dz1) / ($dz2 - $dz1); |
|
500 die("assertion failed\nifac = $ifac (target_dz = $target_dz, dz1 = $dz1, dz2 = $dz2)") |
|
501 unless ($ifac>= -0.5 && $ifac<=2); |
|
502 return $y1 + $ifac*($y2-$y1); |
|
503 } |
|
504 |
|
505 sub binterp($$$$) |
|
506 { |
|
507 my($ADCP,$ens,$target_bin,$ADCP_field) = @_; |
|
508 |
|
509 my($crt) = cos(rad($ADCP->{ENSEMBLE}[$ens]->{TILT})); # calc center depth of target bin |
|
510 my($target_dz) = ($ADCP->{DISTANCE_TO_BIN1_CENTER} + $target_bin*$ADCP->{BIN_LENGTH}) * $crt; |
|
511 |
|
512 return (&binterp1($ADCP,$ens,$target_dz,$ADCP_field,0), # interpolate all four beams |
|
513 &binterp1($ADCP,$ens,$target_dz,$ADCP_field,1), |
|
514 &binterp1($ADCP,$ens,$target_dz,$ADCP_field,2), |
|
515 &binterp1($ADCP,$ens,$target_dz,$ADCP_field,3)); |
|
516 } |
|
517 |
|
518 sub earthVels($$$) |
|
519 { |
|
520 my($ADCP,$ens,$bin) = @_; |
|
521 if ($RDI_Coords::binMapping eq 'linterp') { |
|
522 return velInstrumentToEarth($ADCP,$ens, |
|
523 velBeamToInstrument($ADCP,$ens, |
|
524 binterp($ADCP,$ens,$bin,'VELOCITY'))); |
|
525 } elsif ($RDI_Coords::binMapping eq 'none') { |
|
526 return velInstrumentToEarth($ADCP,$ens, |
|
527 velBeamToInstrument($ADCP,$ens, |
|
528 @{$ADCP->{ENSEMBLE}[$ens]->{VELOCITY}[$bin]})); |
|
529 } else { |
|
530 die("earthVels(): unknown bin mapping '$RDI_Coords::binMapping '\n"); |
|
531 } |
|
532 } |
|
533 |
|
534 sub BPEarthVels($$$) |
|
535 { |
|
536 my($ADCP,$ens,$bin) = @_; |
|
537 if ($RDI_Coords::binMapping eq 'linterp') { |
|
538 return velBeamToBPEarth($ADCP,$ens,binterp($ADCP,$ens,$bin,'VELOCITY')); |
|
539 } elsif ($RDI_Coords::binMapping eq 'none') { |
|
540 return velBeamToBPEarth($ADCP,$ens,@{$ADCP->{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}); |
|
541 } else { |
|
542 die("BPEarthVels(): unknown bin mapping '$RDI_Coords::binMapping '\n"); |
|
543 } |
|
544 } |
|
545 |
|
546 #---------------------------------------------------------------------- |
|
547 |
423 1; |
548 1; |