function Pcadp = pcadp2nc(cdfbFile, cdfsFile, ncFile, switches, bursts)
% pcadp2nc - is a script that reads pcadp data from a netcdf file into matlab, 
%  applies rotations and outputs epic compliant .nc output
% you will want to convert from a .adp file to .cdf using adp2cdf.m
% this writes velocity, pressure & temperature data only
% using the batch capability will save lots of interactive input
%
% usage: Pcadp = pcadp2nc(cdfbFile, cdfsFile, outfile, switches, bursts)
%   cdfbFile = raw PCADP data netCDF input file of burst data *b.cdf
%   cdfsFile = raw PCADP data netCDF input file of statistics data *s.cdf
%   outfile = file name for the netCDF output file, to which b.nc is appended 
%       for the burst data, which is not EPIC compliant, shaped [time, profile, depth]; 
%       and *s-cal.nc for the statistics, in EPIC compliant format
% switches = structure of controls as to how data are processed, defaults shown
%   switches.ecorr_showplots = 0;
%   switches.resolve_ambiguity = 0 do not resolve, 1 resolve
%   switches.amb.showplots = 0 suppress plots, 1 show plots (this is useful for 1
%           or two bursts that you are trying to clean up)
%   switches.amb.timecorrect = 1 use time correction, 0 suppress (use this
%           if your profiles are more than 5 sec apart in time)
%   switches.amb.cut_factor = 0.7, change at your own risk...;
%   switches.amb = []; use all defaults for ambiguity correction
%   switches.range_correct = 1; adjust distance to bottom measured by beam
%       range for pitch and roll
%   The following can be a per burst vector omitted entirely
%       Vector length must match the number of bursts this function is
%       asked to process
%       switches.override_heading = []; 
%       switches.override_pitch = []; 
%       switches.override_roll = []; 
% bursts = [startburst endburst], as Sontek burst numbers-
%   note this is NOT the burst index!  Use browsehydraburst to find the
%   first good burst, then add from there, if you just want to do a few.
% Pcadp = output structure of data using Lacy format, the last burst read
%

warning off MATLAB:divideByZero

% 13 Oct 09 (MM) found some nc. attribute references where that infor
% hadn't been defined yet, fixed that and some m-lint stuff & the status
% variable, which is no longer used.
% 11 Aug 08
% -- added bindist to the variables
% -- fixed missing SN for all PCADP variables, xducer -> transducer
% -- changed corr counts to 'percent'
% 4 Nov 07
% -- trap headings of < 0 after declination correction and fix.
% 30 May 07 
% -- fix typo "Pcadp.cdf" which would have resulted in missind conductivity data
% THanks Kurt Rosenberger!
% 21 Nov 05
% -- change the input bursts to reference Sontek burst numbers, not file
% indeces
% 25 oct 05
% -- add acoustic backscatter intensity output to nc files
% -- add bottom range correction for tilt, stats file only
% 24 oct 05 version 2.2
% -- add resolution velocity output to burst nc file
% 06 oct 05
% -- remove Qa output, now handled by quality file
% 19 sep 05
% -- remove and make private subfunctions that are the same for ADV and
% PCADP
% 7 Jul 2005 version 2.1
% -- Ures is stored in .cdf file in cm/s, make sure that is OK with this code
% -- fix bug where sed was not getting written when coef were supplied
% 27 Apr 2005
% -- replace pcoord with uv2polar which comes with the toolbox
% 7 Apr 2005
% -- fix some typos in the internal pressure code
% 29 Mar 2005
% -- add code for calibrated OBS data
% 25 Mar 2005
% -- replace variable attribute serial_number with serial to be consistent 
% with .cdf file, fix typo that was causing missing serial numbers
% 22 Mar 2005
% -- fix depth and height units, all as m
% -- Add Rachel's initial_sensor_height attributes
% 3 Mar 2005
% -- fix start and stop time computation
% 12/15/04
% -- add pressure mask option
% 15 Oct 2004
% -- after a data meeting finalaize variable names (again)
% 10/12/04 fix EPIC issues: temperature names, history storage
% this involved many variable name changes
% 9/28/04 got ambiguity correction working.  It will not work if certain
% variables in the stats file are not available if there are bad bursts,
% these bursts are therefore will not have ambiguity correction applied.
% 9/27/04 found and fixed some bugs
% 9/26/04
% -- check depth and height calculation, extend height calculateion to a
% height time seriea for each cell
% 8/25/04
% -- take height, range to bed from the stats file instead.  This is likely
% to be corrected by the user in between running adp2cdf anf pcadp2nc.
% -- put the correct expression of depth and height in the right places
% 8/20/04 Add median heading, pitch and roll to stats output
%       use gmax, gmin, gmean, gmedian and gstd throughout
% 6/23/04 Misc small fixes
% -- figure out where the real blank info comes from... done & tested using
% read data.  When the user sets the blank, it is stored in ProfBlankLenCoh
% 4/27/04 add pressure.  CTD will not be added here.  
%       Only PCADP sensors, and those one might require for burst statistics
%       heading, pitch, roll, amp, cor are not copied over to save disk space
%       only include them if they are needed for scientific analysis
% 4/6/04 change to read burst by burst and match the new shape of the data
% 12/15/03 got rotations reading right and quickly
% 12/03/03 figuring out Sontek terminology for PCADP settings
% 11/20/03-ish removed perl dependencies, now reads from a netcdf raw file
% this software is too path dependent.  Trying to make it work as a toolbox.
% Installation of the toolbox...
% take this toolbox and put it in a dirctory and put that directory ont he
% matlab path
% 11/12/03 first version
% Adapted by Marinna Martini for the U.S. Geological Survey 
% Coastal & Marine Program Woods Hole Field Center, Woods Hole, MA
% http://woodshole.er.usgs.gov/ Please report bugs to mmartini@usgs.gov
% based on code written by Jessie Lacy USGS Santa Cruz Science Center

% prerequisites & dependents
% netCDF toolbox
% downlooking PCADP data only
% xyz2enu.m         rotate velocities
% hydrapress2db.m   handler for frequency Druck or Paros sensors
% adp2cdf.m         convert raw Sontek binary adp files to raw cdf (engineering units).  No processing done
% readdruck.m       read pressure calibration file *.drk
% pcadpambiguity.m  apply ambiguity correction
% gmean.m           compute the mean while ignoring non-finite values
% add_vardesc.m

% tools which can see this data
% advbrowsegui.m for the burst data
% ncbrowser.m for the stats data

mversion = 'SVN $Revision$';

% ----- resolve inputs

if ~exist('cdfbFile','var'),
    [fname, pathname] = uigetfile('*b.cdf', 'Pick a netcdf file containing PCADP data');
    if isequal(fname,0) || isequal(pathname,0)
        disp('User pressed cancel')
        %status = 0;
        return
    end
    cdfbFile = fullfile(pathname, fname);
    disp(['User selected ', cdfbFile])
end
if ~exist('cdfsFile','var'),
    [fname, pathname] = uigetfile('*s.cdf', 'Pick a netcdf file containing PCADP data');
    if isequal(fname,0) || isequal(pathname,0)
        disp('User pressed cancel')
        %status = 0;
        return
    end
    cdfsFile = fullfile(pathname, fname);
    disp(['User selected ', cdfsFile])
end
if ~exist('ncFile','var') || isempty(ncFile),    
    prompt  = {'Enter an output file name to which b-cal.nc and vp-cal.nc will be appended:'}; 
    def = {'outputname'}; 
    title   = 'Output file name';
    lineNo  = 1;
    result = inputdlg(prompt,title,lineNo,def);
    ncFile = result{1};
end

% make sure this is a burst PCADP data file
cdfb = netcdf(cdfbFile,'nowrite');
if ~strcmp(cdfb.DESCRIPT(:),'Sontek PCADP raw data burst file'),
    disp('File is not a Sontek PCADP raw data netCDF file containing all profiles (*b.cdf)')
    close (cdfb)
    %status = -1;
    return
end
% make sure this is a statistics PCADP data file
cdfs = netcdf(cdfsFile,'nowrite');
if ~strcmp(cdfs.DESCRIPT(:),'Sontek PCADP raw data statistics file'),
    disp('File is not a Sontek PCADP raw data netCDF file containing statistics (*s.cdf)')
    close (cdfs)
    %status = -1;
    return
end
if ~exist('switches','var'),
    switches.resolve_ambiguity = 1;
    button = questdlg('Do you want to resolve ambiguity?',...
        'Ambiguity');
    if strcmp(button,'No'), switches.resolve_ambiguity = 0; end
    switches.show_amb_plots = 1;
    button = questdlg('Do you want to see the ambiguity plots?',...
        'Ambiguity');
    if strcmp(button,'No'), switches.show_amb_plots = 0; end
else
    if ~isfield(switches, 'resolve_ambiguity'), switches.resolve_ambiguity = 1; end
    if ~isfield(switches, 'show_amb_plots'), switches.show_amb_plots = 1; end
    if ~isfield(switches, 'range_correct'), switches.range_correct = 1; end
end

% end of user input section

% more basic checks
% 1) the data must be recorded in BEAM for this processing technique to
% work
if ~strcmp(cdfb.PCADPUserSetupCoordSystem(:),'Beam'),
    disp('File is not in beam coordinates.  This technique only works on beam data.')
    close(cdfb)
    %status = -1;
    return
end
% 2) the data must be collected in Pulse Coherent Mode, but don't terminate
if ~strcmp(cdfb.PCADPUserSetupPCMode(:),'Yes'),
    disp('Data was not collected in Pulse Coherent Mode.  This script written for PCMode = yes')
    disp('Aborting because corellations from PCMode are needed to process data using this technique')
    close(cdfb)
    %status = -1;
    return
end

% the user wants a limited number of bursts
if ~exist('bursts','var') || isempty(bursts),
    bursts = [1 Inf];
end
if find(bursts < 1), % no negative indeces
    disp('Cannot use negative indeces')
    close (cdfb)
    return
end
% translate the burst numbers given by the user into burst indeces into the
% netCDF file
idx = find(bursts(1) == cdfb{'burst'}(:),1,'first');
if isempty(idx), bursts(1) = 1; 
else bursts(1) = idx;
end
if isinf(bursts(2)) || (bursts(2) > length(recdim(cdfb))), 
    bursts(2) = length(recdim(cdfb)); 
else
    idx = find(bursts(2) == cdfb{'burst'}(:),1,'first');
    if isempty(idx), bursts(2) = length(recdim(cdfb)); 
    else bursts(2) = idx;
    end
end

fprintf('Processing bursts %d to %d in file %s\n',...
    cdfb{'burst'}(bursts(1)),cdfb{'burst'}(bursts(2)),name(cdfb))
if ~isfield(switches,'resolve_ambiguity'), switches.resolve_ambiguity = 0; end
if switches.resolve_ambiguity, disp('Resolving ambiguity'),
    if ~isfield(switches,'amb'), switches.amb = []; end
    if ~isempty(switches.amb),
        if ~isfield(switches.amb,'showplots'), amb.showplots = switches.amb.showplots; end
        if ~isfield(switches.amb,'timecorrect'), amb.timecorrect = switches.amb.timecorrect; end
        if ~isfield(switches.amb,'cut_factor'), amb.cut_factor = switches.amb.cut_factor; end
    else
        switches.amb.showplots = 0;
        switches.amb.timecorrect = 1;
        switches.amb.cut_factor = 0.7;
    end
else disp('Not resolving ambiguity')
end    
if ~isfield(switches,'ecorr_showplots'), switches.ecorr_showplots = 0; end

