function epDataFile = adcp2epburst(adcpFile, epDataFile, ADCPtype, dlgFile, settings)

% ADCP2EP is used to translate RDI ADCP data into earth coordinates
% and create an epic compatible data file.
% If the data is in Beam coordinates it will be transformed by runbm2g.m into
% Earth coordinates.  This transformation can be run on workhorse and broad
% band data, but the ADCP type must be specified.
%
% Magnetic Declination
%		If a magnetic declination was provided to the insturment prior to
%	deployment or to rdi2cdf in post-processing it will be applied at this
%	time for both Earth and Beam coordinat data
%
%function epDataFile = adcp2ep(adcpFile, epDataFile, ADCPtype, dlgFile, settings)
%Inputs:
%	adcpFile = the ADcp data file in beam coordinates
%		(Note:if running routines in sequence it should be the trimFile.)
%	epDataFile = the new Epic compatable file that will be created
%	ADCPtype = WH or BB; will default to WH if not specified
%		WH = workhorse, BB = broad band
%		note: if BB, do not need a dlgFile
%	dlgFile = the dialog file that was created when the ADCP was "deployed"
%   settings = a structure with the metadata inputs
%
%     settings.epDataFile = '8351wh.nc'; % final output file name
%     settings.experiment = 'Monterey Canyon 2007'; % your metadata
%     settings.project = 'MC2007'; % your metadata
%     settings.descript = 'Downward Facing Workhorse in MSC Sand wave field'; % your metadata, station or site number, for example
%     settings.SciPi = 'JPX'; % your metadata, principal investigator
%     settings.cmnt = 'barnacles were found on the transducers';  % your metadata
%     settings.water_mass = ' '; % an EPIC requirement
%     settings.adcp2ep.magnetic_variation_at_site = gatt.magnetic_variation;
%     % your data will be rotated by this value!
%     settings.adcp2ep.magnetic_variation_applied = settings.adcp2ep.magnetic_variation_at_site;
%     % the dialog file for your ADCP, there is one in the Demo directory.
%     % the toolbox use the instrument elevation and azimuth specific to each ADCP
%     % which can vary slightly as the manufacturing process does not line up the
%     % transducers exactly, and RDI has a correction for this.
%     % if you want to use the default matrix, use 'default' for the file name
%     settings.dlgFile = 'junk.dlg'; % generated by the ADCP using the PS3 command
%     settings.ADCPtype = 'WH'; % workhorse or BB for broadband
%
%   if you want to generate a burst file with this program, do the
%   following
%   settings.interval = interval between start of each burst
%   settings.nsamples = number of samples per burst
%   settings.startindeces = array of start indeces for the bursts,
%       or [1 Inf] will simply use [1:interval:end]

%	Note: If the names of the files are not given, they will be requested.
%
% Output:
%	epDataFile = same as input


%%% START USGS BOILERPLATE -------------%
% Use of this program is described in:
%
% Acoustic Doppler Current Profiler Data Processing System Manual
% Jessica M. Côté, Frances A. Hotchkiss, Marinna Martini, Charles R. Denham
% Revisions by: Andrée L. Ramsey, Stephen Ruane
% U.S. Geological Survey Open File Report 00-458
% Check for later versions of this Open-File, it is a living document.
%
% Program written in Matlab v7.1.0 SP3
% Program updated in Matlab 7.2.0.232 (R2006a)
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty,
% expressed or implied, is made by the USGS or the United States
% Government as to the accuracy and functioning of the program
% and related program material nor shall the fact of distribution
% constitute any such warranty, and no responsibility is assumed
% by the USGS in connection therewith."
%
%%% END USGS BOILERPLATE --------------

% Written by Jessica M. Cote
% for the U.S. Geological Survey
% Coastal and Marine Geology Program
% Woods Hole, MA
% http://woodshole.er.usgs.gov/
% Please report bugs to jcote@usgs.gov
%
% 24-mar-2010 (MM) convert adcp2ep to a version that writes a burst file

% get the current SVN version- the value is automatically obtained in svn
% is the file's svn.keywords which is set to "Revision"
fprintf('MATLAB %s\n',version)
rev_info = 'SVN $Revision: 1801 $';
fprintf('%s %s running\n',mfilename,rev_info)
diagnostic = 0;

if nargin < 1, help(mfilename), adcpFile=''; end
if nargin < 2, epDataFile=''; end
if nargin < 3, ADCPtype = ''; end
if nargin < 4, dlgFile=''; end

if isempty(adcpFile), adcpFile = '*'; end
if isempty(epDataFile), epDataFile = '*'; end
if isempty(ADCPtype), settings.ADCPtype = 'WH'; end
if isempty(dlgFile), dlgFile = '*'; end


% Open ADCP beam file.
if any(adcpFile == '*')
    [theFile, thePath] = uigetfile(adcpFile, 'Select ADCP File:');
    if ~any(theFile), return, end
    if thePath(end) ~= filesep, thePath(end+1) = filesep; end
    adcpFile = [thePath theFile];
end

[path,name]=fileparts(adcpFile);
suggest=[name(1:end-1) '.nc'];
%create ADCP Geographic coordinates File
if any(epDataFile == '*')
    [theFile, thePath] = uiputfile(suggest, 'Save Data in geographic coordinates as:');
    if ~any(theFile), return, end
    if thePath(end) ~= filesep, thePath(end+1) = filesep; end
    epDataFile = [thePath theFile];
end
%
if ~exist(adcpFile,'file'),
    fprintf('%s: input file %s not found - aborting\n',mfilename,adcpFile)
    epDataFile = [];
    return,
end
B = netcdf(adcpFile );
if isempty(B), return, end

% are we generating bursts?
if isfield(settings,'interval'), interval = settings.interval;
else     interval = 1024;
end
if isfield(settings,'nsamples'), nsamples = settings.nsamples;
else     nsamples = interval;
end
if isfield(settings,'goodens'), goodens = settings.goodens;
else     goodens = [1 Inf];
end
% sort out the indeces
[r,c] = size(goodens);
if r>2 || c>2, % we have  a vector of starting indeces
    interval = NaN;
    % calculate how many complete bursts we can do
    nbursts = length(goodens);
else % things are normal
    if isinf(goodens(2)), 
        goodens(2) = length(B{'TIM'}); 
    end
    if goodens(1) <= 0, goodens(1) = 1; end
    if goodens(2) > length(B{'TIM'}), goodens(2) = length(B{'TIM'}); end
    % calculate how many complete bursts we can do
    nbursts = floor((goodens(2)-goodens(1)+1)./interval);
end
fprintf('%d bursts will be generated\n',nbursts)

% now get the rest of the attributes that vary with experiment/site
% fields and def must have the same length and values must line up
fields = {'experiment', 'descript', 'project', 'cmnt',...
    'origin', 'SciPi', 'conventions', 'magnetic_variation_applied','magnetic_variation_at_site'};
def = {'','','','none',...
    '','','EPIC/PMEL','0','0'};
