diff --git a/Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m b/Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m new file mode 100644 index 000000000..f53065f94 --- /dev/null +++ b/Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m @@ -0,0 +1,60 @@ +function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_nortek3beam(roll, pitch, distAlongBeams, beamAngle, number_of_beams) +% function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_nortek3beam(roll, pitch, distAlongBeams, beamAngle, number_of_beams) +% +% Adjusts height above sensor for each beam in adcpBinMapping.m +% Can read 3 beam Nortek instruments. +% Uses the beam angle and the roll/pitch of the instrument. +% No adjustement when tilt is zero, +% nonMappedHeightAboveSensorBeam will be the same than distAlongBeams. +% +% Inputs: +% +% roll - 1D single precision array (Ntime), the internal roll angles in degrees +% pitch - 1D single precision array (Ntime), the internal pitch angles in degrees +% distAlongBeams - 1D array (Nbins), the internal distance along beam for each bin +% beamAngle - scalar, beam angle from the instrument +% number_of_beams - interger, number of beams on the ADCP +% +% Outputs: +% +% nonMappedHeightAboveSensorBeam - 3D single precision array (Ntime, Nbins, Nbeams), +% tilt adjusted height above sensor for +% each of the 3 beams +% +% Example: +% +% roll_ini = [-180:10:180]'; +% pitch_ini = [-180:10:180]'; +% distAlongBeams = [20:20:600]; +% beamAngle = 0.3; +% number_of_beams = 3; +% +% x = NortekADCP.adcp_adjusted_distances_nortek3beam(roll_ini,pitch_ini,distAlongBeams,beamAngle,number_of_beams); +% +% idx0 = find(roll_ini == 0); +% for i=1:number_of_beams; assert(isequal(x(idx0,:,i), distAlongBeams)); assert(~isequal(x(~idx0,:,i), distAlongBeams)); end +% +% + +narginchk(5,5) + +number_of_beams = 3; +nBins = length(distAlongBeams); + +% set to single precision +nonMappedHeightAboveSensorBeam = nan(length(pitch), length(distAlongBeams), number_of_beams, 'single'); + +nonMappedHeightAboveSensorBeam(:,:,1) = (cos(beamAngle - pitch) / cos(beamAngle)) * distAlongBeams; +nonMappedHeightAboveSensorBeam(:,:,1) = repmat(cos(roll), 1, nBins) .* nonMappedHeightAboveSensorBeam(:,:,1); + +beamAngleX = atan(tan(beamAngle) * cos(60 * pi / 180)); % beams 2 and 3 angle projected on the X axis +beamAngleY = atan(tan(beamAngle) * cos(30 * pi / 180)); % beams 2 and 3 angle projected on the Y axis + +nonMappedHeightAboveSensorBeam(:,:,2) = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams; +nonMappedHeightAboveSensorBeam(:,:,2) = repmat(cos(beamAngleY + roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam(:,:,2); + +nonMappedHeightAboveSensorBeam(:,:,3) = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams; +nonMappedHeightAboveSensorBeam(:,:,3) = repmat(cos(beamAngleY - roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam(:,:,3); + +end + diff --git a/Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m b/Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m new file mode 100644 index 000000000..6cf726231 --- /dev/null +++ b/Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m @@ -0,0 +1,75 @@ +function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_rdi4beam(roll, pitch, pitchSign, distAlongBeams, beamAngle, number_of_beams) +% function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_rdi4beam(roll, pitch, pitchSign, distAlongBeams, beamAngle, number_of_beams) +% +% Adjusts height above sensor for each beam in adcpBinMapping.m +% Can read 4 beam RDI instruments. +% Uses the beam angle and the roll/pitch of the instrument. +% No adjustement when tilt is zero, +% nonMappedHeightAboveSensorBeam will be the same than distAlongBeams. +% Deals with up and down facing ADCPs +% +% Inputs: +% +% roll - 1D single precision array (Ntime), the internal roll angles in degrees +% pitch - 1D single precision array (Ntime), the internal pitch angles in degrees +% pitchSign - integer (-1 or 1), pitch sign based on the ADCP's orientation +% distAlongBeams - 1D array (Nbins), the internal distance along beam for each bin +% beamAngle - scalar, beam angle from the instrument +% number_of_beams - interger, number of beams on the ADCP +% +% Outputs: +% +% nonMappedHeightAboveSensorBeam - 3D single precision array (Ntime, Nbins, Nbeams), +% tilt adjusted height above sensor for +% each of the 4 beams +% +% Example: +% +% %testing for up and down facing ADCPs - no change for beam 1 and 2 (roll) +% %changes only for beam 3 and 4 (pitch) +% +% roll_ini = [-180:10:180]'; +% pitch_ini = [-180:10:180]'; +% pitchSign_up = 1; +% pitchSign_down = -1; +% distAlongBeams = [20:20:600]; +% beamAngle = 0.3; +% number_of_beams = 4; +% +% x_up = TeledyneADCP.adcp_adjusted_distances_rdi4beam(roll_ini,pitch_ini,pitchSign_up,distAlongBeams,beamAngle,number_of_beams); +% x_down = TeledyneADCP.adcp_adjusted_distances_rdi4beam(roll_ini,pitch_ini,pitchSign_down,distAlongBeams,beamAngle,number_of_beams); +% +% assert(isequal(x_up(:,:,1), x_down(:,:,1))); +% assert(isequal(x_up(:,:,2), x_down(:,:,2))); +% assert(~isequal(x_up(:,:,3), x_down(:,:,3))); +% assert(~isequal(x_up(:,:,4), x_down(:,:,4))); +% +% idx0 = find(roll_ini == 0); +% for i=1:number_of_beams; assert(isequal(x_up(idx0,:,i), distAlongBeams)); assert(isequal(x_down(idx0,:,i), distAlongBeams)); end +% for i=1:number_of_beams; assert(~isequal(x_up(~idx0,:,i), distAlongBeams)); assert(~isequal(x_down(~idx0,:,i), distAlongBeams)); end +% +% + +narginchk(6,6) + +nmh=1; +if pitchSign < 0 + % adjust nonMappedHeightAboveSensorBeam sign when downfacing + nmh = -1; + % adjust roll when downfacing + roll = roll-pi; +end + +CP = cos(pitchSign * pitch); +CR = cos(roll); + +% set to single precision +nonMappedHeightAboveSensorBeam = nan(length(pitch), length(distAlongBeams), number_of_beams, 'single'); + +nonMappedHeightAboveSensorBeam(:,:,1) = nmh * (CP .* (cos(beamAngle + roll) / cos(beamAngle) .* distAlongBeams)); +nonMappedHeightAboveSensorBeam(:,:,2) = nmh * (CP .* (cos(beamAngle - roll) / cos(beamAngle) .* distAlongBeams)); +nonMappedHeightAboveSensorBeam(:,:,3) = nmh * (CR .* (cos(beamAngle - pitchSign * pitch) / cos(beamAngle) .* distAlongBeams)); +nonMappedHeightAboveSensorBeam(:,:,4) = nmh * (CR .* (cos(beamAngle + pitchSign * pitch) / cos(beamAngle) .* distAlongBeams)); + +end + diff --git a/Preprocessing/adcpBinMappingPP.m b/Preprocessing/adcpBinMappingPP.m index feef80220..6016cae1a 100644 --- a/Preprocessing/adcpBinMappingPP.m +++ b/Preprocessing/adcpBinMappingPP.m @@ -69,9 +69,9 @@ % RDI 4 beams ADCPs: % It is assumed that the beams are in a convex configuration such as beam 1 % and 2 (respectively 3 and 4) are aligned on the pitch (respectively roll) - % axis. When pitch is positive beam 3 is closer to the surface while beam 4 - % gets further away. When roll is positive beam 2 is closer to the surface - % while beam 1 gets further away. + % axis. For an upward looking ADCP, when pitch is positive beam 3 is closer + % to the surface while beam 4 gets further away. When roll is positive beam + % 2 is closer to the surface while beam 1 gets further away. % % Nortek 3 beams ADCPs: % It is assumed that the beams are in a convex configuration such as beam 1 @@ -81,11 +81,16 @@ % while beams 2 and 3 get further away. When roll is positive beam 3 is % closer to the surface while beam 2 gets further away. % + + % Set pitch sign to deal with downward facing ADCPs + pitchSign = 1; + distAlongBeams = sample_data{k}.dimensions{distAlongBeamsIdx}.data'; if all(diff(distAlongBeams)<0) %invert distAlongBeams so we have increasing values % this is required for the interpolation function. distAlongBeams = distAlongBeams*-1; + pitchSign = -1; end pitch = sample_data{k}.variables{pitchIdx}.data * pi / 180; roll = sample_data{k}.variables{rollIdx}.data * pi / 180; @@ -93,35 +98,16 @@ %TODO: the adjusted distances should be exported % as variables to enable further diagnostic plots and debugging. - %TODO: the adjusted distances should be a Nx4 array or Nx3 array. - if isRDI% RDI 4 beams + if isRDI && absic4Idx ~= 0 % RDI 4 beams number_of_beams = 4; - %TODO: block function. - CP = cos(pitch); - % H[TxB] = P[T] x (+-R[T] x B[B]) - nonMappedHeightAboveSensorBeam1 = CP .* (cos(beamAngle + roll) / cos(beamAngle) .* distAlongBeams); - nonMappedHeightAboveSensorBeam2 = CP .* (cos(beamAngle - roll) / cos(beamAngle) .* distAlongBeams); - - % H[TxB] = R[T] x (-+P[T] x B[B]) - CR = cos(roll); - nonMappedHeightAboveSensorBeam3 = CR .* (cos(beamAngle - pitch) / cos(beamAngle) .* distAlongBeams); - nonMappedHeightAboveSensorBeam4 = CR .* (cos(beamAngle + pitch) / cos(beamAngle) .* distAlongBeams); - else - number_of_beams = 3; - nBins = length(distAlongBeams); - %TODO: block function, include tests. - % Nortek 3 beams - nonMappedHeightAboveSensorBeam1 = (cos(beamAngle - pitch) / cos(beamAngle)) * distAlongBeams; - nonMappedHeightAboveSensorBeam1 = repmat(cos(roll), 1, nBins) .* nonMappedHeightAboveSensorBeam1; - beamAngleX = atan(tan(beamAngle) * cos(60 * pi / 180)); % beams 2 and 3 angle projected on the X axis - beamAngleY = atan(tan(beamAngle) * cos(30 * pi / 180)); % beams 2 and 3 angle projected on the Y axis + nonMappedHeightAboveSensorBeam = TeledyneADCP.adcp_adjusted_distances_rdi4beam(roll, pitch, pitchSign, distAlongBeams, beamAngle, number_of_beams); + + elseif isNortek && absic4Idx == 0 % Nortek 3 beams + number_of_beams = 3; - nonMappedHeightAboveSensorBeam2 = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams; - nonMappedHeightAboveSensorBeam2 = repmat(cos(beamAngleY + roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam2; + nonMappedHeightAboveSensorBeam = NortekADCP.adcp_adjusted_distances_nortek3beam(roll, pitch, distAlongBeams, beamAngle, number_of_beams); - nonMappedHeightAboveSensorBeam3 = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams; - nonMappedHeightAboveSensorBeam3 = repmat(cos(beamAngleY - roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam3; end nSamples = length(pitch); @@ -136,13 +122,13 @@ switch n case 1 - dvar = nonMappedHeightAboveSensorBeam1; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); case 2 - dvar = nonMappedHeightAboveSensorBeam2; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); case 3 - dvar = nonMappedHeightAboveSensorBeam3; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); case 4 - dvar = nonMappedHeightAboveSensorBeam4; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); otherwise errormsg('Beam number %s not supported', num2str(n)) end