tic
% get the user setup things of import, translate them into Jessie's structure
% was read_ctl
disp('Getting metadata')
Pcadp.magnetic_variation = cdfb.magnetic_variation(:);
Pcadp.def_temp = cdfb.PCADPUserSetupTemp(:); % ok
Pcadp.def_sal = cdfb.PCADPUserSetupSal(:); % ok
Pcadp.def_csound = cdfb.PCADPUserSetupCW(:); % ok
Pcadp.ncells = cdfb.PCADPUserSetupNcells(:); % ok
Pcadp.nbeams = cdfb.PCADPProbeNbeams(:); % total number of beams.  There are 4 beam adps out there
Pcadp.cellsize = cdfb.PCADPUserSetupCellSize(:); % ok in m
Pcadp.blank_dist = cdfb.PCADPProbeProfBlankLenCoh(:)./100; % ok cm convert to m
%Pcadp.blank_dist = cdfb.PCADPProbeBlankDist(:); %the same as above
Pcadp.sensor_depth = cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:); % ok
Pcadp.avg_int = cdfb.PCADPUserSetupAvgInterval(:); % ok
Pcadp.p_int = cdfb.PCADPUserSetupProfileInterval(:); % ok
Pcadp.b_int = cdfb.PCADPUserSetupBurstInterval(:); % ok
Pcadp.ppb = cdfb.PCADPUserSetupProfilesPerBurst(:); % ok
Pcadp.pulse_length = cdfb.PCADPProbeProfPlenCoh(:)./100; % ok
% Lacy definition
if cdfb.cpuSoftWareVerNum(:) > 16.0,
    Pcadp.max_range = (Pcadp.cellsize*Pcadp.ncells + Pcadp.blank_dist); % ok
else
    Pcadp.max_range = (Pcadp.cellsize*Pcadp.ncells.*10 + Pcadp.blank_dist); % not sure
end

%Pcadp.max_range
% ProfLag and ResLag are the ranges (distance) for the profiling and resolution pulses respectively, 
% not the blanking distances.
Pcadp.user_pulse_lag = cdfb.PCADPProbeProfPlagCoh(:)./100; % ok, what the user requested
Pcadp.user_profile_lag = Pcadp.user_pulse_lag; % ok, what the user requested, better name at Jessie's request
Pcadp.system_lag = cdfb.PCADPProbeProfActualPlagCoh(:)./100; % ok, the the PCADP can really do
Pcadp.user_res_lag = cdfb.PCADPProbeResPlagCoh(:)./100; % ok, what the user requested
Pcadp.system_res_lag = cdfb.PCADPProbeActualResPlagCoh(:)./100; % ok, the the PCADP can really do

% this section checked OK 12/3/03 as per Lacy email
% from adpdatapc.c by Sontek
% tmp =   0.01*(AdpOperation.ResPlagCoh);
% fprintf(Fp,"\nUser Res Lag ----------- (m) ------> %.2f", tmp);
% tmp =   0.01*(AdpOperation.ActualResPlagCoh);
% fprintf(Fp,"\nSystem Res Lag --------- (m) ------> %.2f", tmp);
% tmp = 2* 1.25e-6 * UserSetup.Cw *UserSetup.Cw /
%      (tmp*ADP_FREQ[AdpConf.AdpType]);
% fprintf(Fp,"\nMaxVertVel (Res Lag) --- (m/s) ----> %.2f", tmp);
% tmp  = tmp / sin(DegToRad(0.1*AdpConf.SlantAngle));
% fprintf(Fp,"\nMaxHorizVel (Res Lag) -- (m/s) ----> %.2f", tmp);
% Lacy:  These correspond to  MaxVertVel(ResLag) -- (m/s) and MaxHorizVel(ResLag) - (m/s)  
% These are the ambiguity velocities based on the range of the resolution pulse, 
% not the profiling pulse, which is the same as the 2x the max. resolved velocities 
% after ambiguity error correction.  My original conclusion remains unchanged:
% Pcadp.max_vert_vel = MaxVertVel(ResLag)
% Pcadp.max_horiz_vel = MaxHorizVel(ResLag)
Pcadp.max_vert_vel = (2.*1.25e-6.*Pcadp.def_csound.^2)./...
     (Pcadp.system_res_lag.*cdfb.PCADPProbeFrequency(:)./100000) ; % from adpdatpc.c
% (2.*(1.25e-6).*(1508.1.^2))./(0.54.*1500000/100000) works
Pcadp.max_horiz_vel = Pcadp.max_vert_vel./sin(cdfb.SlantAngle(:)*pi/180); % from adpdatpc.c
% ((2.*(1.25e-6).*(1508.1.^2))./(0.54.*1500000/100000))/(sin(15*pi/180)) works
Pcadp.min_corr_level = cdfb.PCADPProbeProfMinCorrLevel(:); % OK

% These correspond to  MaxVertVel(ProfLag) -- (m/s) and MaxHorizVel(ProfLag) - (m/s)  
% not used by Jessie's code
Pcadp.max_vert_Prof_vel = (2.*1.25e-6.*Pcadp.def_csound.^2)./...
     (Pcadp.system_lag.*cdfb.PCADPProbeFrequency(:)./100000) ; % from adpdatpc.c
% % (2.*(1.25e-6).*(1508.1.^2))./(1.21.*1500000/100000) works
Pcadp.max_horiz_Prof_vel = Pcadp.max_vert_Prof_vel./sin(cdfb.PCADPProbeSlantAngle(:)*pi/180); % from adpdatpc.c
% % ((2.*(1.25e-6).*(1508.1.^2))./(1.21.*1500000/100000))/(sin(15*pi/180)) works


Pcadp.av = 2.*Pcadp.max_horiz_vel; % Lacy: Along beam ambiguity velocity (2X max. unambiguous vel.) for the resolution velocities, cm/s
Pcadp.cell_start = cdfb.PCADPProbeResBlankLenCoh(:)./100; % OK?, Blanking Distance for the resolution cell (actual), m
Pcadp.cell_length = cdfb.PCADPProbeResPlenCoh(:)./100; % Resolution Pulse Length, not controlled by the user, m, should be 0.1 m
% but I think this is the correct one, because ResPlenCoh reads out of the binary file as 0.02 m
Pcadp.cell_length = cdfb.PCADPProbeResCellSize(:)./100; % seems more straightforward and = 0.1 m
Pcadp.pr_plag = Pcadp.system_lag; % redundant
Pcadp.avpl = 2.*Pcadp.max_horiz_Prof_vel; % Lacy: Along beam ambiguity velocity (2X max. unambiguous vel.) for the profile velocities, cm/s
Pcadp.prange_start = Pcadp.blank_dist; % Profile blanking distance

% Jessie's Ambiguity routine wants >400 profiles per burst...
if switches.resolve_ambiguity && (cdfb.PCADPUserSetupProfilesPerBurst(:) <= 400),
    fprintf('There are %d profiles per burst, less than the 401 desired by pcadpambiguity.\n', cdfb.PCADPUserSetupProfilesPerBurst(:))
    disp('Ambiguity correction will be performed anyway, but may not be as accurate.')
end

% things needed for pressure
% p is a separate struct with pressure specific info to be passed to the pressure normalization function
% build the static part of the pressure data structure needed for hydrapress2db
if ~isempty(cdfb{'extpressfreq'}) || ~isempty(cdfb{'extpress'}) || ~isempty(cdfb{'pressure'}),
    P.ExtPressInstalled = cdfb.ExtPressInstalled(:); % 'No' | 'Paros' | 'Druck' | 'ParosFreq'
    P.PressOffset = cdfb.PressOffset(:);
    P.PressScale = cdfb.PressScale(:);
    P.PressScale_2 = cdfb.PressScale_2(:);
    P.PressFreqOffset = cdfb.PressFreqOffset(:);
    P.INST_TYPE = cdfb.INST_TYPE(:);
    % frequency sensor case
    if strcmp(P.ExtPressInstalled,'Druck') || strcmp(P.ExtPressInstalled,'ParosFreq'),
        P.presscals = get_calvalues(cdfb{'extpressfreq'}); 
        if isempty(P.presscals), 
            disp('Expected to find frequency pressure calibrations and did not.')
            disp('Pressure will not be processed')
        end
    else P.presscals = [];
    end
else
    P = [];
end

% things needed for external sensors
% build the static part of the optical data
if ~isempty(cdfb{'extsensor1'}) 
    extmeta{1} = getmetafromatts(cdfb{'extsensor1'});
end
if ~isempty(cdfb{'extsensor2'}) 
    extmeta{2} = getmetafromatts(cdfb{'extsensor2'});
end

% need to improve how height off bottom and depth are done
% fix the elevations, using all data present in the stats file no matter
% how many bursts are being worked on in the burst file
% range data only in PCMode header... so often not useful

%this is the original method
for beam = 1:3, % beams 1 - 3
    varname = sprintf('MeanRange%1d',beam);
    Pcadp.range(:,beam) = cdfs{varname}(bursts(1):bursts(2))./1000; % range to bottom (m) converted from mm
    Pcadp.range(:,beam) = ecorr(Pcadp.range(:,beam),switches.ecorr_showplots);
end

%   switches.range_correct = 1; adjust distance to bottom measured by beam
if switches.range_correct,
    Pcadp.range = vert_bd(Pcadp.range, cdfs{'MedianPitch'}(bursts(1):bursts(2)), cdfs{'MedianRoll'}(bursts(1):bursts(2)));
end
Pcadp.meanrange = gmean(Pcadp.range);

% meta data to pass through
lon = cdfb.longitude(:);
lat = cdfb.latitude(:); 

% define the .nc file
ncs = definencsfile(ncFile, cdfb, cdfs); % define the file which will contain burst stats
ncs{'lon'}(:) = lon;
ncs{'lat'}(:) = lat; 

ncb = definencbfile(ncFile, cdfb, bursts); % define the file which will contain the bursts
ncb{'lon'}(:) = lon;
ncb{'lat'}(:) = lat; 
clear lat lon

% ------------------ now process everything else .... operate burst by burst --------------
%
%
% **************** I have not used autonaning here yet, may want to
% make holders for things to increase speed
Pcadp.nbursts = length(recdim(cdfb)); 
Pcadp.nprof = Pcadp.nbursts.*Pcadp.ppb;   % total number of profiles
% pre define the large arrays for speed
Pcadp.num = ones(Pcadp.ppb,1).*NaN;
Pcadp.jt = Pcadp.num;
Pcadp.nsamp = Pcadp.num;
Pcadp.c = Pcadp.num;
Pcadp.hdg = Pcadp.num;
Pcadp.pitch = Pcadp.num;
Pcadp.roll = Pcadp.num;
Pcadp.t = Pcadp.num;
Pcadp.p = Pcadp.num;
Pcadp.hdg_sd = Pcadp.num;
Pcadp.pitch_sd = Pcadp.num;
Pcadp.roll_sd = Pcadp.num;
Pcadp.t_sd = Pcadp.num;
Pcadp.p_sd = Pcadp.num;
Pcadp.battV = Pcadp.num;
Pcadp.ad1 = Pcadp.num;
Pcadp.ad2 = Pcadp.num;
Pcadp.ad1_sd = Pcadp.num;
Pcadp.ad2_sd = Pcadp.num;
Pcadp.extpress = Pcadp.num;
Pcadp.pfreq = Pcadp.num;
Pcadp.bot_dist = ones(Pcadp.ppb,Pcadp.nbeams).*NaN;
Pcadp.res_vel = Pcadp.bot_dist;
Pcadp.res_cor = Pcadp.bot_dist;
nandata = ones(Pcadp.ppb,Pcadp.ncells).*NaN;
for beam = 1:Pcadp.nbeams, 
    eval(sprintf('Pcadp.v%1d = nandata;',beam)); 
    eval(sprintf('Pcadp.co%1d = nandata;',beam));
    eval(sprintf('Pcadp.sn%1d = nandata;',beam));
end
Pcadp.binhts = zeros(Pcadp.ppb, Pcadp.ncells);
Pcadp.instht = gmean(Pcadp.meanrange); % Pcadp.meanrange has already been corrected by ecorr
Pcadp.vx = ones(size(Pcadp.v1)).*NaN; % MATLAB's happier when you allocate large arrays ahead of time
Pcadp.vy = Pcadp.vx; Pcadp.vz = Pcadp.vx;
clear nandata
    
% deal with the CTD now, load it at once to avoid extra read in for loop
if ~isempty(cdfs{'CTD_temp'}),
    Pcadp.ctd(:,1) = cdfs{'CTD_temp'}(:);
    Pcadp.ctd(:,2) = cdfs{'CTD_cond'}(:);
    Pcadp.ctd(:,3) = cdfs{'CTD_press'}(:);
    Pcadp.ctd(:,4) = cdfs{'CTD_sal'}(:);
end


% ------------------------- main loop