if exist('settings','var')
    % deal out the settings
    for ifield = 1:length(fields),
        if isfield(settings,fields{ifield}),
            eval(sprintf('%s = settings.%s;',fields{ifield},fields{ifield}))
            fprintf('%s = settings.%s;\n',fields{ifield},fields{ifield})
        else eval(sprintf('%s = '''';',fields{ifield}))
        end
    end
else
    titletext   = ['Data Collection Information: ',B.MOORING(:)];
    lineNo  = 1;
    dlgresult  = inputdlg(fields,titletext,lineNo,def,'on');
    for ifield = 1:length(fields),
        if isempty(dlgresult{ifield}), dlgresult{ifield} = ' '; end % MATLAB netCDf hates null strings
        eval(sprintf('%s = dlgresult{ifield};',fields{ifield}))
    end
end

%get some information
    
theFillValue = fillval(B{'vel1'});
nbins = size(B('bin'),1);
nbeams = 4;
transducer_off = B{'D'}.transducer_offset_from_bottom (1);
coord = B.transform(:);
orientation = B.orientation(:);
bin_size = B{'D'}.bin_size(:);
center_first_bin = B{'D'}.center_first_bin(:);

% define the new file
G = define_epFile(epDataFile, B, nbursts, nsamples);
G.DATA_ORIGIN = ncchar(origin);
G.Conventions = ncchar(conventions);
G.EXPERIMENT = ncchar(experiment);
G.PROJECT = ncchar(project);
G.DESCRIPT = ncchar(descript);
G.DATA_CMNT = ncchar(cmnt);
G.SciPi = ncchar(SciPi);
% these are getting into the file
%Need to calculate or derive the following
% TODO - makes ure wehave the right sense of the fill flga here
G.FILL_FLAG=nclong(0);
G.COMPOSITE=nclong(0);
% this is added at the very end by add_vardesc.m
%G.VAR_DESC='u:v:w:Werr:AGC:PGd:Tx';

%get the velocity data and set the new velocity variable
%dependent on being in Beam or Earth coordinates
fprintf('Data is in %s coordinates\n',B.transform(:))
fprintf('Heading alignment applied by ADCP (EH) %f\n',B{'Hdg'}.heading_alignment(:))
fprintf('Heading bias applied by ADCP or wavesmon (EB) %f\n',B{'Hdg'}.heading_bias(:))
if exist('magnetic_variation_at_site','var'),
    if isempty(magnetic_variation_at_site), magnetic_variation_at_site = 0; end;
    if ischar(magnetic_variation_at_site), magnetic_variation_at_site = str2double(magnetic_variation_at_site); end
    fprintf('User has provided declination metadata of %f \n',magnetic_variation_at_site)
else
    fprintf('User has not provided declination metadata\n')
    magnetic_variation_at_site = [];
end
if exist('magnetic_variation_applied','var'),
    if isempty(magnetic_variation_applied), magnetic_variation_applied = 0; end;
    if ischar(magnetic_variation_applied), magnetic_variation_applied = str2double(magnetic_variation_applied); end
    fprintf('User has requested an additional rotation of %f \n',magnetic_variation_applied)
else
    fprintf('User has not provided an additional rotation value\n')
    magnetic_variation_applied = 0;
end
G.magnetic_variation_at_site = ncfloat(magnetic_variation_at_site);
G.magnetic_variation_applied = ncfloat(magnetic_variation_applied);
if magnetic_variation_applied,
    fprintf('velocities are being rotated by %f\n', magnetic_variation_applied)
end
fprintf('Heading is being rotated by %f\n',magnetic_variation_applied);

if strcmpi(coord,'BEAM'), 
    [theElevations, theAzimuths] = setupbm2g(dlgFile, settings); 
end

% we will accumulate this for later calculations
pressmean = ones(nbursts, 1).*NaN;
brangemean = ones(nbursts, 1).*NaN;

%get the velocity data and set the new velocity variable
vel = cell(4,1);
cur = cell(4,1);
for ibeam = 1:4;
    vel{ibeam} = ones(nsamples,nbins).*NaN;   % Input ADCP velocities.
    cur{ibeam} = ones(nsamples,nbins).*NaN; % preallocate cur
end

% test read for velocity...
for ibeam = 1:4;
    vel{ibeam} = B{['vel' int2str(ibeam)]};
    % this is very necessary - weed out bad data before doing rotations
    vel{ibeam} = autonan(vel{ibeam},1);
end
iens = 1;
p = [vel{1}(iens, :); vel{2}(iens, :); vel{3}(iens, :); vel{4}(iens, :)];
% in R2008a and earlier, the line above was fine.  But the new native
% netCDF interface in R2009a onwards, shapes are different.
if size(p,2) == 1, flipflag = 1; else flipflag = 0; end

tic
fprintf('Working on %d bursts\n',nbursts);
%% Start the main processing loop
for iburst = 1:nbursts, % chunking it
    if ~isnan(interval), S = goodens(1)+(iburst-1)*interval;
    else S = goodens(iburst);
    end
    E = S + nsamples-1;

    G{'burst'}(iburst) = iburst;
    G{'profile'}(iburst,:) = B{'Rec'}(S:E);

    % time used to be adjusted here to shift the time stamp to the midle of
    % the actual set of pings (with RDI one has control of the spacing)
    % this will be corrected after the file is written, just read and
    % convert time for now
    % TODO shift time properly
    tj = B{'TIM'}(S:E);
    G{'time'}(iburst,:)= floor(tj);
    G{'time2'}(iburst,:) = (tj-floor(tj)).*(3600*24*1000);

    if ~isempty(B{'brange'}(S:E)),
        brange = B{'brange'}(S:E);
        brangemean(iburst) = gmean(brange);
        G{'brange'}(iburst,:) = brange;
    else
        brangemean(iburst) = NaN;        
    end
    
    % Do NOT Average echo intensity & percent good.
    for ibeam = 1:nbeams,
        G{sprintf('AGC%1d_1202',ibeam)}(iburst,:,:) = ...
            B{sprintf('AGC%1d',ibeam)}(S:E,:);
        G{sprintf('PGd%1d_1203',ibeam)}(iburst,:,:) = ...
            B{sprintf('PGd%1d',ibeam)}(S:E,:);
    end
    
    % while vectors are rotated, actual heading was never corrected
    Hdg = B{'Hdg'}(S:E)+magnetic_variation_applied;
    Hdg( Hdg >= 360 ) = Hdg(Hdg >= 360) - 360;
    Hdg( Hdg < 0 ) = Hdg(Hdg < 0) + 360;
    G{'Hdg_1215'}(iburst,:) = Hdg;
    G{'Ptch_1216'}(iburst,:) = B{'Ptch'}(S:E);
    G{'Roll_1217'}(iburst,:) = B{'Roll'}(S:E);
    
    switch coord
        case 'BEAM'
            isample = 1;
            for iens=S:E, % Ow, must find a way to get rid of this for loop
                if flipflag, % in R2009a onwards
                    p = [vel{1}(iens, :)'; vel{2}(iens, :)'; vel{3}(iens, :)'; vel{4}(iens, :)'];
                else
                    p = [vel{1}(iens, :); vel{2}(iens, :); vel{3}(iens, :); vel{4}(iens, :)];
                end
                %disp(sprintf('size p at %s line 188 = [%d %d]',mfilename,size(p)))
                %using autonan %p(p == theFillValue) = nan;
                
                %convert from beam to geo, note the transposes!
                %disp(sprintf('size p.'' at %s line 188 = [%d %d]',mfilename,size(p.')))
                q = bm2geo(p.', theElevations, theAzimuths, Hdg(isample), B{'Ptch'}(iens), B{'Roll'}(iens), orientation).';
                q = q./10; % mm/s to cm/s - faster to do it here
                
                if diagnostic,
                    subplot(1,2,1)
                    plot(p,1:nbins)
                    set(gca,'xlim',[-1000 1000])
                    xlabel(sprintf('BEAM %s',vel{1}.units(:)))
                    title(sprintf('Burst %d, Ens %d',iburst,B{'Rec'}(iens)))
                    
                    subplot(1,2,2)
                    plot(q,1:nbins)
                    set(gca,'xlim',[-1000 1000])
                    xlabel('EARTH')
                    title(sprintf('%s',datestr(datenum(gregorian(tj(isample))))))
                    pause(0.25)
                end
                
                % using autonan %q(isnan(q)) = theFillValue;
                % this is just slightly faster than a for loop
                cur{1}(isample, :) = q(1, :); 
                cur{2}(isample, :) = q(2, :); 
                cur{3}(isample, :) = q(3, :); 
                cur{4}(isample, :) = q(4, :); 
                isample = isample+1;
            end              
        otherwise
            %u = vel{1}(:);
            %v = vel{2}(:);
            % Change by MM 12-feb-2007 prevent double rotation of EARTH data
            % ignore heading bias settings, as these come from the EB command
            % changes by the ADCP or by wavesmon and have already been applied
            % to the heading and thus will be in the vectors already
            % if user supplied a user_requested_heading_correction
            % heading was rotated in rdi2cdf, but since the beam transformation
            % was performed by the ADCP, the vectors are not rotated.  Do it here.
            % this correction will usually be declination, but as this is written,
            % SHIP or INST coordinates can be adjusted by a heading alignment
            if magnetic_variation_applied,
                % is the sign correct for West = negative correction sense?
                % this is correct - F. Lightsom 9-feb-2007
                theta = -1*magnetic_variation_applied;
                %[ur,vr] = uv_rotate(u,v,theta);
                [ur,vr] = uv_rotate(vel{1}(S:E,:),vel{2}(S:E,:),theta);
                %need to reset theFillValue
                urmax = max(max(ur));  vrmax = max(max(vr));
                ur(ur == urmax) = theFillValue;
                vr(vr == vrmax) = theFillValue;
                %fill in rotated velocity values
                [cur{1},cur{2}] = deal(ur,vr);
            end
            cur{3} = vel{3}(S:E,:);
            cur{4} = vel{4}(S:E,:);
    end
       
    fprintf('On burst %d ens# %d %s in %5.2f min\n',...
        iburst,B{'Rec'}(S),...
        datestr(datenum(gregorian(tj(1)))),toc/60)
    G{'u_1205'}(iburst, :,:) = cur{1}(:,:);
    G{'v_1206'}(iburst, :,:) = cur{2}(:,:);
    G{'w_1204'}(iburst, :,:) = cur{3}(:,:);
    G{'Werr_1201'}(iburst, :,:) = cur{4}(:,:);
    
    G{'Tx_1211'}(iburst,:) = B{'Tx'}(S:E);
    if strcmpi(B.conductivity_sensor(:),'YES'),
        G{'S_40'}(iburst,:) = B{'S'}(S:E);
    end
    G{'SV_80'}(iburst,:) = B{'sv'}(S:E);

    if ~isempty(B{'Pressure'}),
        % for P_1294 need to convert pascals back to deca-pascals
        G{'P_1294'}(iburst,:) = B{'Pressure'}(S:E)./10;
        pressmean(iburst) = gmean(B{'Pressure'}(S:E));
    else
        pressmean(iburst) = NaN;
    end
    if ~isempty(B{'PressVar'}),
        G{'SDP_850'}(iburst,:) = B{'PressVar'}(S:E);
    end
end

%Depth Calculation: (added 04-Apr-2003)
%   1.  Use pressure sensor if available
%   2.  Use values from findsurface.m calculated using Trimbins.m (only for upward looking)
%   3.  Use user input values calculated using Trimbins.m for upward looking or user input values asked for now for downward looking
% this water_depth gets changed by trimbins to MSL
% but it is possible for the user to bypass trimbins - add traps
% if trimbins is bypassed, then water_depth is from rdi2cdf, the user's input
wdepth = B{'D'}.WATER_DEPTH(:);
wdepth_source = B{'D'}.WATER_DEPTH_source(:);
wdepth_datum = B{'D'}.datum(:);

% first, see if trimbins was run - or we really know the source of water_depth
if isempty(wdepth) || (wdepth<=0) || isempty(wdepth_source) || ~isempty(findstr(wdepth_source, 'rdi2cdf')),
    % we need to double check
    fprintf('%s: water depth was taken from %s\n', mfilename, wdepth_source)
    fprintf('water_depth = %f, re-evaluating\n', wdepth)
    if ~isempty(B{'Pressure'}), % pressure data is present
        depth_head = gmean(pressmean)/9806.65;   % depth at transducer head in meters
        wdepth = depth_head+transducer_off;
        wdepth_source = 'water depth = MSL from pressure sensor, by adcp2ep';
        %dnote = 'Depth values were calculated using the ADCP Pressure
        %Sensor, assuming 9806.65 Pascals per meter.';
    elseif ~isempty(B{'brange'}) && strcmp(orientation,'UP'), % pressure is preferable, but we can use ADCP range to surface
        depth_head = gmean(brangemean);
        wdepth = depth_head+transducer_off;
        wdepth_source = 'water depth = MSL from acoustic range, by adcp2ep';
        %dnote = 'Depth values were calculated using acoustic range to surface.';
    else % we have nothing
        % TODO - this dialog box was preventing the command line from
        % updating and it loked like MATLAB was hung.  No fix known as
        % yet.
        %             prompt  = {'Enter the mean sea level value in meters:'};
        %             def     = {'0'};
        %             titletext   = 'User must input the nominal water depth information';
        %             lineNo  = 1;
        %             dlgresult  = inputdlg(prompt,titletext,lineNo,def);
        %             MSL = str2double(dlgresult{1});
        if isempty(wdepth), % if compute in trimbins, keep that value
            wdepth = input('Enter the mean sea level value in meters: ');
            fprintf('User input value for depth: %f meters\n',wdepth)
            wdepth_source = 'water depth = MSL from user input by adcp2ep';
        end
        depth_head = wdepth - transducer_off;
        %dnote = 'Depth values were calculated using user input values';
    end %end if pressure/else
else
    depth_head = wdepth - transducer_off;
end
disp(wdepth_source)
% fix up global variable dependent on the dta
G.nominal_sensor_depth = ncfloat(wdepth - B{'D'}.transducer_offset_from_bottom(1));
G.nominal_sensor_depth_note = 'inst_depth = (water_depth - inst_height); nominal depth below surface, meters';
G.WATER_DEPTH = ncfloat(wdepth); % water depth computed from surace detect or pressure sensor
% these we have already determined
G.WATER_DEPTH_source = wdepth_source;
G.WATER_DEPTH_datum = wdepth_datum;
G.VAR_FILL = ncfloat(theFillValue);

% bindist is range to bin from head
disp(orientation)
fprintf('Center_first_bin = %f\n',center_first_bin);
fprintf('bin_size = %f\n',bin_size);
fprintf('depth_head = %f\n',depth_head);
bindist = center_first_bin:bin_size:(center_first_bin+((nbins-1)*bin_size));
%bindist = (0:nbins-1).*bin_size+center_first_bin;
fprintf('%d bin ranges for %d bins\n',length(bindist), nbins);

% now compute the bin depths
switch orientation
    case 'UP'
        depth = depth_head - bindist;
        %dnote = 'uplooking bin depths = ADCP depth - range to bin';
        dnote = 'uplooking bin depths = water_depth - transducer_offset_from_bottom - bindist';
    case 'DOWN'
        depth = depth_head + bindist;
        dnote = 'downlooking bin depths = water_depth - transducer_offset_from_bottom + bindist';
end

% switch orientation
%     case 'UP'
%         % compute the depths of the bins relative tot he sea surface
%         depth = wdepth - B{'D'}(:); % D contains bin depths relative to the sea bed
%         dnote = 'uplooking bin depths = water_depth - bin distance to sea bed';
%     case 'DOWN'
%         bin1 = depth_head+center_first_bin;
%         binEnd = ((nbins-1)*bin_size)+bin1;
%         depth = (bin1:bin_size:binEnd)';
%         dnote = 'downlooking bin depths = ADCP depth + range to bin';
% end %end switch
disp(dnote)
fprintf('Depths are from %6.2f at the head to %6.2f\n',depth(1),depth(end))
G{'depth'}(1:nbins) = depth;
% add new attributes
G{'depth'}.NOTE = ncchar(dnote);
% remove old attributes
delete(G{'depth'}.NOTE_3);
% update attributes
% duplicating as per CRS request 8/5/09
G{'depth'}.WATER_DEPTH = ncfloat(G.WATER_DEPTH);
G{'depth'}.WATER_DEPTH_source = G.WATER_DEPTH_source(:);
G{'depth'}.WATER_DEPTH_datum = G.WATER_DEPTH_datum(:);

G{'bindist'}(1:nbins) = bindist;
vnames = ncnames(var(G));
for ivar = 1:length(vnames),
    if ~isempty(G{vnames{ivar}}.sensor_depth),
        G{vnames{ivar}}.sensor_depth = ncfloat(wdepth - transducer_off);
    end
end

% % we want to shift the time stamp tothe
% % middle of the set of pings
% %let's deal with the time
% ppens=B.pings_per_ensemble(:);
% tp=B.time_between_ping_groups(:);
% %in order to ping as fast as possible tp may be set to 0,
% %but in reality it is approximately 0.288 sec per ping.
% if isequal(tp,0);
%     tp = 0.288;
% end
% 
% tens=ppens*tp;
% 
% % TODO - this is memory intensive, this could be done better
% %Tim is the beginning of the ensemble, make it the middle
% disp('Converting "TIM" to "time" and "time2"');
% m = size(Mid_Ens_time, 1);
% alltime = (ones(m,2).*NaN);
% for ii=1:m;
%     alltime(ii,:)=ep_time(Mid_Ens_time(ii,:));
%     if (ii<1000 && ~rem(ii,100)) || (ii>1000 && ~rem(ii,1000)),
%         fprintf('%d converted\n',ii),
%     end
% end
% 
% Time=alltime(:,1);
% Time2=alltime(:,2);
% 
% clear alltime
% tsec=gtD(:,4)*3600 + gtD(:,5)*60 + gtD(:,6);
% tmid=tsec+tens/2;
% disp(['TIM was corrected by ' num2str(tens/2) ' sec = half the ensemble time']);
% disp(' ')
% sec=rem(tmid,60);
% hmt=(tmid-sec)/60;
% minu=rem(hmt,60);
% hr=(hmt-minu)/60;
% Mid_Ens_time=[gtD(:,1) gtD(:,2) gtD(:,3) hr minu sec];
% 
% allTIM=B{'TIM'}(:);
% % calculate time step for output file, in seconds (FH 10 May 2000)
% timsecs=B{'TIM'}(:)*24*3600;
% delta=mean(diff(timsecs));
% %clear timsecs
% 
% gtD=gregorian(allTIM);
% start_time = datestr(datenum(gtD(1,1),gtD(1,2),gtD(1,3),gtD(1,4),gtD(1,5),gtD(1,6)),0);
% stop_time = datestr(datenum(gtD(end,1),gtD(end,2),gtD(end,3),gtD(end,4),gtD(end,5),gtD(end,6)),0);
% 
% %calculate dt** this section removed FH 10 May 2000
% %tim=B{'TIM'};
% %gt=gregorian(tim(:));
% %dt=diff(gt);
% %delta=(dt(1,4).*60)+(dt(1,5))+(dt(1,6)./60);
% 
% G.DELTA_T=ncchar(num2str(delta));
% G.start_time = start_time;
% G.stop_time = stop_time;

% add minimums and maximums
% G=add_minmaxvalues(G);
% add the variable description string
add_vardesc(G);

% fix one little bit of sematics
close(B)
close(G)

if exist('bm2gcomment','var'),
    thecomment = sprintf('%s\n Written to an EPIC standard data file by %s %s; \n',...
        bm2gcomment, mfilename, rev_info);
else
    thecomment = sprintf('Written to an EPIC standard data file by %s %s; \n',...
        mfilename, rev_info);
end
history(epDataFile,thecomment);

return

function [theElevations, theAzimuths, bm2gcomment]...
    = setupbm2g(dlgFile, settings)

        % bm2geo can also prompt for the dialog file, but doing it here
        % so history can be updated properly without messing with
        % bm2geo too much
        if any(dlgFile == '*') || isempty(dlgFile),
            [theFile, thePath] = uigetfile({'*.dlg';'*.log';'*.txt';'*.*'},...
                'Select ADCP Deployment Log File:');
            if ~any(theFile),
                dlgFile = 'default';
            else
                if thePath(end) ~= filesep, thePath(end+1) = filesep; end
                dlgFile = [thePath theFile];
            end
        end
        if ~strcmp(dlgFile,'default'),
            bm2gcomment = sprintf('Transformed to earth coordinates by bm2geo, transformation matrix in %s;\n',dlgFile);
        else
            % for someone who just wants to use the default
            dlgFile = 'default';
            % these defaults are set in runbm2geo
            %bm2gsettings.theElevations = [-70 -70 -70 -70];
            %bm2gsettings.theAzimuths = [270 90 0 180];
            bm2gcomment = sprintf('Transformed to earth coordinates by bm2geo with default transformation matrix;\n');
        end
        fprintf('%s\n',bm2gcomment);

        % code from runbm2g

        % now get the rest of the attributes that vary with experiment/site
        fprintf('The ADCP is a %s\n',settings.ADCPtype);

        dlg=fopen(dlgFile);
        disp(['Obtaining Beam configuration information from ' dlgFile])
        while 1
            line=fgetl(dlg);
            
            names=findstr(line,'Elevation');
            if ~isempty(names)
                for ii=1:4;
                    line=fgetl(dlg);
                    theBeams(ii)=str2double(line(3));
                    theElevations(ii)=str2num(line(13:20));
                    theAzimuths(ii)=str2num(line(23:30));
                end
                break;
            end
        end
        fclose(dlg);

        %for broad band (BB) has a perfect beam configuration
        % or for someone who just wants to use the default and didn't supply it...
        if ~exist('theElevations','var') || isempty(theElevations),
            theElevations = [-70 -70 -70 -70];
            fprintf('theElevations are missing: Using default elevations\n');
        end
        if ~exist('theAzimuths','var') || isempty(theAzimuths),
            theAzimuths = [270 90 0 180];
            fprintf('theAzimuths are missing: Using default azimuths\n');
        end
        fprintf('The Elevations are %4.1f %4.1f %4.1f %4.1f\n',theElevations);
        fprintf('The Azimuths are %4.1f %4.1f %4.1f %4.1f\n',theAzimuths);



return

function G = define_epFile(epDataFile, B, nbursts, nsamples)
% **********************************************************************
% Create the new output file.
G = netcdf(epDataFile,'clobber');

nbeams = 4;
theFillValue = fillval(B{'vel1'});

%copy the globals from the Beam File
disp('Copying the global attributes');
A=att(B);
% mlint doesn't like this line.  Don't touch it.
G < A;

%change and add the following attributes and variables
G.CREATION_DATE = datestr(now,0);
G.initial_instrument_height = ncfloat(B{'D'}.transducer_offset_from_bottom(1));
G.initial_instrument_height_note = 'height in meters above bottom: accurate for tripod mounted instruments';
% this one is added for readability, the old one is being kept for
% backwards compatability
G.transducer_offset_from_bottom = ncfloat(G.transducer_offset_from_bottom(1));

%RENAME A FEW
pulse = G.transmit_pulse_length;
rename(pulse,'transmit_pulse_length_cm');
pulse = G.ADCP_serial_number;
rename(pulse,'serial_number');

%Define some new ones for Epic compatability - these are already known
% these dont' change
G.transform = 'EARTH';
% this is OK, keep ADCP as generic designation, acoustic Doppler current
% profile would apply to other instruments Nortek Aquadopps, Sontek ADPs
% that use other designations to differentiate themselves from TRDI ADCP
G.DATA_TYPE = 'ADCP';
G.DATA_SUB_TYPE = 'MOORED';
G.COORD_SYSTEM = 'GEOGRAPHIC';
G.WATER_MASS = ncchar('?');
G.POS_CONST = nclong(0);  %1 if consistent
G.DEPTH_CONST = nclong(0);  %1 if consistent
G.DRIFTER = nclong(0);


%define the dimensions
G('time') = 0; %nbursts;
G('depth') = length(B{'D'});
G('lon') = 1;
G('lat') = 1;
G('profile') = nsamples;

%Variables and Attributes
disp('Defining variables and their attributes')

% 624:TIM:EPIC SYSTEM TIME         :time: : :2 integers, jul day, millisec since midnight
G{'time'} = nclong('time','profile');
G{'time'}.FORTRAN_format = ncchar('F10.2');
G{'time'}.units = ncchar('True Julian Day');
G{'time'}.type = ncchar('UNEVEN');
G{'time'}.epic_code = nclong(624);
G{'time'}.datum = ncchar('Time (UTC) in True Julian Days: 2440000 = 0000 h on May 23, 1968');
G{'time'}.NOTE = ncchar('Decimal Julian day [days] = time [days] + ( time2 [msec] / 86400000 [msec/day] )');
% no fill value for time - it breaks code.
%G{'time'}.FillValue_ = theFillValue;
G{'time2'} = nclong('time','profile') ;
G{'time2'}.FORTRAN_format = ncchar('F10.2');
G{'time2'}.units = ncchar('msec since 0:00 GMT');
G{'time2'}.type = ncchar('UNEVEN');
G{'time2'}.epic_code = nclong(624);
G{'time2'}.datum = ncchar('Time (UTC) in True Julian Days: 2440000 = 0000 h on May 23, 1968');
G{'time2'}.NOTE = ncchar('Decimal Julian day [days] = time [days] + ( time2 [msec] / 86400000 [msec/day] )');
% CRS 8/17/09 do not include with time2
%G{'time2'}.datum = ncchar('Time (UTC) in USGS convention: 2440000 = 0000 h on May 23, 1968');
%G{'time2'}.FillValue_ = theFillValue;

% 3:D  :DEPTH (M)                :depth:m:f10.1:
G{'depth'} = ncfloat('depth'); ncobj = G{'depth'};
% copy old attributes, omitting some
theAtts = att(B{'D'});
for iAtt = 1:length(theAtts),
    copy(theAtts{iAtt},ncobj);
end

disp('...Defining "bindist"')
G{'bindist'} = ncfloat('depth');
ncobj = G{'bindist'};
ncobj.FORTRAN_format = ncchar('f10.2');
ncobj.units = ncchar('m');
ncobj.type = ncchar('EVEN');
ncobj.epic_code = nclong(0);
ncobj.long_name = ncchar('bin distance from instrument (m)');
ncobj.name = ncchar('bindist');
ncobj.generic_name = ncchar('bindist');
ncobj.bin_size = B{'D'}.bin_size(:);
ncobj.transducer_offset_from_bottom = ncfloat(B{'D'}.transducer_offset_from_bottom(1));
ncobj.center_first_bin = ncfloat(B{'D'}.center_first_bin(1));
ncobj.blanking_distance = ncfloat(B{'D'}.blanking_distance(1));
ncobj.sensor_depth = ncfloat(theFillValue);
ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
ncobj.FillValue_ = theFillValue;
ncobj.NOTE = ncchar('distance is along profile from instrument head to center of bin');

disp('...Defining "lon"')
G{'lon'} = ncfloat('lon'); %% 1 element.
G{'lon'}.FORTRAN_format = ncchar('f10.4');
G{'lon'}.units = ncchar(B.longitude_units(:));
G{'lon'}.type = ncchar('EVEN');
G{'lon'}.epic_code = nclong(502);
G{'lon'}.name = ncchar('LON');
G{'lon'}.long_name = ncchar('LONGITUDE');
G{'lon'}.generic_name = ncchar('lon');
G{'lon'}.datum = ncchar('NAD83');
%G{'lon'}.FillValue_ = theFillValue;
if ischar(G.longitude(:)),
    G{'lon'}(1) = str2double(G.longitude(:));
else
    G{'lon'}(1) = G.longitude(:);
end

disp('...Defining "lat"')
G{'lat'} = ncfloat('lat'); %% 1 element.
G{'lat'}.FORTRAN_format = ncchar('F10.2');
G{'lat'}.units = ncchar(B.latitude_units(:));
G{'lat'}.type = ncchar('EVEN');
G{'lat'}.epic_code = nclong(500);
G{'lat'}.name = ncchar('LAT');
G{'lat'}.long_name = ncchar('LATITUDE');
G{'lat'}.generic_name = ncchar('lat');
G{'lat'}.datum = ncchar('NAD83');
%G{'lat'}.FillValue_ = theFillValue;
if ischar(G.latitude(:)),
    G{'lat'}(1) = str2double(G.latitude(:));
else
    G{'lat'}(1) = G.latitude(:);
end

G{'burst'} = ncfloat('time'); 
G{'burst'}.units = ncchar('count');
G{'burst'}.name = ncchar('burst');
G{'burst'}.generic_name = ncchar('record');
G{'burst'}.epic_code = nclong(1207);
G{'burst'}.long_name = ncchar('Burst Number');

% replaced by profile to match other burst files
% disp('...Defining "Rec"')
% G{'Rec'} = ncfloat('time','profile');
% ncobj = G{'Rec'};
% ncobj.FORTRAN_format = ncchar('f10.2');
% ncobj.units = ncchar('count');
% ncobj.type = ncchar('EVEN');
% ncobj.epic_code = nclong(0);
% ncobj.long_name = ncchar('ensemble number');
% ncobj.name = ncchar('ensemble');
% ncobj.generic_name = ncchar('ensemble');
% ncobj.FillValue_ = theFillValue;

G{'profile'} = ncfloat('time','profile'); ncobj = G{'profile'};
ncobj.units = ncchar('count');
ncobj.generic_name = ncchar('record');
ncobj.epic_code = nclong(1207);
ncobj.long_name = ncchar('Profile Number'); % this is the same sense as Sontek

disp('...Defining "u_1205"')
G{'u_1205'} = ncfloat('time', 'profile', 'depth');
G{'u_1205'}.name = ncchar('u');
G{'u_1205'}.long_name = ncchar('Eastward Velocity');
G{'u_1205'}.generic_name = ncchar('u');
G{'u_1205'}.FORTRAN_format = ncchar(' ');
G{'u_1205'}.units = ncchar('cm/s');
G{'u_1205'}.epic_code = nclong(1205);
G{'u_1205'}.sensor_type = B.INST_TYPE(:);
G{'u_1205'}.sensor_depth = ncfloat(theFillValue);
G{'u_1205'}.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
G{'u_1205'}.height_depth_units = ncchar('m');
G{'u_1205'}.serial_number = B.ADCP_serial_number(:);
G{'u_1205'}.minimum = ncfloat(theFillValue);
G{'u_1205'}.maximum = ncfloat(theFillValue);
%Not sure this is the right valid range, but what Fran had
G{'u_1205'}.valid_range = ncfloat([-1000 1000]);
G{'u_1205'}.FillValue_ = theFillValue;

disp('...Defining "v_1206"')
G{'v_1206'} = ncfloat('time', 'profile', 'depth');
G{'v_1206'}.name = ncchar('v');
G{'v_1206'}.long_name = ncchar('Northward Velocity');
G{'v_1206'}.generic_name = ncchar('v');
G{'v_1206'}.FORTRAN_format = ncchar(' ');
G{'v_1206'}.units = ncchar('cm/s');
G{'v_1206'}.epic_code = nclong(1206);
G{'v_1206'}.sensor_type = B.INST_TYPE(:);
G{'v_1206'}.sensor_depth = ncfloat(theFillValue);
G{'v_1206'}.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
G{'u_1206'}.height_depth_units = ncchar('m');
G{'v_1206'}.serial_number = B.ADCP_serial_number(:);
G{'v_1206'}.minimum = ncfloat(theFillValue);
G{'v_1206'}.maximum = ncfloat(theFillValue);
G{'v_1206'}.valid_range = ncfloat([-1000 1000]);
G{'v_1206'}.FillValue_ = theFillValue;

disp('...Defining "w_1204"')
G{'w_1204'} = ncfloat('time', 'profile', 'depth');
G{'w_1204'}.name = ncchar('w');
G{'w_1204'}.long_name = ncchar('Vertical Velocity');
G{'w_1204'}.generic_name = ncchar('w');
G{'w_1204'}.FORTRAN_format = ncchar(' ');
G{'w_1204'}.units = ncchar('cm/s');
G{'w_1204'}.epic_code = nclong(1204);
G{'w_1204'}.sensor_type = B.INST_TYPE(:);
G{'w_1204'}.sensor_depth = ncfloat(theFillValue);
G{'w_1204'}.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
G{'u_1204'}.height_depth_units = ncchar('m');
G{'w_1204'}.serial_number = B.ADCP_serial_number(:);
G{'w_1204'}.minimum = ncfloat(theFillValue);
G{'w_1204'}.maximum = ncfloat(theFillValue);
G{'w_1204'}.valid_range = ncfloat([-1000 1000]);
G{'w_1204'}.FillValue_ = theFillValue;

disp('...Defining "Werr_1201"')
G{'Werr_1201'} = ncfloat('time', 'profile', 'depth');
G{'Werr_1201'}.name = ncchar('Werr');
G{'Werr_1201'}.long_name = ncchar('Error Velocity');
G{'Werr_1201'}.generic_name = ncchar('w');
G{'Werr_1201'}.FORTRAN_format = ncchar('F8.1');
G{'Werr_1201'}.units = ncchar('cm/s');
G{'Werr_1201'}.epic_code = nclong(1201);
G{'Werr_1201'}.sensor_type = B.INST_TYPE(:);
G{'Werr_1201'}.sensor_depth = ncfloat(theFillValue);
G{'Werr_1201'}.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
G{'Werr_1201'}.height_depth_units = ncchar('m');
G{'Werr_1201'}.serial_number = B.ADCP_serial_number(:);
G{'Werr_1201'}.minimum = ncfloat(theFillValue);
G{'Werr_1201'}.maximum = ncfloat(theFillValue);
G{'Werr_1201'}.valid_range = B.error_velocity_threshold(:);
G{'Werr_1201'}.FillValue_ = theFillValue;

disp('...Defining "AGC_1202"')
for ibeam = 1:nbeams,
    vname = sprintf('AGC%1d_1202',ibeam);
    G{vname} = ncfloat('time','profile', 'depth');
    ncobj = G{vname};
    ncobj.name = ncchar(sprintf('AGC%1d',ibeam));
    ncobj.long_name = ncchar(sprintf('Beam %1d Average Echo Intensity (AGC)',ibeam));
    ncobj.generic_name = ncchar('AGC');
    ncobj.FORTRAN_format = ncchar('F5.1');
    ncobj.units = ncchar('counts');
    ncobj.epic_code = nclong(1202);
    ncobj.sensor_type = B.INST_TYPE(:);
    ncobj.sensor_depth = ncfloat(theFillValue);
    ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
    ncobj.height_depth_units = ncchar('m');
    ncobj.serial_number = B.ADCP_serial_number(:);
    ncobj.norm_factor = B{sprintf('AGC%1d',ibeam)}.norm_factor(:);
    ncobj.NOTE = ncchar('normalization to db');
    ncobj.minimum = ncfloat(theFillValue);
    ncobj.maximum = ncfloat(theFillValue);
    ncobj.valid_range = B.false_target_reject_values(:);
    ncobj.FillValue_ = theFillValue;
end

%1203:PGd :Percent Good Pings       : : : :ADCP Percent Good Pings
for ibeam = 1:nbeams,
    vname = sprintf('PGd%1d_1203',ibeam);
    G{vname} = ncfloat('time','profile', 'depth');
    ncobj = G{vname};
    ncobj.name = ncchar(sprintf('PGd%1d',ibeam));
    ncobj.long_name = ncchar(sprintf('Beam %1d Percent Good Pings',ibeam));
    ncobj.generic_name = ncchar('PGd');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.units = ncchar('percent'); % Workhorse Command and Output Data Format manual
    ncobj.epic_code = nclong(1203);
    ncobj.sensor_type = B.INST_TYPE(:);
    ncobj.sensor_depth = ncfloat(theFillValue);
    ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
    ncobj.height_depth_units = ncchar('m');
    ncobj.serial_number = B.ADCP_serial_number(:);
    ncobj.minimum = ncfloat(theFillValue);
    ncobj.maximum = ncfloat(theFillValue);
    ncobj.valid_range = B.minmax_percent_good(:);
    ncobj.FillValue_ = theFillValue;
end

if ~isempty(B{'brange'}(:)),
    disp('...Defining "brange"')
    G{'brange'} = ncfloat('time', 'profile');
    ncobj = G{'brange'};
    ncobj.name = ncchar('brange');
    ncobj.long_name = ncchar('range from instrument head to boundary');
    ncobj.generic_name = ncchar('brange');
    ncobj.FORTRAN_format = ncchar('f10.2');
    ncobj.units = ncchar('m');
    ncobj.epic_code = nclong(0);
    ncobj.sensor_type = B.INST_TYPE(:);
    ncobj.sensor_depth = ncfloat(theFillValue);
    ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
    ncobj.height_depth_units = ncchar('m');
    ncobj.minimum = ncfloat(theFillValue);
    ncobj.maximum = ncfloat(theFillValue);
    ncobj.serial_number = B.ADCP_serial_number(:);
    switch B.frequency(1), % this is more realistic
        case 1200
            ncobj.valid_range = ncfloat([0 50]);
        case 600
            ncobj.valid_range = ncfloat([0 100]);
        case 300
            ncobj.valid_range = ncfloat([0 200]);
        otherwise
            ncobj.valid_range = ncfloat([0 1000]);
    end
    ncobj.FillValue_ = theFillValue;
    ncobj.NOTE = B{'brange'}.NOTE(:);
else
    disp('No range to boundary data, skipping brange')
end

% 1211:Tx  :ADCP Transducer Temp.    :temp:deg. C:F10.2:ADCP Transducer Temp
disp('...Defining "Tx_1211"')
G{'Tx_1211'} = ncfloat('time', 'profile');
G{'Tx_1211'}.name = ncchar('Tx');
G{'Tx_1211'}.long_name = ncchar('instrument Transducer Temp.');
G{'Tx_1211'}.generic_name = ncchar('temp');
G{'Tx_1211'}.units = ncchar('degrees.C');
G{'Tx_1211'}.epic_code = nclong(1211);
G{'Tx_1211'}.sensor_type = B.INST_TYPE(:);
G{'Tx_1211'}.sensor_depth = ncfloat(theFillValue);
G{'Tx_1211'}.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
G{'Tx_1211'}.height_depth_units = ncchar('m');
G{'Tx_1211'}.serial_number = B.ADCP_serial_number(:);
G{'Tx_1211'}.minimum = ncfloat(theFillValue);
G{'Tx_1211'}.maximum = ncfloat(theFillValue);
G{'Tx_1211'}.valid_range = B{'Tx'}.valid_range(:);
G{'Tx_1211'}.FillValue_ = theFillValue;

% consider if salinity is constant and a sensor is present
if strcmpi(B.conductivity_sensor(:),'NO'),
    disp('No conductivity sensor attached.')
    fprintf('Setting salinity_set_by_user to %f\n',B{'S'}(1));
    G.salinity_set_by_user = ncfloat(B{'S'}(1));
    G.salinity_set_by_user_units = ncchar('PPT');
elseif strcmpi(B.conductivity_sensor(:),'YES'),
    % 40:S  :SALINITY (PPT)           :sal:PPT:f10.2:From SALIN1
    disp('Conductivity sensor attached.')
    disp('...Defining "S_40"')
    G{'S_40'} = ncfloat('time', 'profile');
    ncobj = G{'S_40'};
    ncobj.name = ncchar('S');
    ncobj.long_name = ncchar('SALINITY (PPT)');
    ncobj.generic_name = ncchar('sal');
    ncobj.units = ncchar('PPT');
    ncobj.epic_code = nclong(40);
    ncobj.sensor_type = B.INST_TYPE(:);
    ncobj.sensor_depth = ncfloat(theFillValue);
    ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
    ncobj.height_depth_units = ncchar('m');
    ncobj.serial_number = B.ADCP_serial_number(:);
    ncobj.minimum = ncfloat(theFillValue);
    ncobj.maximum = ncfloat(theFillValue);
    ncobj.valid_range = B{'S'}.valid_range(:);
    ncobj.FillValue_ = theFillValue;
end

% 80:SV :SOUND VELOCITY (M/S)     : :m s-1: :
disp('...Defining "SV_80"')
G{'SV_80'} = ncfloat('time',  'profile');
ncobj = G{'SV_80'};
ncobj.name = ncchar('SV');
ncobj.long_name = ncchar('SOUND VELOCITY (M/S)');
ncobj.generic_name = ncchar(' ');
ncobj.units = ncchar('m s-1');
ncobj.epic_code = nclong(80);
ncobj.sensor_type = B.INST_TYPE(:);
ncobj.sensor_depth = ncfloat(theFillValue);
ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
ncobj.height_depth_units = ncchar('m');
ncobj.serial_number = B.ADCP_serial_number(:);
ncobj.minimum = ncfloat(theFillValue);
ncobj.maximum = ncfloat(theFillValue);
ncobj.valid_range = ncfloat(B{'sv'}.valid_range(:));
ncobj.FillValue_ = theFillValue;

%TODO
% copy(v, myVariable, copyData, copyAttributes) copies the
%           contents of variable v into variable myVariable, as
%           directed by the flags (defaults = 0).
disp('...Defining "Hdg_1215"')
% 1215:Hdg :INST Heading             :hdg:degrees:F10.2:ADCP heading
G{'Hdg_1215'} = ncfloat('time', 'profile');
ncobj = G{'Hdg_1215'};
ncobj.name = ncchar('Hdg');
ncobj.long_name = ncchar('INST Heading             ');
ncobj.generic_name = ncchar('hdg');
ncobj.epic_code = nclong(1215);
ncobj.units = ncchar('degrees');
ncobj.FORTRAN_format = ncchar('f10.2');
ncobj.sensor_type = B.INST_TYPE(:);
ncobj.sensor_depth = ncfloat(theFillValue);
ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom(1));
ncobj.height_depth_units = ncchar('m');
ncobj.serial = B.ADCP_serial_number(:);
ncobj.minimum = ncfloat(theFillValue);
ncobj.maximum = ncfloat(theFillValue);
ncobj.valid_range = ncfloat([0 360]);
ncobj.FillValue_ = theFillValue;
ncobj.heading_alignment_EH = ncfloat(B{'Hdg'}.heading_alignment(1));
ncobj.heading_bias_EB = ncfloat(B{'Hdg'}.heading_bias(1));
ncobj.NOTE = ncchar(B{'Hdg'}.NOTE_9(:));

disp('...Defining "Ptch_1216"')
% 1216:Ptch:INST Pitch               :ptch:degrees:F10.2:ADCP pitch
G{'Ptch_1216'} = ncfloat('time', 'profile');
ncobj = G{'Ptch_1216'};
ncobj.name = ncchar('Ptch');
ncobj.long_name = ncchar('INST Pitch             ');
ncobj.generic_name = ncchar('ptch');
ncobj.epic_code = nclong(1216);
ncobj.units = ncchar('degrees');
ncobj.FORTRAN_format = ncchar('f10.2');
ncobj.sensor_type = B.INST_TYPE(:);
ncobj.sensor_depth = ncfloat(theFillValue);
ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
ncobj.height_depth_units = ncchar('m');
ncobj.serial = B.ADCP_serial_number(:);
ncobj.minimum = ncfloat(theFillValue);
ncobj.maximum = ncfloat(theFillValue);
ncobj.valid_range = ncfloat([-180 180]);
ncobj.FillValue_ = theFillValue;

disp('...Defining "Roll_1217"')
% 1217:Roll:INST Roll                :roll:degrees:F10.2:ADCP roll
G{'Roll_1217'} = ncfloat('time','profile');
ncobj = G{'Roll_1217'};
ncobj.name = ncchar('Roll');
ncobj.long_name = ncchar('INST Roll             ');
ncobj.generic_name = ncchar('roll');
ncobj.epic_code = nclong(1217);
ncobj.units = ncchar('degrees');
ncobj.FORTRAN_format = ncchar('f10.2');
ncobj.sensor_type = B.INST_TYPE(:);
ncobj.sensor_depth = ncfloat(theFillValue);
ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
ncobj.height_depth_units = ncchar('m');
ncobj.serial = B.ADCP_serial_number(:);
ncobj.minimum = ncfloat(theFillValue);
ncobj.maximum = ncfloat(theFillValue);
ncobj.valid_range = ncfloat([-180 180]);
ncobj.FillValue_ = theFillValue;

% 1294:P:ADCP Pressure:pres:deca-pascals:F15.2:ADCP Pressure at Transducer Head
if ~isempty(B{'Pressure'}),
    disp('...Defining "P_1294"')
    G{'P_1294'} = ncfloat('time', 'profile');
    ncobj = G{'P_1294'};
    ncobj.name = ncchar('P');
    ncobj.long_name = ncchar('Pressure at Transducer Head');
    ncobj.generic_name = ncchar('pres');
    ncobj.units = ncchar('deca-pascals');
    ncobj.FORTRAN_format = ncchar('f15.2');
    ncobj.epic_code = nclong(1294);
    ncobj.sensor_type = B.INST_TYPE(:);
    ncobj.sensor_depth = ncfloat(theFillValue);
    ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
    ncobj.height_depth_units = ncchar('m');
    ncobj.serial_number = B.ADCP_serial_number(:);
    ncobj.minimum = ncfloat(theFillValue);
    ncobj.maximum = ncfloat(theFillValue);
    ncobj.valid_range = B{'Pressure'}.valid_range(:);
    ncobj.FillValue_ = theFillValue;
    if ~isempty(B{'Pressure'}.AZ_applied(:)),
        ncobj.AZ_applied = ncchar(B{'Pressure'}.AZ_applied(:));
        if strcmpi(B{'Pressure'}.AZ_applied(:),'yes'), 
            ncobj.NOTE = ncchar('sensor zeroed before deployment');
            ncobj.comment = ncchar('pressure does not include atmospheric pressure');
        else
            ncobj.NOTE = ncchar('sensor not zeroed before deployment');
            ncobj.comment = ncchar('pressure includes atmospheric pressure and unknown sensor drift');
        end
    else ncobj.AZ_applied = ncchar('unknown');
        ncobj.NOTE = ncchar('unknown if pressure was zeroed before deployment');
        ncobj.comment = ncchar('Pressure may or may not include atmospheric pressure and sensor drift');
    end
else
    disp('No pressure data...skipping "P_1294"')
end

% 850:SDP:STAND. DEV. (PRESS)     :pres:mbar:f10.5:std. deviation of burst pressures
if ~isempty(B{'PressVar'}),
    disp('...Defining "SDP_850"')
    % 850:SDP:STAND. DEV. (PRESS)     :pres:mbar:f10.5:std. deviation of burst pressures
    G{'SDP_850'} = ncfloat('time','profile');
    ncobj = G{'SDP_850'};
    ncobj.epic_code = nclong(850);
    ncobj.name = ncchar('SDP');
    ncobj.long_name = ncchar('STAND. DEV. (PRESS)     ');
    ncobj.generic_name = ncchar('pres');
    ncobj.units = ncchar('mbar');
    ncobj.sensor_type = B.INST_TYPE(:);
    ncobj.sensor_depth = ncfloat(theFillValue);
    ncobj.initial_sensor_height = ncfloat(B{'D'}.transducer_offset_from_bottom (1));
    ncobj.height_depth_units = ncchar('m');
    ncobj.serial_number = B.ADCP_serial_number(:);
    ncobj.minimum = ncfloat(theFillValue);
    ncobj.maximum = ncfloat(theFillValue);
    ncobj.valid_range = ncfloat([0 5000]);
    ncobj.FillValue_ = theFillValue;
end

endef(G)
return