# VMT_MBBathy

## PURPOSE

Computes the multibeam bathymetry from the four beams of the ADCP

## SYNOPSIS

function [A] = VMT_MBBathy(z,A,savefile,beamAng,magVar,wsedata,saveaux)

## DESCRIPTION

     Computes the multibeam bathymetry from the four beams of the ADCP
     using a script by D.Mueller (USGS). Beam depths are computed for each
     transect prior to any averaging or mapping.

     Added the ability to export timestamps, pitch, roll, and heading
     (2/1/10)

     Added an identifier for each individual transect to the csv output
    (FEL, 6/14/12)

     V3 Adds capability to handle timeseries of WSE, PRJ 8-7-12

     P.R. Jackson, USGS, 8-7-12

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- VMT --- THE VELOCITY MAPPING TOOLBOX ---

## SOURCE CODE

    0001 function [A] = VMT_MBBathy(z,A,savefile,beamAng,magVar,wsedata,saveaux)
    0002 % Computes the multibeam bathymetry from the four beams of the ADCP
    0003 % using a script by D.Mueller (USGS). Beam depths are computed for each
    0004 % transect prior to any averaging or mapping.
    0005 %
    0006 % Added the ability to export timestamps, pitch, roll, and heading
    0007 % (2/1/10)
    0008 %
    0009 % Added an identifier for each individual transect to the csv output
    0010 %(FEL, 6/14/12)
    0011 %
    0012 % V3 Adds capability to handle timeseries of WSE, PRJ 8-7-12
    0013 %
    0014 % P.R. Jackson, USGS, 8-7-12
    0015
    0016 %% Start
    0017 try
    0018     %disp('Computing corrected beam depths')
    0019     if isstruct(wsedata)
    0020         if length(wsedata.elev) == 1
    0021             %disp('WSE is a constant value')
    0022             wsefiletype = 0;
    0023         else
    0024             %disp('WSE is a timeseries')
    0025             wsefiletype = 1;
    0026         end
    0027     else
    0028         %disp('WSE is a constant value')
    0029         wsedata.elev = wsedata;
    0030         wsefiletype = 0;
    0031     end
    0032
    0033     %% Step through each transect in the given set
    0034     %figure(1); clf
    0035     for zi = 1 : z
    0036
    0037         % Work on the WSE if a vector
    0038         %WSE vector must have a value for each ensemble, so interpolate given
    0039         %values to ensemble times
    0040
    0041         if wsefiletype  %only process as vector if loaded file rather than single value
    0042             %Build an ensemble time vector
    0043             enstime = datenum([A(zi).Sup.year+2000 A(zi).Sup.month A(zi).Sup.day...
    0044                 A(zi).Sup.hour A(zi).Sup.minute (A(zi).Sup.second+A(zi).Sup.sec100./100)]);
    0045             %Had to add 2000 to year--will not work for years < 2000
    0046             %Check the times (for debugging)
    0047             if 0
    0048                 obs_start = datestr(wsedata.obstime(1))
    0049                 obs_end = datestr(wsedata.obstime(end))
    0050
    0051                 ens_start = datestr(enstime(1))
    0052                 ens_end = datestr(enstime(end))
    0053             end
    0054
    0055             % Interpolate the WSE values to the ENS times
    0056             wse = interp1(wsedata.obstime,wsedata.elev,enstime);
    0057             % Plot for debugging
    0058             if 0
    0059                 figure(1); hold on
    0060                 plot(enstime,wse,'k-')
    0061                 datetick('x',13)
    0062                 ylabel('WSE, in meters')
    0063             end
    0064         else
    0065             wse = wsedata.elev; %Single value
    0066         end
    0067
    0068         % Compute position and elevation of each beam depth
    0069         [exyz] = depthxyz(A(zi).Nav.depth,A(zi).Sup.draft_cm,...
    0070             A(zi).Sensor.pitch_deg,A(zi).Sensor.roll_deg,....
    0071             A(zi).Sensor.heading_deg,beamAng,...
    0072             'm',A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,wse,A(zi).Sup.ensNo);  %magVar,removed 4-8-10
    0073
    0074         %Build the auxillary data matrix
    0075         if saveaux
    0076             auxmat = [A(zi).Sup.year+2000 A(zi).Sup.month A(zi).Sup.day...
    0077                 A(zi).Sup.hour A(zi).Sup.minute (A(zi).Sup.second+A(zi).Sup.sec100./100) ...
    0078                 A(zi).Sensor.heading_deg A(zi).Sensor.pitch_deg A(zi).Sensor.roll_deg ...
    0079                 repmat(zi,size(A(zi).Sup.year))]; % Added transect #(zi) FLE 6/14/12    %Had to add 2000 to year--will not work for years < 2000
    0080             auxmat2 = [];
    0081             for i = 1:length(A(zi).Sup.ensNo);
    0082                 dum = repmat(auxmat(i,:),4,1);
    0083                 auxmat2 = cat(1,auxmat2,dum);
    0084             end
    0085             clear auxmat dum enstime wse
    0086         end
    0087
    0088         % Store results
    0089         idxmbb = find(~isnan(exyz(:,4))& ~isnan(exyz(:,2)));
    0090         if zi==1
    0091             zmbb=[exyz(idxmbb,1) exyz(idxmbb,2)...
    0092                 exyz(idxmbb,3) exyz(idxmbb,4)];
    0093             if saveaux
    0094                 auxmbb = auxmat2(idxmbb,:);
    0095             end
    0096         else
    0097             zmbb=cat(1,zmbb,[exyz(idxmbb,1)...
    0098                 exyz(idxmbb,2) exyz(idxmbb,3) exyz(idxmbb,4)]);
    0099             if saveaux
    0100                 auxmbb = cat(1,auxmbb,auxmat2(idxmbb,:));
    0101             end
    0102         end
    0103
    0104         A(zi).Comp.exyz = exyz(idxmbb,:);
    0105
    0106
    0107         clear idxmbb exyz;
    0108         %pack;
    0109     end
    0110
    0111
    0112     %% Save the data
    0113
    0114     if 1
    0115         %disp('Exporting multibeam bathymetry')
    0116         %disp([savefile(1:end-4) '_mbxyz.csv'])
    0117         outfile = [savefile(1:end-4) '.csv'];
    0118         if saveaux
    0119             outmat = [zmbb auxmbb];
    0120             ofid   = fopen(outfile, 'wt');
    0121             outcount = fprintf(ofid,'EnsNo,     Easting_WGS84_m,    Northing_WGS84_m,  Elev_m,  Year,  Month,  Day,  Hour,  Minute,  Second,  Heading_deg,  Pitch_deg,  Roll_deg,  Transect\n'); % Modified to output transect # FLE 6/14/12
    0122             outcount = fprintf(ofid,'%6.0f, %14.2f, %14.2f, %8.2f, %4.0f, %2.0f, %2.0f, %2.0f, %2.0f, %2.2f, %3.3f, %3.3f, %3.3f, %3.0f\n',outmat');
    0123             fclose(ofid);
    0124         else
    0125             outmat = zmbb;
    0126             ofid   = fopen(outfile, 'wt');
    0127             outcount = fprintf(ofid,'EnsNo,     Easting_WGS84_m,    Northing_WGS84_m,  Elev_m\n');
    0128             outcount = fprintf(ofid,'%6.0f, %14.2f, %14.2f, %8.2f\n',outmat');
    0129             fclose(ofid);
    0130         end
    0131         %dlmwrite([savefile(1:end-4) '_mbxyz.csv'],outmat,'precision',15);
    0132     end
    0133 catch err
    0134      if isdeployed
    0135         errLogFileName = fullfile(pwd,...
    0136             ['errorLog' datestr(now,'yyyymmddHHMMSS') '.txt']);
    0137         msgbox({['An unexpected error occurred. Error code: ' err.identifier];...
    0138             ['Error details are being written to the following file: '];...
    0139             errLogFileName},...
    0140             'VMT Status: Unexpected Error',...
    0141             'error');
    0142         fid = fopen(errLogFileName,'W');
    0143         fwrite(fid,err.getReport('extended','hyperlinks','off'));
    0144         fclose(fid);
    0145         rethrow(err)
    0146     else
    0147         msgbox(['An unexpected error occurred. Error code: ' err.identifier],...
    0148             'VMT Status: Unexpected Error',...
    0149             'error');
    0150         rethrow(err);
    0151     end
    0152 end