ncburstidx = 1; % write index for output files
for cdfburstidx = bursts(1):bursts(2),
    % These were checked on 12/10/03 and look like I got the right thing in the right place
    % reading the contents of the PCrecordType data structure, what was performed in read_pc, see adpdata.h
    Pcadp.num  = (cdfb{'profile'}(cdfburstidx,:))';      % profile number
    % assumes all profiles are the same (to save space), so take the mean of these values
    % Lacy: Profile lag for resolution pulse (name should not include 'av', this is not an ambiguity velocity)
    Pcadp.av_plag = cdfs{'MeanResLag'}(cdfburstidx)./1000;    % ambiguity resolution pulse lag (m) converted from mm
    Pcadp.av = cdfs{'MeanResUa'}(cdfburstidx)./10;            % ambiguity velocity corresponding to the
    % resolution lag (velocity range) (cm/s) converted from mm/s
    Pcadp.cell_start = cdfs{'MeanResStart'}(cdfburstidx)./1000;       % beginning position of the resolution cell (m) converted from mm
    Pcadp.cell_length = cdfs{'MeanResLength'}(cdfburstidx)./1000;     % length of the resolution cell (m) converted from mm
    Pcadp.pr_plag = cdfs{'MeanPrfLag'}(cdfburstidx)./1000;    % profile pulse lag (m) converted from mm
    Pcadp.avpl = cdfs{'MeanPrfUa'}(cdfburstidx)./10;          % ambiguity velocity corresponding to the profile lag (cm/s) converted from mm/s
    Pcadp.prange_start = cdfs{'MeanPrfStart'}(cdfburstidx)./1000;     % beginning of the profile range (m) converted from mm
    Pcadp.prange_length = cdfs{'MeanPrfLength'}(cdfburstidx)./1000;   % length of the profile range (m) converted from mm
    for beam = 1:3, % beams 1 - 3
        varname = sprintf('Range%1d',beam);
        Pcadp.bot_dist(:,beam) = (cdfb{varname}(cdfburstidx,:)./1000)'; % range to bottom (m) converted from mm
        varname = sprintf('Ures%1d',beam);
        Pcadp.res_vel(:,beam) = (cdfb{varname}(cdfburstidx,:))'; % velocity from the resolution pulse in cm/s
        varname = sprintf('Cres%1d',beam);
        Pcadp.res_cor(:,beam) = (cdfb{varname}(cdfburstidx,:))'; % correlation from the resolution pulse (%)
    end
    
    % reading the contents of the ProfileHdrType_2 data structure, what was performed in read_hdr, see adpdata.h
    Pcadp.jt = (cdfb{'time'}(cdfburstidx,:)+cdfb{'time2'}(cdfburstidx,:)./(24*3600*1000))';
    Pcadp.nsamp = (cdfb{'Npings'}(cdfburstidx,:))';  % number of samples averaged for this profile
    Pcadp.c = cdfs{'MeanSoundspd'}(cdfburstidx); % sound velocity (m/s)
    Pcadp.hdg = (cdfb{'heading'}(cdfburstidx,:))';   % mean heading (degrees magnetic)
    Pcadp.pitch = (cdfb{'pitch'}(cdfburstidx,:))';   % rotation about Y axis, deg
    Pcadp.roll = (cdfb{'roll'}(cdfburstidx,:))';     % rotation about X axis, deg
    Pcadp.t = (cdfb{'temperature'}(cdfburstidx,:))'; % temperature, deg C
    Pcadp.pfreq = (cdfb{'extpressfreq'}(cdfburstidx,:))';  % pressure frequency, Hz
    Pcadp.extpress = (cdfb{'extpress'}(cdfburstidx,:))';    % pressure, counts, least significant bits
    P.MeanPress = (cdfb{'pressure'}(cdfburstidx,:))'; 
    Pcadp.hdg_sd = cdfb{'StdHeading'}(cdfburstidx); % sd mean heading (degrees magnetic)
    Pcadp.pitch_sd = cdfb{'StdPitch'}(cdfburstidx); % sd pitch, deg
    Pcadp.roll_sd = cdfb{'StdRoll'}(cdfburstidx);   % sd roll, deg
    Pcadp.t_sd = cdfs{'StdTemperature'}(cdfburstidx);   % sd temperature, deg C
    Pcadp.battV = (cdfs{'MeanBattery'}(cdfburstidx,:))';     % battery voltage
    Pcadp.ad1 = (cdfb{'extsensor1'}(cdfburstidx,:))';    % analog input ch 1, a/d volts
    Pcadp.ad2 = (cdfb{'extsensor2'}(cdfburstidx,:))';    % analog input ch 2, a/d volts
    Pcadp.ad1_sd = cdfs{'StdExtsensor1'}(cdfburstidx); % sd ch1, a/d volts
    Pcadp.ad2_sd = cdfs{'StdExtsensor2'}(cdfburstidx); % sd ch2, a/d volts

    for beam = 1:3, % beams 1 - 3
        % varname = sprintf('Range%1d',beam);
        % reading the contents of the Profile3BeamType data structure, what was performed in read_v1-3, see adpdata.h
        %disp(sprintf('Getting beam velocities for beam %1d',beam))
        eval(sprintf('Pcadp.v%1d = cdfb{''Vel%1d''}(cdfburstidx,:,:);',beam,beam)); % [profile, cell] in cm/s
        % reading the correlation from the PCrecordType data structure, what was performed in read_co1-3, see adpdata.h
        % GADPCOR extracts the variable VelStd, which is corellation in pulse coherent mode and STD in regular mode
        %disp(sprintf('Getting correlations for beam %1d',beam))
        eval(sprintf('Pcadp.cor%1d = cdfb{''Cor%1d''}(cdfburstidx,:,:);',beam,beam)); % [profile, cell] in %
        %disp(sprintf('Getting signal amplitudes for beam %1d',beam))
        eval(sprintf('Pcadp.sn%1d = cdfb{''Amp%1d''}(cdfburstidx,:,:);',beam,beam)); % [profile, cell] in counts
    end
    
    % find number of complete bursts and indices of first profiles in each burst
    % make this work for Lacy for mat with 1 burst
    Pcadp.nburst = 1; % not to be confused with nbursts in the file given by Pcadp.nbursts
    Pcadp.bur_ind = 1:Pcadp.ppb; %(1:Pcadp.ppb:(Pcadp.nburst-1)*Pcadp.ppb+1);
    
    % calculate height above bottom of cells
    % note: negative values are below bed
    %b.cellsize = b.cellsize/10; %original Lacy code
    for icell = 1:Pcadp.ncells;
        % as defined by Sontek
        Pcadp.binhts(:,icell) = Pcadp.instht - Pcadp.blank_dist - (icell)*Pcadp.cellsize;
    end
    
    %if cdfb{'burst'}(cdfburstidx) > 121, keyboard; end
    % resolve ambiguity before applying rotations
    if switches.resolve_ambiguity,
        if isnan(Pcadp.cell_start),
            fprintf('Skipping ambiguity resolution of burst %d because of NaN in stats file\n',...
                cdfb{'burst'}(cdfburstidx))
        else
            [rv1c,rv2c,rv3c, v1c, v2c,v3c, bc1, bc2, bc3] = pcadpambiguity(Pcadp, switches.amb);
            % get nout1-3 back for qa info
            Pcadp.v1 = v1c;
            Pcadp.v2 = v2c;
            Pcadp.v3 = v3c;
        end
    end

    % ----------------------- Do Rotations
    % going to do it how Sontek does it, using the stored matrix
    % the stored matrix is simply the sines and cosines, but in case a beam misaligment must be accounted for, it's taken care of here.
    % in adpdata.c ReadAdpConf calls ComputeBeamToXyzMat in adpdata1.c
    % get the conversion matrix for 123 to XYZ using technique ComputeBeamToXyzMat in adpdata1.c
    Pcadp.nbeams = cdfb.PCADPProbeNbeams(:);
    Pcadp.beamtoxyzmat = cdfb.PCADPProbeXformMat.*(1.0/1024.0); % this is a 1x9 array of [x x x y y y z z z] transforms
    Pcadp.beamtoxyzmat = round(Pcadp.beamtoxyzmat.*1000)./1000; % there are really only 3 significant digits, MATLAB adds more
    Pcadp.beamtoxyzmat = reshape(Pcadp.beamtoxyzmat,Pcadp.nbeams,Pcadp.nbeams)'; % [x x x; y y y; z z z], i.e. [axis][beam]
    Pcadp.orientation = cdfb.SensorOrientation(:);
    if strcmp(Pcadp.orientation,'down'),
        Pcadp.beamtoxyzmat(2,:) = Pcadp.beamtoxyzmat(2,:).*(-1); % reverse Y axis
        Pcadp.beamtoxyzmat(3,:) = Pcadp.beamtoxyzmat(3,:).*(-1); % reverse Z axis
    end
    % do it the MATLAB way, adapted from ConvertVelPrfToOutCoord in adpdata1.c
    % results checked agains GADPVEL -CXYZ output 12/16/03
    % get the data into the right shape for a Matlab Matrix Multiplication
    % this is the beam alignment matrix
    %Pcadp.vx = ones(size(Pcadp.v1)).*NaN; % MATLAB's happier when you allocate large arrays ahead of time
    %Pcadp.vy = Pcadp.vx; Pcadp.vz = Pcadp.vx;
    % TODO can we vectorize this for loop?
    for profile = 1:Pcadp.ppb,
        PrfVel123 = [Pcadp.v1(profile,:); Pcadp.v2(profile,:); Pcadp.v3(profile,:);]; % axis by cell
        PrfVelXYZ = Pcadp.beamtoxyzmat*PrfVel123;
        Pcadp.vx(profile,:) = PrfVelXYZ(1,:);
        Pcadp.vy(profile,:) = PrfVelXYZ(2,:);
        Pcadp.vz(profile,:) = PrfVelXYZ(3,:);
        ResVel123 = Pcadp.res_vel(profile,:)'; % axis
        ResVelXYZ(profile,:) = Pcadp.beamtoxyzmat*ResVel123; % profile by beam
    end
    % override orientation if user wants to
    if isfield(switches,'override_heading') && ~isempty(switches.override_heading),
        Pcadp.hdg = switches.override_heading(cdfburstidx);
    end
    if isfield(switches,'override_pitch') && ~isempty(switches.override_pitch),
        Pcadp.pitch = switches.override_pitch(cdfburstidx);
    end
    if isfield(switches,'override_roll') && ~isempty(switches.override_roll),
        Pcadp.roll = switches.override_roll(cdfburstidx);
    end
    % correct for declination
    Pcadp.hdg = Pcadp.hdg + Pcadp.magnetic_variation;
    idx = find(Pcadp.hdg>=360.);
    if ~isempty(idx), Pcadp.hdg(idx)=Pcadp.hdg(idx)-360.; end
    idx = find(Pcadp.hdg<0);
    if ~isempty(idx), Pcadp.hdg(idx)=Pcadp.hdg(idx)+360.; end
    % Jessie's way, apply the sines and cosines as such.  Probably faster.
    %[Pcadp.vx, Pcadp.vy, Pcadp.vz] = bm2xyz(Pcadp.v1,Pcadp.v2,Pcadp.v3);
    % note that xyz2enu does NOT check for noise in the orientation data, if given a vector, it chooses the median value
    [Pcadp.ve, Pcadp.vn, Pcadp.vu] = xyz2enu(Pcadp.vx,Pcadp.vy,Pcadp.vz,Pcadp.hdg,Pcadp.pitch,Pcadp.roll,Pcadp.ppb);
    [res.ve, res.vn, res.vu] = xyz2enu(ResVelXYZ(:,1),ResVelXYZ(:,2),ResVelXYZ(:,3),Pcadp.hdg,Pcadp.pitch,Pcadp.roll,Pcadp.ppb);

    % calculate horizontal speed and direction
    Pcadp.hspd = zeros(size(Pcadp.ve));
    for icell = 1:Pcadp.ncells,
        % replace pcoord with uv2polar which comes with the toolbox
        %[Pcadp.hspd(:,icell), Pcadp.hdir(:,icell)] = pcoord(Pcadp.ve(:,icell),Pcadp.vn(:,icell));
        [Pcadp.hdir(:,icell),Pcadp.hspd(:,icell)] = uv2polar(Pcadp.ve(:,icell),Pcadp.vn(:,icell));
    end

    %   switches.range_correct = 1; adjust distance to bottom measured by beam
    %if switches.range_correct,
    %    Pcadp.bot_dist = vert_bd(Pcadp.bot_dist, Pcadp.pitch, Pcadp.roll);
    %end

    % process pressure
    % build the pressure data structure needed
    P.temp = Pcadp.t; % temperature, deg C
    P.extpress = Pcadp.extpress;    % pressure, counts, least significant bits
    P.stdpress = Pcadp.p_sd; % std of pressure, for strain guage or smart paros
    P.extpressfreq = Pcadp.pfreq;
    P = hydrapress2db(P);
    Pcadp.p = P.pressdb;
    
    % process the external sensors
    for i=1:2,
        varname = sprintf('extsensor%d',i);
        cdfobj = cdfb{varname};
        if ~isempty(cdfobj)
            eval(sprintf('optic{i} = dohydraoptic(extmeta{i},Pcadp.ad%d);',i));
            if isempty(optic{i}), optic{i} = zeros([Pcadp.ppb, 1]); end
        end
    end
    
    % write the burst data as shape [profile, time, depth, lat, lon] | [time, depth, lat, lon]
    ncb{'burst'}(ncburstidx) = cdfb{'burst'}(cdfburstidx);
    ncb{'profile'}(ncburstidx,:) = cdfb{'profile'}(cdfburstidx,:);
    ncb{'time'}(ncburstidx,:) = cdfb{'time'}(cdfburstidx,:);
    ncb{'time2'}(ncburstidx,:) = cdfb{'time2'}(cdfburstidx,:);
    ncb{'u_1205'}(ncburstidx,:) = Pcadp.ve;
    ncb{'v_1206'}(ncburstidx,:) = Pcadp.vn;
    ncb{'w_1204'}(ncburstidx,:) = Pcadp.vu;
    ncb{'ResU'}(ncburstidx,:) = res.ve;
    ncb{'ResV'}(ncburstidx,:) = res.vn;
    ncb{'ResW'}(ncburstidx,:) = res.vu;
    ncb{'AGC1_1221'}(ncburstidx,:) = Pcadp.sn1;
    ncb{'AGC2_1222'}(ncburstidx,:) = Pcadp.sn2;
    ncb{'AGC3_1223'}(ncburstidx,:) = Pcadp.sn3;
    % add corrs here
    ncb{'cor1_1285'}(ncburstidx,:) = Pcadp.cor1;
    ncb{'cor2_1286'}(ncburstidx,:) = Pcadp.cor2;
    ncb{'cor3_1287'}(ncburstidx,:) = Pcadp.cor3;
    
    if ~isempty(ncb{'P_4023'}),
        ncb{'P_4023'}(ncburstidx,:) = (Pcadp.p).*100; % convert from db to mbar
    end
    ncb{'Tx_1211'}(ncburstidx,:) = Pcadp.t;
    
    % write the external sensors for both burst and stats
    for i=1:2,
        cdfobj = cdfb{sprintf('extsensor%d',i)};
        if ~isempty(cdfobj)
            if ~isempty(findstr('TRANS',cdfobj.serial(:))),
                ncvarname = sprintf('ATTN%d_55',i);
                eval(sprintf('ncb{''tran%d_4010''}(ncburstidx,:) = Pcadp.ad%d;',i,i))
                eval(sprintf('ncs{''tran%d_4010''}(ncburstidx) = gmean(Pcadp.ad%d);',i,i))
            elseif ~isempty(findstr('OBS',cdfobj.serial(:))),
                ncvarname = sprintf('Sed%d_981',i);
                eval(sprintf('ncb{''NEP%d_56''}(ncburstidx,:) = Pcadp.ad%d;',i,i))
                eval(sprintf('ncs{''NEP%d_56''}(ncburstidx) = gmean(Pcadp.ad%d);',i,i))
            else
                ncvarname = sprintf('extsensor%d',i);
            end
            if ~isempty(ncb{ncvarname}),
                ncb{ncvarname}(ncburstidx,:) = optic{i};
                ncs{ncvarname}(ncburstidx) = gmean(optic{i});
            end
        end
    end
    
    % write the stats data
    ncs{'burst'}(ncburstidx) = cdfb{'burst'}(cdfburstidx);
    % time for the stats is the beginning of the burst
    ncs{'time'}(ncburstidx) = cdfb{'time'}(cdfburstidx,1); 
    ncs{'time2'}(ncburstidx) = cdfb{'time2'}(cdfburstidx,1);
    uMean = gmean(Pcadp.ve);
    ncs{'u_1205'}(ncburstidx,:,1,1) = uMean;
    ncs{'USTD_4097'}(ncburstidx,:,1,1) = gstd(Pcadp.ve);
    ncs{'u_1205min'}(ncburstidx,:,1,1) = gmin(Pcadp.ve);
    ncs{'u_1205max'}(ncburstidx,:,1,1) = gmax(Pcadp.ve);
    vMean = gmean(Pcadp.vn);
    ncs{'v_1206'}(ncburstidx,:,1,1) = vMean;
    ncs{'VSTD_4098'}(ncburstidx,:,1,1) = gstd(Pcadp.vn);
    ncs{'v_1206min'}(ncburstidx,:,1,1) = gmin(Pcadp.vn);
    ncs{'v_1206max'}(ncburstidx,:,1,1) = gmax(Pcadp.vn);
    ncs{'w_1204'}(ncburstidx,:,1,1) = gmean(Pcadp.vu);
    ncs{'WSTD_4099'}(ncburstidx,:,1,1) = gstd(Pcadp.vu);
    ncs{'w_1204min'}(ncburstidx,:,1,1) = gmin(Pcadp.vu);
    ncs{'w_1204max'}(ncburstidx,:,1,1) = gmax(Pcadp.vu);
    [dMean, sMean] = uv2polar(uMean,vMean);
    ncs{'CS_300'}(ncburstidx,:,1,1) = sMean;
    ncs{'CD_310'}(ncburstidx,:,1,1) = dMean;
    ncs{'AGC1_1221'}(ncburstidx,:,1,1) = gmean(Pcadp.sn1);
    ncs{'AGC2_1222'}(ncburstidx,:,1,1) = gmean(Pcadp.sn2);
    ncs{'AGC3_1223'}(ncburstidx,:,1,1) = gmean(Pcadp.sn3);
    % add corrs here
    ncs{'cor1_1285'}(ncburstidx,:) = gmean(Pcadp.cor1);
    ncs{'cor2_1286'}(ncburstidx,:) = gmean(Pcadp.cor2);
    ncs{'cor3_1287'}(ncburstidx,:) = gmean(Pcadp.cor3);
    ncs{'Hdg_1215'}(ncburstidx,1,1,1) = gmedian(Pcadp.hdg);
    ncs{'HSD_1218'}(ncburstidx,1,1,1) = gstd(Pcadp.hdg);
    ncs{'Ptch_1216'}(ncburstidx,1,1,1) = gmedian(Pcadp.pitch);
    ncs{'PSD_1219'}(ncburstidx,1,1,1) = gstd(Pcadp.pitch);
    ncs{'Roll_1217'}(ncburstidx,1,1,1) = gmedian(Pcadp.roll);
    ncs{'RSD_1220'}(ncburstidx,1,1,1) = gstd(Pcadp.roll);
    if ~isempty(ncs{'P_4023'}), 
        ncs{'P_4023'}(ncburstidx,1,1,1) = gmean(Pcadp.p).*100; % change from db to mbar
        ncs{'SDP_850'}(ncburstidx,1,1,1) = gstd(Pcadp.p).*100; % should be in mbar
        Pcadp.p_sd = gstd(Pcadp.p);  % sd pressure, dBars
    end
    ncs{'Tx_1211'}(ncburstidx,1,1,1) = gmean(Pcadp.t);    
   
    if cdfburstidx<50 && ~rem(cdfburstidx,10), fprintf('Finished burst %d, %f mintues elapsed\n',cdfb{'burst'}(cdfburstidx),toc./60), end
    if cdfburstidx>50 && ~rem(cdfburstidx,100), fprintf('Finished burst %d, %f mintues elapsed\n',cdfb{'burst'}(cdfburstidx),toc./60), end

    ncburstidx = ncburstidx +1;
end % end of for burstidx = 1:Pcadp.nbursts

% ------------------------- wrap up

% rotate and write the ambiguity (resolution bin) velocities now
for beam = 1:3,
    varname = sprintf('MeanUres%1d',beam);
    res.vel(:,beam) = cdfs{varname}(bursts(1):bursts(2)); % velocity from the resolution pulse (cm/s) 
end
for profile = 1:length(res.vel),
    PrfVel123 = res.vel(profile,:)'; % axis by one cell
    PrfVelXYZ = Pcadp.beamtoxyzmat*PrfVel123;
    res.vx(profile,:) = PrfVelXYZ(1,:);
    res.vy(profile,:) = PrfVelXYZ(2,:);
    res.vz(profile,:) = PrfVelXYZ(3,:);
end
% override orientation if user wants to
if isfield(switches,'override_heading') && ~isempty(switches.override_heading),
    res.hdg = gmean(switches.override_heading);
else res.hdg = cdfs{'MedianHeading'}(bursts(1):bursts(2));
end
if isfield(switches,'override_pitch') && ~isempty(switches.override_pitch),
    res.pitch = gmean(switches.override_pitch);
else res.pitch = cdfs{'MedianPitch'}(bursts(1):bursts(2));
end
if isfield(switches,'override_roll') && ~isempty(switches.override_roll),
    res.roll = gmean(switches.override_roll);
else res.roll = cdfs{'MedianRoll'}(bursts(1):bursts(2));
end
% correct for declination
res.hdg = res.hdg + Pcadp.magnetic_variation;
idx = find(res.hdg>=360.);
if ~isempty(idx), res.hdg(idx)=res.hdg(idx)-360.; end
% Jessie's way, apply the sines and cosines as such.  Probably faster.
%[Pcadp.vx, Pcadp.vy, Pcadp.vz] = bm2xyz(Pcadp.v1,Pcadp.v2,Pcadp.v3);
% note that xyz2enu does NOT check for noise in the orientation data, if given a vector, it chooses the median value
[res.ve, res.vn, res.vu] = xyz2enu(res.vx,res.vy,res.vz,res.hdg,res.pitch,res.roll,length(res.vel));
ncs{'ResU'}(:) = res.ve;
ncs{'ResV'}(:) = res.vn;
ncs{'ResW'}(:) = res.vu;

% put in the bin depth based on WATER_DEPTH, instrument height and
% configuration- 
% per sontek, you do not add a half cell size to get distance to
% measurement- it's just (blanking_distance + bin*cell_size).  This works
% because the measurement volume is actually 2*cell_size, and there's a 25%
% overlap with each adjacent cell.
%
% anyhow- don't overthink this and try to add back the .5 cell_size-
% they've already factored it in.
%
% previous methods
% r1082
% ncb{'depth'}(icell) = cdfb.WATER_DEPTH(:) - Pcadp.instht + Pcadp.blank_dist + (icell).*Pcadp.cellsize;
% r1649
% ncb{'depth'}(icell) = cdfb.WATER_DEPTH(:) - Pcadp.instht + Pcadp.blank_dist + (Pcadp.cellsize/2) + (icell-1).*Pcadp.cellsize;
% r1714 - note the difference in parentheses placement- it's critical!!!
% ncb{'depth'}(icell) = cdfb.WATER_DEPTH(:) - (Pcadp.instht + Pcadp.blank_dist + icell.*Pcadp.cellsize);
%
% note the difference in parentheses placement- it's critical!!!
for icell = 1:Pcadp.ncells;
    ncs{'depth'}(icell) = cdfb.WATER_DEPTH(:) - Pcadp.instht + Pcadp.blank_dist + (icell.*Pcadp.cellsize);
    ncb{'depth'}(icell) = cdfb.WATER_DEPTH(:) - Pcadp.instht + Pcadp.blank_dist + (icell.*Pcadp.cellsize);
end

% add bindist to show the transducer-centric view
bindist= Pcadp.blank_dist+Pcadp.cellsize:Pcadp.cellsize:cdfb.PCADPProbeProfRange;
 % trim if it ends up longer than needed
 if length(bindist) > Pcadp.ncells
     bindist=(bindist(1:Pcadp.ncells));
 end
 ncs{'bindist'}(1:length(bindist))=bindist;
 ncb{'bindist'}(1:length(bindist))=bindist;

%write the range for the three beams as [time, beam, lat, lon]
if strcmp(cdfb.PCADPUserSetupPCMode(:),'Yes'),
    ncs{'brange'}(:,:,1,1) = Pcadp.range; % range to bottom (m)
end

% write the CTD all at once
if ~isempty(cdfs{'CTD_temp'}),
    ncs{'CTDTMP_4211'}(:,1,1,1) = Pcadp.ctd(:,1);
    ncs{'CTDCON_4218'}(:,1,1,1) = Pcadp.ctd(:,2);
    if ~isempty(ncs{'CTDPRS_4203'}), 
        ncs{'CTDPRS_4203'}(:,1,1,1) = Pcadp.ctd(:,3)-10.1325;
    end
    ncs{'CTDSAL_4214'}(:,1,1,1) = Pcadp.ctd(:,4);
end    

% add to the history
ncs.history = sprintf('velocities rotated by %s; %s: %s', mfilename, mversion,cdfs.history(:));
ncb.history = sprintf('velocities rotated by %s; %s: %s',mfilename, mversion,cdfb.history(:));

if switches.resolve_ambiguity,
    ncs.history = sprintf('ambiguity resolution applied, %s',ncs.history(:));
    ncb.history = sprintf('ambiguity resolution applied, %s',ncb.history(:));
end

% a couple more EPIC attributes
tjul = ncs{'time'}(1)+ncs{'time2'}(1)./(3600*1000*24);
ncs.start_time = ncchar(datestr(datenum(gregorian(tjul))));
tjul = ncs{'time'}(end)+ncs{'time2'}(end)./(3600*1000*24);
ncs.stop_time = ncchar(datestr(datenum(gregorian(tjul))));
tjul = ncb{'time'}(1,1)+ncb{'time2'}(1,1)./(3600*1000*24);
ncb.start_time = ncchar(datestr(datenum(gregorian(tjul))));
tjul = ncb{'time'}(end,end)+ncb{'time2'}(end,end)./(3600*1000*24);
ncb.stop_time = ncchar(datestr(datenum(gregorian(tjul))));

ncs=add_minmaxvalues(ncs);
%  Take this outside the script and use add_minmax2any
%add_minmaxvalues(ncb);
remove_PCADPglobal_atts(ncs)
remove_PCADPglobal_atts(ncb)

close(ncs)
close(ncb)
close(cdfb)
close(cdfs)

fprintf('Finished writing data, %4.2f mintues elapsed\n',toc./60)

return

% ------------------- sub-functions

function nc = definencbfile(ncFile,cdfb, bursts)

disp('Now writing the new netCDF burst file ... this may take a while')

% the new burst file
[pathstr,filename,fileext,versn] = fileparts(name(cdfb));
ncFile = fullfile(pathstr,[ncFile 'b-cal.nc' versn]);
nc = netcdf(ncFile, 'clobber');

% copy all the global attributes for now.  Adjust this later by doing it by
% name as the variables are done below
ncobjs = att(cdfb);
for n = 1:length(ncobjs), copy(ncobjs{n},nc); end

nc.DESCRIPT = ncchar('Sontek PCADP calibrated data burst file');
nc.CREATION_DATE = datestr(now);

% new dimensions, going from profile number to time
nc('time') = 0; % the record dimension, set fixed as we are in define mode
nc('depth') = cdfb.Ncells(:);
nc('lon') = 1;
nc('lat') = 1;
nc('profile') = cdfb.ProfilesPerBurst(:);

% copy some of the variables, do the 1D ones first
% basic indeces
nc{'burst'} = ncfloat('time'); 
nc{'burst'}.units = ncchar('count');
nc{'burst'}.name = ncchar('burst');
nc{'burst'}.generic_name = ncchar('record');
nc{'burst'}.epic_code = nclong(1207);
nc{'burst'}.long_name = ncchar('Burst Number');

nc{'profile'} = ncfloat('time','profile'); 
nc{'profile'}.units = ncchar('count');
nc{'profile'}.generic_name = ncchar('record');
nc{'profile'}.epic_code = nclong(1207);
nc{'profile'}.long_name = ncchar('Profile Number');

nc{'time'} = nclong('time','profile'); 
nc{'time'}.FORTRAN_format = ncchar('F10.2');
nc{'time'}.units = ncchar('True Julian Day');
nc{'time'}.type = ncchar('EVEN');
nc{'time'}.name = ncchar('time');
nc{'time'}.generic_name = ncchar('time');
nc{'time'}.epic_code = nclong(624);
nc{'time'}.comment = ncchar('UT Julian days that begin at midnight; 1968-05-23 = 2440000');%old: nc{'time'}.comment = ncchar('UT Julian days using USGS convention where 1968-05-23 00:00:00 UT = 2440000'); FLL

nc{'time2'} = nclong('time','profile'); 
nc{'time2'}.FORTRAN_format = ncchar('F10.2');
nc{'time2'}.type = ncchar('EVEN');
nc{'time2'}.name = ncchar('time2');
nc{'time2'}.generic_name = ncchar('time');
nc{'time2'}.epic_code = nclong(624);
nc{'time2'}.units = ncchar('msec since 0:00 GMT');

nc{'depth'} = ncfloat('depth'); 
nc{'depth'}.FORTRAN_format = ncchar('F10.2');
nc{'depth'}.units = ncchar('m');
nc{'depth'}.type = ncchar('EVEN');
nc{'depth'}.name = ncchar('D');
nc{'depth'}.long_name = ncchar('DEPTH (m)');
nc{'depth'}.generic_name = ncchar('depth');
nc{'depth'}.epic_code = nclong(3);
nc{'depth'}.blanking_distance = ncdouble(cdfb.PCADPProbeBlankDistance(:));
nc{'depth'}.bin_size = ncdouble(cdfb.CellSize(:));
nc{'depth'}.transducer_offset_from_bottom = ncdouble(cdfb.PCADPProbeHeight(:));

nc{'lon'} = ncfloat('lon'); 
nc{'lon'}.FORTRAN_format = ncchar('f10.4');
nc{'lon'}.units = ncchar('degree_east');
nc{'lon'}.type = ncchar('EVEN');
nc{'lon'}.epic_code = nclong(502);
nc{'lon'}.name = ncchar('LON');
nc{'lon'}.long_name = ncchar('LONGITUDE');
nc{'lon'}.generic_name = ncchar('lon');

nc{'lat'} = ncfloat('lat'); 
nc{'lat'}.FORTRAN_format = ncchar('F10.2');
nc{'lat'}.units = ncchar('degree_north');
nc{'lat'}.type = ncchar('EVEN');
nc{'lat'}.epic_code = nclong(500);
nc{'lat'}.name = ncchar('LAT');
nc{'lat'}.long_name = ncchar('LATITUDE');
nc{'lat'}.generic_name = ncchar('lat');

names  = {'u','v','w'};
long_names = {'Eastward','Northward','Vertical'};
epcodes = [1205 1206 1204];
% write the new 2D rotated velocities
for i = 1:length(names),
    ncname = sprintf('%s_%d',names{i},epcodes(i));    
    nc{ncname} = ncfloat('time', 'profile', 'depth'); ncobj = nc{ncname};
    ncobj.name = ncchar(names{i});
    ncobj.long_name = ncchar([long_names{i} ' Velocity']);
    ncobj.generic_name = ncchar(names{i});
    ncobj.units = ncchar('cm/s');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-1000 1000]);
end

names  = {'ResU','ResV','ResW'};
long_names = {'Eastward','Northward','Vertical'};
%epcodes = [0 0 0];
% write the new 2D rotated velocities
for i = 1:length(names),
    ncname = names{i};    
    nc{ncname} = ncfloat('time', 'profile'); ncobj = nc{ncname};
    ncobj.name = ncchar(names{i});
    ncobj.long_name = ncchar([long_names{i} ' Resolution Velocity']);
    ncobj.generic_name = ncchar(names{i});
    ncobj.units = ncchar('cm/s');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-1000 1000]);
end

% include the signal amplitude of each axis
% 1221:AGC1:Echo Intensity (AGC) Beam 1:AGC:counts: :ADCP Beam 1 AGC
% 1222:AGC2:Echo Intensity (AGC) Beam 2:AGC:counts: :ADCP Beam 2 AGC
% 1223:AGC3:Echo Intensity (AGC) Beam 3:AGC:counts: :ADCP Beam 3 AGC
% Signal strength, recorded for each ADV receiver, is a measure of the intensity of the reflected
% acoustic signal. This is recorded as raw signal amplitude in internal logarithmic units of counts;
% one count is equal to 0.43 dB. With the ADV software, signal strength can be accessed either as
% signal amplitude in counts or as signal-to-noise ratio (SNR) in dB. SNR is derived from signal
% amplitude by subtracting the ambient noise level and converting to units of dB.
epcodes = [1221 1222 1223];
for i = 1:length(epcodes),
    ncname = sprintf('AGC%1d_%d',i,epcodes(i));    
    nc{ncname} = ncfloat('time', 'profile', 'depth'); ncobj = nc{ncname};
    ncobj.name = ncchar(sprintf('AGC%1d',i));
    ncobj.long_name = ncchar(sprintf('Echo Intensity (AGC) Beam %1d',i));
    ncobj.generic_name = ncchar('AGC');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.units = ncchar('counts');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.valid_range = [0 1000]; % max range not known
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
end
% include the signal correlation of each axis
% 1285:cor1:Beam 1 correlation: :counts:i10:ADCP correlation of beam 1
% 1286:cor2:Beam 2 correlation: :counts:i10:ADCP correlation of beam 2
% 1287:cor3:Beam 3 correlation: :counts:i10:ADCP correlation of beam 3
epcodes = [1285 1286 1287];
for i = 1:length(epcodes),
    ncname = sprintf('cor%1d_%d',i,epcodes(i));    
    nc{ncname} = ncfloat('time', 'profile', 'depth'); ncobj = nc{ncname};
    ncobj.name = ncchar(sprintf('cor%1d',i));
    ncobj.long_name = ncchar(sprintf(' Beam %1d Correlation',i));
    ncobj.generic_name = ncchar('Cor');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.units = ncchar('percent');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = cdfb.PCADPProbeSerial(:);
    ncobj.valid_range = [0 100]; % max range not known
end

% modelled loosely after this key
%4023:P   :INTERVAL PRESSURE        :pres:mbar:f10.3:interval pressure measurement, Seadata tripod
if ~isempty(cdfb{'extpressfreq'}) || ~isempty(cdfb{'extpress'}) || ~isempty(cdfb{'pressure'}),
    nc{'P_4023'} = ncfloat('time', 'profile');  ncobj = nc{'P_4023'};
    ncobj.name = ncchar('P');
    ncobj.long_name = ncchar('PRESSURE            ');
    ncobj.epic_code = nclong(4023);
    ncobj.generic_name = ncchar('pres');
    ncobj.FORTRAN_format = ncchar('f10.3');
    ncobj.units = ncchar('mbar');
    if strcmp(cdfb.ExtPressInstalled(:),'No'),
        ncobj.sensor_type = ncchar('Sontek PCADP Internal');
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
        ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
        ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    else
        ncobj.sensor_type = ncchar(cdfb.ExtPressInstalled(:));
    end
    if strcmp(cdfb.ExtPressInstalled(:),'Druck') || strcmp(cdfb.ExtPressInstalled(:),'ParosFreq'),
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb{'extpressfreq'}.height(:));
        ncobj.initial_sensor_height = ncdouble(cdfb{'extpressfreq'}.height(:));
        ncobj.serial = cdfb{'extpressfreq'}.serial(:);
    % doesn't seem to get here- had some internal pressure that eneded up with no
    % depth and height attributes...
    elseif strcmp(cdfb.PressInstalled,'Yes'), % internal pressure
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
        ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
        ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    end
    %ncobj.serial = ncchar(cdfb.pressureSerialNum(:));
    ncobj.valid_range = ncfloat([0 100000]);
end

% 1211:Tx  :ADCP Transducer Temp.    :temp:deg. C:F10.2:ADCP Transducer Temp
nc{'Tx_1211'} = ncfloat('time','profile'); ncobj = nc{'Tx_1211'};
ncobj.name = ncchar('Tx');
ncobj.long_name = ncchar('Transducer Temp.    ');
ncobj.generic_name = ncchar('temp');
ncobj.epic_code = nclong(1211);
ncobj.units = ncchar('deg. C');
ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
ncobj.valid_range = ncfloat([-5 40]);

% bindist  :distance to bin center from transducer head.  not an EPIC variable
nc{'bindist'} = ncfloat('depth'); ncobj = nc{'bindist'};
ncobj.name = ncchar('bindist');
ncobj.long_name = ncchar('Distance from transducer head to bin center');
ncobj.generic_name = ncchar('bindist');
ncobj.units = ncchar('m');
ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
ncobj.valid_range = ncfloat([0 5]);

for i=1:2,
    cdfobj = cdfb{sprintf('extsensor%d',i)};
    if ~isempty(cdfobj),
        if findstr('TRANS',cdfobj.serial(:)), % a transmissometer
            %  55:ATTN:ATTENUATION             :attn:m-1:f7.5:added for r2d2 ctd data
            ncvarname = sprintf('ATTN%d_55',i);
            nc{ncvarname} = ncfloat('time','profile'); ncobj = nc{ncvarname};
            transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('attn%d',i));
            ncobj.long_name = ncchar(sprintf('ATTENUATION #%d',i));
            ncobj.generic_name = ncchar('attn');
            ncobj.epic_code = nclong(55);
            ncobj.units = ncchar('m-1');
            ncobj.sensor_type = ncchar('transmissometer');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.FORTRAN_format = ncchar('F7.5');
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 10000]);
            ncobj.varcomment.ATTN_55 = 'ATTN = -(1/focal_length) .* log(tran_4010./(.95*Vair))';
            % 4010:tran:TRANSMISSION (VOLTS)     :trans:volts:f10.3:basic measurement, transmissometer
            ncvarname = sprintf('tran%d_4010',i);
            nc{ncvarname} = ncfloat('time','profile'); ncobj = nc{ncvarname};
            transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('tran%d',i));
            ncobj.long_name = ncchar(sprintf('TRANSMISSION (VOLTS) #%d',i));
            ncobj.generic_name = ncchar('transmissometer');
            ncobj.epic_code = nclong(4010);
            ncobj.units = ncchar('volts');
            ncobj.FORTRAN_format = ncchar('F10.3');
            ncobj.sensor_type = ncchar('trans');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 5]);
        elseif findstr('OBS',cdfobj.serial(:)), % an OBS
            ncvarname = sprintf('NEP%d_56',i);
            %  56:NEP:BACKSCATTER INTENSITY    :nephelometer:v:f10.6: added for Dave Pashinski
            nc{ncvarname} = ncfloat('time','profile'); ncobj = nc{ncvarname};
            ncals = transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('nep%d',i));
            ncobj.long_name = ncchar(sprintf('BACKSCATTER INTENSITY #%d',i));
            ncobj.generic_name = ncchar('nephelometer');
            ncobj.epic_code = nclong(56);
            ncobj.units = ncchar('v');
            ncobj.sensor_type = ncchar('nephelometer');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.FORTRAN_format = ncchar('F10.6');
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 5]);
            % now, if there's cal info to compute concentration, then set up a variable
            if (ncals > 0) && ~isempty(cdfobj.cal_coef(:)),
            % 981:Sed:Sediment concentration   :sed:g/l:f10.2:sediment concentration
                ncvarname = sprintf('Sed%d_981',i);
                nc{ncvarname} = ncfloat('time','profile'); ncobj = nc{ncvarname};
                transfer_calvalues(cdfobj,ncobj);
                ncobj.name = ncchar(sprintf('sed%d',i));
                ncobj.long_name = ncchar(sprintf('Sediment concentration %d',i));
                ncobj.generic_name = ncchar('sed');
                ncobj.epic_code = nclong(981);
                ncobj.units = ncchar('g/l');
                ncobj.sensor_type = ncchar('nephelometer');
                ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
                ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
                ncobj.serial = ncchar(cdfobj.serial(:));
                ncobj.minimum = ncfloat(0);
                ncobj.maximum = ncfloat(0);
                ncobj.valid_range = ncfloat([0 10000]);
            else
                fprintf('No coefficients for ext sensor %d, Sed variable not written\n',i)
            end
        else
            ncvarname = sprintf('extsensor%d',i);
            nc{ncvarname} = ncfloat('time','profile'); ncobj = nc{ncvarname};
            transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('ext%d',i));
            ncobj.long_name = ncchar(sprintf('external sensor #%d',i));
            ncobj.generic_name = ncchar(' ');
            ncobj.epic_code = nclong(0);
            ncobj.units = ncchar('v');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 5]);
        end
    end
end

% speed and direction are omitted here to save disk space
% heading, pitch, roll, cor and amp are also omitted.  All of these can be
% found in the .cdf raw file

add_fillvalues(nc,1e35)
add_vardesc(nc)

return

function nc = definencsfile(ncFile, cdfb, cdfs)

fprintf('Defining the new netCDF average file ... %d mintues elapsed so far\n',toc./60)

% the new average file
[pathstr,filename,fileext,versn] = fileparts(name(cdfb));
ncfile = fullfile(pathstr,[ncFile 's-cal.nc' versn]);
nc = netcdf(ncfile, 'clobber');

% copy all the global attributes for now.  Adjust this later by doing it by
% name as the variables are done below
ncobjs = att(cdfb);
for n = 1:length(ncobjs), 
    copy(ncobjs{n},nc); 
end

nc.DESCRIPT = ncchar('Sontek PCADP calibrated data statistics file');
nc.CREATION_DATE = datestr(now);

% new dimensions, going from profile number to time
nc('time') = 0; % the record dimension, set fixed as we are in define mode
nc('depth') = cdfb.Ncells(:);
nc('lon') = 1;
nc('lat') = 1;
nc('beam') = 3; % for the three beam ranges

% define first, then load, as we must do statistics, etc.

% basic indeces
nc{'burst'} = ncfloat('time'); 
nc{'burst'}.units = ncchar('count');
nc{'burst'}.name = ncchar('burst');
nc{'burst'}.long_name = ncchar('Burst Number');

nc{'time'} = nclong('time'); ncobj = nc{'time'};
ncobj.units = ncchar('Julian Day');
ncobj.FORTRAN_format = ncchar('F10.2');
ncobj.long_name = ncchar('EPIC SYSTEM TIME');
ncobj.epic_code = nclong(624);
ncobj.type = ncchar('EVEN');
ncobj.valid_min = nclong(0);
ncobj.comment = ncchar('UT Julian days that begin at midnight; 1968-05-23 = 2440000');%old: nc{'time'}.comment = ncchar('UT Julian days using USGS convention where 1968-05-23 00:00:00 UT = 2440000'); FLL
ncobj.comment1 = ncchar('time, taken from the instrument, fixed at beginning of burst');
% copy all the global attributes for now.  Adjust this later by doing it by
% name as the variables are done below
ncobjs = att(cdfb);
for n = 1:length(ncobjs), 
    copy(ncobjs{n},nc); 
end

nc{'time2'} = nclong('time');  ncobj = nc{'time2'};
ncobj.FORTRAN_format = ncchar('F10.2');
ncobj.units = ncchar('msec');
ncobj.type = ncchar('EVEN');
ncobj.long_name = ncchar('msec since 0:00 GMT');
ncobj.epic_code = nclong(624);
ncobj.valid_min = nclong(0);

nc{'depth'} = ncfloat('depth'); 
nc{'depth'}.FORTRAN_format = ncchar('F10.2');
nc{'depth'}.units = ncchar('m');
nc{'depth'}.type = ncchar('EVEN');
nc{'depth'}.epic_code = nclong(3);
nc{'depth'}.name = ncchar('D');
nc{'depth'}.generic_name = ncchar('depth');
nc{'depth'}.long_name = ncchar('DEPTH (m)');
nc{'depth'}.blanking_distance = ncdouble(cdfb.PCADPProbeBlankDistance(:));
nc{'depth'}.bin_size = ncdouble(cdfb.CellSize(:));
nc{'depth'}.transducer_offset_from_bottom = ncdouble(cdfb.PCADPProbeHeight(:));

nc{'lon'} = ncfloat('lon'); 
nc{'lon'}.FORTRAN_format = ncchar('f10.4');
nc{'lon'}.units = ncchar('degree_east');
nc{'lon'}.type = ncchar('EVEN');
nc{'lon'}.epic_code = nclong(502);
nc{'lon'}.name = ncchar('LON');
nc{'lon'}.long_name = ncchar('LONGITUDE');
nc{'lon'}.generic_name = ncchar('lon');

nc{'lat'} = ncfloat('lat'); 
nc{'lat'}.FORTRAN_format = ncchar('F10.2');
nc{'lat'}.units = ncchar('degree_north');
nc{'lat'}.type = ncchar('EVEN');
nc{'lat'}.epic_code = nclong(500);
nc{'lat'}.name = ncchar('LAT');
nc{'lat'}.long_name = ncchar('LATITUDE');
nc{'lat'}.generic_name = ncchar('lat');

% Velocities - mean for burst
names  = {'u','v','w'};
long_names = {'Eastward','Northward','Vertical'};
epcodes = [1205 1206 1204];
% write the new 2D rotated velocities
for i = 1:length(names),
    ncname = sprintf('%s_%d',names{i},epcodes(i));    
    nc{ncname} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{ncname};
    ncobj.epic_code = nclong(epcodes(i));
    ncobj.name = ncchar(names{i});
    ncobj.long_name = ncchar(['Mean ' long_names{i} ' Velocity']);
    ncobj.generic_name = ncchar(names{i});
    ncobj.units = ncchar('cm/s');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-1000 1000]);
end

% Std. Dev. of Velocities 
% 4097:USTD:STD Eastward velocity     :ustd:cm s-1:f8.2:standard deviation
%                                                        of east component
% 4098:VSTD:STD Northward velocity    :vstd:cm s-1:f8.2:standard deviation
%                                                        of north component
% 4099:WSTD:STD Vertical velocity     :wstd:cm s-1:f8.2:standard deviation
%   
names  = {'USTD','VSTD','WSTD'};
long_names = {'east','north','upward'};
epcodes = [4097 4098 4099]; % just used to name the variable
namestub = {'ustd','vstd','wstd'};
for i = 1:length(names),
    ncname = sprintf('%s_%d',names{i},epcodes(i));    
    nc{ncname} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{ncname};
    ncobj.epic_code = nclong(epcodes(i));
    ncobj.name = ncchar(names{i});
    ncobj.long_name = ncchar(['Std. Dev. of ' long_names{i} ' component']);
    ncobj.generic_name = ncchar(namestub{i});
    ncobj.units = ncchar('cm s-1');
    ncobj.FORTRAN_format = ncchar('f8.2');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-1000 1000]);
end

% Min of Velocities 
names  = {'u','v','w'};
long_names = {'Eastward','Northward','Vertical'};
epcodes = [1205 1206 1204];
namestub = 'min';
for i = 1:length(names),
    ncname = sprintf('%s_%d%s',names{i},epcodes(i),namestub);    
    nc{ncname} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{ncname};
    % No epic code defined for this yet
    %ncobj.epic_code = nclong(epcodes(i));
    ncobj.name = ncchar(sprintf('%s_%s',names{i},namestub));
    ncobj.long_name = ncchar(['Minimum ' long_names{i} ' Velocity']);
    ncobj.generic_name = ncchar(sprintf('%s_%s',names{i},namestub));
    ncobj.units = ncchar('cm/s');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-1000 1000]);
end

% Max of Velocities 
names  = {'u','v','w'};
long_names = {'Eastward','Northward','Vertical'};
epcodes = [1205 1206 1204];
namestub = 'max';
for i = 1:length(names),
    ncname = sprintf('%s_%d%s',names{i},epcodes(i),namestub);    
    nc{ncname} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{ncname};
    % No epic code defined for this yet
    %ncobj.epic_code = nclong(epcodes(i));
    ncobj.name = ncchar(sprintf('%s_%s',names{i},namestub));
    ncobj.long_name = ncchar(['Maximum ' long_names{i} ' Velocity']);
    ncobj.generic_name = ncchar(sprintf('%s_%s',names{i},namestub));
    ncobj.units = ncchar('cm/s');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-1000 1000]);
end

%300:CS :CURRENT SPEED (CM/S)      :vspd:cm s-1:f8.2:oceanographic (going to)
nc{'CS_300'} = ncfloat('time', 'depth', 'lat', 'lon'); 
nc{'CS_300'}.epic_code = nclong(300);
nc{'CS_300'}.name = ncchar('CS');
nc{'CS_300'}.long_name = ncchar('CURRENT SPEED (CM/S)      ');
nc{'CS_300'}.generic_name = ncchar('vspd');
nc{'CS_300'}.units = ncchar('cm s-1');
nc{'CS_300'}.FORTRAN_format = ncchar('f8.2');
nc{'CS_300'}.sensor_type = ncchar('Sontek PCADP');
nc{'CS_300'}.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
nc{'CS_300'}.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
nc{'CS_300'}.serial = ncchar(cdfb.PCADPProbeSerial(:));
nc{'CS_300'}.minimum = ncfloat(0);
nc{'CS_300'}.maximum = ncfloat(0);
nc{'CS_300'}.valid_range = ncfloat([0 1000]);

%310:CD :CURRENT DIRECTION (T)     :vdir:degrees:f8.2:oceanographic (going to), using true north
nc{'CD_310'} = ncfloat('time', 'depth', 'lat', 'lon'); 
nc{'CD_310'}.epic_code = nclong(310);
nc{'CD_310'}.name = ncchar('CD');
nc{'CD_310'}.long_name = ncchar('CURRENT DIRECTION (T)     ');
nc{'CD_310'}.generic_name = ncchar('vdir');
nc{'CD_310'}.units = ncchar('degrees');
nc{'CD_310'}.FORTRAN_format = ncchar('f8.2');
nc{'CD_310'}.sensor_type = ncchar('Sontek PCADP');
nc{'CD_310'}.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
nc{'CD_310'}.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
nc{'CD_310'}.serial = ncchar(cdfb.PCADPProbeSerial(:));
nc{'CD_310'}.minimum = ncfloat(0);
nc{'CD_310'}.maximum = ncfloat(0);
nc{'CD_310'}.valid_range = ncfloat([0 360]);

% Velocities - mean for burst
names  = {'ResU','ResV','ResW'};
long_names = {'Eastward','Northward','Vertical'};
%epcodes = [1205 1206 1204];
% write the new 2D rotated velocities
for i = 1:length(names),
    ncname = sprintf('%s',names{i});    
    nc{ncname} = ncfloat('time', 'lat', 'lon'); ncobj = nc{ncname};
    %ncobj.epic_code = nclong(epcodes(i));
    ncobj.name = ncchar(names{i});
    ncobj.long_name = ncchar(['Mean Resolution Velocity' long_names{i}]);
    ncobj.generic_name = ncchar(names{i});
    ncobj.units = ncchar('cm/s');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([1000 1000]);
end

% include the signal amplitude of each axis
% 1221:AGC1:Echo Intensity (AGC) Beam 1:AGC:counts: :ADCP Beam 1 AGC
% 1222:AGC2:Echo Intensity (AGC) Beam 2:AGC:counts: :ADCP Beam 2 AGC
% 1223:AGC3:Echo Intensity (AGC) Beam 3:AGC:counts: :ADCP Beam 3 AGC
% Signal strength, recorded for each PCADP receiver, is a measure of the intensity of the reflected
% acoustic signal. This is recorded as raw signal amplitude in internal logarithmic units of counts;
% one count is equal to 0.43 dB. With the PCADP software, signal strength can be accessed either as
% signal amplitude in counts or as signal-to-noise ratio (SNR) in dB. SNR is derived from signal
% amplitude by subtracting the ambient noise level and converting to units of dB.
epcodes = [1221 1222 1223];
for i = 1:length(names),
    ncname = sprintf('AGC%1d_%d',i,epcodes(i));    
    nc{ncname} = ncfloat('time','depth', 'lat', 'lon'); ncobj = nc{ncname};
    ncobj.name = ncchar(sprintf('AGC%1d',i));
    ncobj.long_name = ncchar(sprintf('Echo Intensity (AGC) Beam %1d',i));
    ncobj.generic_name = ncchar('AGC');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.units = ncchar('counts');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.valid_range = [0 1000]; % max range not known
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
end

% include the signal correlation of each axis
% 1285:cor1:Beam 1 correlation: :counts:i10:ADCP correlation of beam 1
% 1286:cor2:Beam 2 correlation: :counts:i10:ADCP correlation of beam 2
% 1287:cor3:Beam 3 correlation: :counts:i10:ADCP correlation of beam 3
epcodes = [1285 1286 1287];
for i = 1:length(epcodes),
    ncname = sprintf('cor%1d_%d',i,epcodes(i));    
    nc{ncname} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{ncname};
    ncobj.name = ncchar(sprintf('cor%1d',i));
    ncobj.long_name = ncchar(sprintf(' Beam %1d Correlation',i));
    ncobj.generic_name = ncchar('Cor');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.units = ncchar('percent');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = cdfb.PCADPProbeSerial(:);
    ncobj.valid_range = [0 100]; % max range not known
end

if strcmp(cdfb.PCADPUserSetupPCMode(:),'Yes'),
    % woefully unepic, but needed
    nc{'brange'} = ncfloat('time','beam', 'lat', 'lon'); ncobj = nc{'brange'};
    ncobj.name = ncchar('brange');
    ncobj.long_name = ncchar('Sensor Range to Boundary');
    ncobj.generic_name = ncchar('brange');
    ncobj.epic_code = nclong(0);
    ncobj.units = ncchar('m');
    ncobj.FORTRAN_format = ncchar(' ');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([0 1000]);
    ncobj.NOTE = ncchar('edited by ecorr');
end

if strcmp(cdfb.CompassInstalled(:), 'Yes'),
    % 1215:Hdg :INST Heading             :hdg:degrees:F10.2:ADCP heading
    nc{'Hdg_1215'} = ncfloat('time', 'lat', 'lon'); 
    ncobj = nc{'Hdg_1215'};
    ncobj.name = ncchar('Hdg');
    ncobj.long_name = ncchar('Median 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 = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([0 360]);
    ncobj.NOTE = ncchar('Median of burst heading data');
    
    % 1216:Ptch:INST Pitch               :ptch:degrees:F10.2:ADCP pitch
    nc{'Ptch_1216'} = ncfloat('time', 'lat', 'lon'); 
    ncobj = nc{'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 = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-180 180]);
    ncobj.NOTE = ncchar('Median of burst pitch data');

    % 1217:Roll:INST Roll                :roll:degrees:F10.2:ADCP roll
    nc{'Roll_1217'} = ncfloat('time', 'lat', 'lon'); 
    ncobj = nc{'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 = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-180 180]);
    ncobj.NOTE = ncchar('Median of burst roll data');

    % 1218:HSD :Heading Std. Dev.        :hdg:degrees: :ADCP Hdg. Std. Dev.
    nc{'HSD_1218'} = ncfloat('time', 'lat', 'lon'); 
    ncobj = nc{'HSD_1218'};
    ncobj.name = ncchar('HSD');
    ncobj.long_name = ncchar('Heading Std. Dev.        ');
    ncobj.generic_name = ncchar('hdg');
    ncobj.epic_code = nclong(1218);
    ncobj.units = ncchar('degrees');
    %ncobj.FORTRAN_format = ncchar('f10.2');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([0 360]);
    ncobj.NOTE = ncchar('Calculated from burst heading data');
    
    % 1219:PSD :Pitch Std. Dev.          :ptch:degrees: :ADCP Pitch Std. Dev.
    nc{'PSD_1219'} = ncfloat('time', 'lat', 'lon'); 
    ncobj = nc{'PSD_1219'};
    ncobj.name = ncchar('PSD');
    ncobj.long_name = ncchar('Pitch Std. Dev.          ');
    ncobj.generic_name = ncchar('ptch');
    ncobj.epic_code = nclong(1219);
    ncobj.units = ncchar('degrees');
    %ncobj.FORTRAN_format = ncchar('f10.2');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-180 180]);
    ncobj.NOTE = ncchar('Calculated from burst pitch data');

    % 1220:RSD :Roll Std. Dev.           :roll:degrees: :ADCP Roll Std. Dev.
    nc{'RSD_1220'} = ncfloat('time', 'lat', 'lon'); 
    ncobj = nc{'RSD_1220'};
    ncobj.name = ncchar('RSD');
    ncobj.long_name = ncchar('Roll Std. Dev.             ');
    ncobj.generic_name = ncchar('roll');
    ncobj.epic_code = nclong(1220);
    ncobj.units = ncchar('degrees');
    %ncobj.FORTRAN_format = ncchar('f10.2');
    ncobj.sensor_type = ncchar('Sontek PCADP');
    ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
    ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
    ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-180 180]);
    ncobj.NOTE = ncchar('Calculated from burst roll data');
end

if ~isempty(cdfb{'extpressfreq'}) || ~isempty(cdfb{'extpress'}) || ~isempty(cdfb{'pressure'}),
    % 4023:P   :AVERAGE BURST PRESSURE   :pres:mbar:f10.3:average of burst pressures
    nc{'P_4023'} = ncfloat('time', 'lat', 'lon'); ncobj = nc{'P_4023'};
    ncobj.name = ncchar('P');
    ncobj.long_name = ncchar('AVERAGE BURST PRESSURE   ');
    ncobj.generic_name = ncchar('pres');
    ncobj.FORTRAN_format = ncchar('f10.3');
    ncobj.epic_code = nclong(4023);
    ncobj.units = ncchar('mbar');
    if strcmp(cdfb.ExtPressInstalled(:),'No'),
        ncobj.sensor_type = ncchar('Sontek PCADP Internal');
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
        ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
        ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    else
        ncobj.sensor_type = ncchar(cdfb.ExtPressInstalled(:));
    end
    if strcmp(cdfb.ExtPressInstalled(:),'Druck') || strcmp(cdfb.ExtPressInstalled(:),'ParosFreq'),
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb{'extpressfreq'}.height(:));
        ncobj.initial_sensor_height = ncdouble(cdfb{'extpressfreq'}.height(:));
        ncobj.serial = cdfb{'extpressfreq'}.serial(:);
    elseif strcmp(cdfb.PressInstalled,'Yes'), % internal pressure
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
        ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
        ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    end
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([0 500000]);
    
    % 850:SDP:STAND. DEV. (PRESS)     :pres:mbar:f10.5:std. deviation of burst pressures
    nc{'SDP_850'} = ncfloat('time', 'lat', 'lon'); ncobj = nc{'SDP_850'};
    ncobj.epic_code = nclong(1);
    ncobj.name = ncchar('SDP');
    ncobj.long_name = ncchar('STAND. DEV. (PRESS)     ');
    ncobj.generic_name = ncchar('pres');
    ncobj.units = ncchar('mbar');
    if strcmp(cdfb.ExtPressInstalled(:),'No'),
        ncobj.sensor_type = ncchar('Sontek PCADP Internal');
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
        ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
        ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    else
        ncobj.sensor_type = ncchar(cdfb.ExtPressInstalled(:));
    end
    if strcmp(cdfb.ExtPressInstalled(:),'Druck') || strcmp(cdfb.ExtPressInstalled(:),'ParosFreq'),
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb{'extpressfreq'}.height(:));
        ncobj.initial_sensor_height = ncdouble(cdfb{'extpressfreq'}.height(:));
        ncobj.serial = cdfb{'extpressfreq'}.serial(:);
    elseif strcmp(cdfb.PressInstalled,'Yes'), % internal pressure
        ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
        ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
        ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
    end
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([0 5000]);
end

% 1211:Tx  :PCADP Transducer Temp.    :temp:deg. C:F10.2:ADCP Transducer Temp
nc{'Tx_1211'} = ncfloat('time', 'lat', 'lon');  ncobj = nc{'Tx_1211'};
ncobj.name = ncchar('Tx');
ncobj.long_name = ncchar('Transducer Temp.');
ncobj.generic_name = ncchar('temp');
ncobj.units = ncchar('degrees C');
ncobj.epic_code = nclong(1211);
ncobj.sensor_type = ncchar('Sontek PCADP');
ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
ncobj.minimum = ncfloat(0);
ncobj.maximum = ncfloat(0);
ncobj.valid_range = ncfloat([-5 40]);

% bindist  :distance to bin center from transducer head.  not an EPIC
% variable
nc{'bindist'} = ncfloat('depth'); ncobj = nc{'bindist'};
ncobj.name = ncchar('bindist');
ncobj.long_name = ncchar('Distance from transducer head to bin center');
ncobj.generic_name = ncchar('bindist');
ncobj.units = ncchar('m');
ncobj.sensor_type = ncchar('Sontek PCADP');
ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfb.PCADPProbeHeight(:));
ncobj.initial_sensor_height = ncdouble(cdfb.PCADPProbeHeight(:));
ncobj.serial = ncchar(cdfb.PCADPProbeSerial(:));
ncobj.minimum = ncfloat(0);
ncobj.maximum = ncfloat(0);
ncobj.valid_range = ncfloat([0 5]);

for i=1:2,
    cdfobj = cdfb{sprintf('extsensor%d',i)};
    if ~isempty(cdfobj),
        if findstr('TRANS',cdfobj.serial(:)), % a transmissometer
            ncvarname = sprintf('ATTN%d_55',i);
            %  55:ATTN:ATTENUATION             :attn:m-1:f7.5:added for r2d2 ctd data
            nc{ncvarname} = ncfloat('time','lat', 'lon'); ncobj = nc{ncvarname};
            transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('attn%d',i));
            ncobj.long_name = ncchar(sprintf('ATTENUATION %d',i));
            ncobj.units = ncchar('m-1');
            ncobj.epic_code = nclong(55);
            ncobj.sensor_type = ncchar('transmissometer');
            ncobj.FORTRAN_format = ncchar('F7.5');
            ncobj.sensor_type = ncchar('transmissometer');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 10000]);
            ncobj.varcomment.ATTN_55 = 'ATTN = -(1/focal_length) .* log(tran_4010./(.95*Vair))';
            % 4010:tran:TRANSMISSION (VOLTS)     :trans:volts:f10.3:basic measurement, transmissometer
            ncvarname = sprintf('tran%d_4010',i);
            nc{ncvarname} = ncfloat('time','lat', 'lon'); ncobj = nc{ncvarname};
            transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('tran%d',i));
            ncobj.long_name = ncchar(sprintf('TRANSMISSION (VOLTS) #%d',i));
            ncobj.sensor_type = ncchar('transmissometer');
            ncobj.units = ncchar('volts');
            ncobj.FORTRAN_format = ncchar('F10.3');
            ncobj.sensor_type = ncchar('transmissometer');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 5]);
        elseif findstr('OBS',cdfobj.serial(:)), % an OBS
            ncvarname = sprintf('NEP%d_56',i);
            %  56:NEP:BACKSCATTER INTENSITY    :nephelometer:v:f10.6: added for Dave Pashinski
            nc{ncvarname} = ncfloat('time','lat', 'lon'); ncobj = nc{ncvarname};
            ncals = transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('nep%d',i));
            ncobj.long_name = ncchar(sprintf('BACKSCATTER INTENSITY %d',i));
            ncobj.units = ncchar('v');
            ncobj.generic_name = ncchar('nephelometer');
            ncobj.epic_code = nclong(56);
            ncobj.FORTRAN_format = ncchar('F10.6');
            ncobj.sensor_type = ncchar('nephelometer');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 5]);
            % now, if there's cal info to compute concentration, then set up a variable
            if (ncals > 0) && ~isempty(cdfobj.cal_coef(:)),
                % 981:Sed:Sediment concentration   :sed:g/l:f10.2:sediment concentration
                ncvarname = sprintf('Sed%d_981',i);
                nc{ncvarname} = ncfloat('time','lat', 'lon'); ncobj = nc{ncvarname};
                transfer_calvalues(cdfobj,ncobj);
                ncobj.name = ncchar(sprintf('sed%d',i));
                ncobj.long_name = ncchar(sprintf('Sediment concentration %d',i));
                ncobj.generic_name = ncchar('sed');
                ncobj.epic_code = nclong(981);
                ncobj.units = ncchar('g/l');
                ncobj.sensor_type = ncchar('nephelometer');
                ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
                ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
                ncobj.serial = ncchar(cdfobj.serial(:));
                ncobj.minimum = ncfloat(0);
                ncobj.maximum = ncfloat(0);
                ncobj.valid_range = ncfloat([0 10000]);
            else
                fprintf('No coefficients for ext sensor %d, Sed variable not written\n',i)
            end
        else
            ncvarname = sprintf('extsensor%d',i);
            nc{ncvarname} = ncfloat('time', 'lat', 'lon'); ncobj = nc{ncvarname};
            transfer_calvalues(cdfobj,ncobj);
            ncobj.name = ncchar(sprintf('ext%d',i));
            ncobj.long_name = ncchar(sprintf('external sensor %d',i));
            ncobj.generic_name = ncchar(' ');
            ncobj.epic_code = nclong(0);
            ncobj.units = ncchar('v');
            ncobj.sensor_depth = ncdouble(cdfb.WATER_DEPTH(:)-cdfobj.height(:));
            ncobj.initial_sensor_height = ncdouble(cdfobj.height(:));
            ncobj.serial = ncchar(cdfobj.serial(:));
            ncobj.minimum = ncfloat(0);
            ncobj.maximum = ncfloat(0);
            ncobj.valid_range = ncfloat([0 5]);
        end
    end
end

% data from the external CTD
if ~isempty(cdfs{'CTD_temp'}),
    %4211:CTDTMP :Temperature, ITS-90    :temp:degC:f9.4:
    nc{'CTDTMP_4211'} = ncfloat('time', 'depth', 'lat', 'lon');  ncobj = nc{'CTDTMP_4211'};
    ncobj.name = ncchar('CTDTMP');
    ncobj.long_name = ncchar('TEMPERATURE (C)          ');
    ncobj.generic_name = ncchar('temp');
    ncobj.FORTRAN_format = ncchar('f10.2');
    ncobj.units = ncchar('C');
    ncobj.epic_code = nclong(4211);
    ncobj.sensor_type = ncchar('Seabird SBE37');
    ncobj.sensor_depth = ncdouble(cdfs.WATER_DEPTH(:)-cdfs{'CTD_temp'}.height(:));
    ncobj.initial_sensor_height = ncdouble(cdfs{'CTD_temp'}.height(:));
    ncobj.serial = ncchar(cdfs{'CTD_temp'}.serial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([-5 40]);
    ncobj.comment = ncchar('ITS-1990 standard');

    %4218:CTDCON:CTD Conductivity       :con:S/m:f10.3:
    nc{'CTDCON_4218'} = ncfloat('time', 'depth', 'lat', 'lon');  ncobj = nc{'CTDCON_4218'};
    ncobj.name = ncchar('CTDCON');
    ncobj.long_name = ncchar('CTD Conductivity       ');
    ncobj.generic_name = ncchar('con');
    ncobj.FORTRAN_format = ncchar('f10.3');
    ncobj.units = ncchar('S/m');
    ncobj.epic_code = nclong(4218);
    ncobj.sensor_type = ncchar('Seabird SBE37');
    ncobj.sensor_depth = ncdouble(cdfs.WATER_DEPTH(:)-cdfs{'CTD_cond'}.height(:));
    ncobj.initial_sensor_height = ncdouble(cdfs{'CTD_cond'}.height(:));
    ncobj.serial = ncchar(cdfs{'CTD_cond'}.serial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([0.0 10.0]);

    % check to see if there's any pressure in the pressure
    if any(cdfs{'CTD_press'}(:)),
        %4203:CTDPRS:Pressure                :depth:dbar:f9.1:
        nc{'CTDPRS_4203'} = ncfloat('time', 'depth', 'lat', 'lon');  ncobj = nc{'CTDPRS_4203'};
        ncobj.epic_code = nclong(4203);
        ncobj.name = ncchar('CTDPRS');
        ncobj.long_name = ncchar('Pressure            ');
        ncobj.generic_name = ncchar('depth');
        ncobj.units = ncchar('dbar');
        ncobj.FORTRAN_format = ncchar('f9.1');
        ncobj.sensor_type = ncchar('Seabird SBE37');
        ncobj.sensor_depth = ncdouble(cdfs.WATER_DEPTH(:)-cdfs{'CTD_press'}.height(:));
        ncobj.initial_sensor_height = ncdouble(cdfs{'CTD_press'}.height(:));
        ncobj.serial = ncchar(cdfs{'CTD_press'}.serial(:));
        ncobj.minimum = ncfloat(0);
        ncobj.maximum = ncfloat(0);
        ncobj.valid_range = ncfloat([0 3000]);
        ncobj.transducer_offset_from_bottom = ncdouble(cdfs{'CTD_press'}.height(:));
    end

    % 4214:CTDSAL :CTD Salinity, PSS-78   :sal:PSU:f9.4:
    nc{'CTDSAL_4214'} = ncfloat('time', 'depth', 'lat', 'lon');  ncobj = nc{'CTDSAL_4214'};
    ncobj.name = ncchar('CTDSAL');
    ncobj.long_name = ncchar('CTD Salinity, PSS-78');
    ncobj.generic_name = ncchar('sal');
    ncobj.FORTRAN_format = ncchar('f9.4');
    ncobj.units = ncchar('PSU');
    ncobj.epic_code = nclong(4214);
    ncobj.sensor_type = ncchar('Seabird SBE37');
    ncobj.sensor_depth = ncdouble(cdfs.WATER_DEPTH(:)-cdfs{'CTD_sal'}.height(:));
    ncobj.initial_sensor_height = ncdouble(cdfs{'CTD_sal'}.height(:));
    ncobj.serial = ncchar(cdfs{'CTD_sal'}.serial(:));
    ncobj.minimum = ncfloat(0);
    ncobj.maximum = ncfloat(0);
    ncobj.valid_range = ncfloat([0.0 40.0]);
    ncobj.comment = ncchar('Practical Salinity Units');
end

add_fillvalues(nc,1e35)
add_vardesc(nc)

return

function remove_PCADPglobal_atts(cdf)
% remove global attibutes that are not really needed in the BBV
% Hydra specific

% PPCADP specific
delete(cdf.Nbeams);
delete(cdf.Ncells);

delete(cdf.LISSTInstalled);
delete(cdf.LISSTnBytes);
delete(cdf.PCADPProbeType);
delete(cdf.PCADPProbeBoardRev);
delete(cdf.PCADPProbeBeamGeometry);
delete(cdf.PCADPProbePowerSaveMode);
delete(cdf.PCADPProbeNPingPerBeam);
delete(cdf.PCADPProbeLag);
delete(cdf.PCADPProbePulseLength);
delete(cdf.PCADPProbeRecLength);
delete(cdf.PCADPProbeMinBlankLength);
delete(cdf.PCADPProbeOperatingRange);
delete(cdf.PCADPProbePingDelay);
delete(cdf.PCADPProbeAutoFilter);
delete(cdf.PCADPProbeProfFastPing);
delete(cdf.PCADPProbeProfBlankLenCoh);
delete(cdf.PCADPProbeProfPlenCoh);
delete(cdf.PCADPProbeProfSampLenCoh);
delete(cdf.PCADPProbeProfPlagCoh);
delete(cdf.PCADPProbeProfActualPlagCoh);
delete(cdf.PCADPProbeProfMinCorrLevel);
delete(cdf.PCADPProbeResBlankLenCoh);
delete(cdf.PCADPProbeResPlenCoh);
delete(cdf.PCADPProbeResCellSize);
delete(cdf.PCADPProbeResSampLenCoh);
delete(cdf.PCADPProbeResPlagCoh);
delete(cdf.PCADPProbeResActualPlagCoh);
delete(cdf.PCADPProbeProfilesPerBurst);

delete(cdf.PCADPUserSetupTemp);
delete(cdf.PCADPUserSetupSal);
delete(cdf.PCADPUserSetupCW);
delete(cdf.PCADPUserSetupSensorDepth);
delete(cdf.PCADPUserSetupTempMode);

delete(cdf.PCADPUserSetupCoordSystem);
delete(cdf.PCADPUserSetupOutMode);
delete(cdf.PCADPUserSetupOutFormat);
delete(cdf.PCADPUserSetupRecorderEnabled);
delete(cdf.PCADPUserSetupRecorderMode);
delete(cdf.PCADPUserSetupDeploymentMode);
delete(cdf.PCADPUserSetupDeploymentName);
delete(cdf.PCADPUserSetupBeginDeployment);
delete(cdf.PCADPUserSetupADPComments);
delete(cdf.PCADPUserSetupBottomTrack);
delete(cdf.PCADPUserSetupDeclination);
delete(cdf.PCADPUserSetupRecordPseries);
delete(cdf.PCADPUserSetupRecordPseriesNProfiles);
delete(cdf.PCADPUserSetupRecordPseriesSpectra);
delete(cdf.PCADPUserSetupRecordPseriesSamplingRate);
delete(cdf.PCADPUserSetupRecordPseriesLength);
delete(cdf.PCADPUserSetupRecordPseriesOutMode);
delete(cdf.PCADPUserSetupSlaveMode);
delete(cdf.PCADPUserSetupSyncMode);
return