% metaPCADPExample - an example script for processing a single PCADP data file.  
% The steps here follow the steps outlined in the manual
%
% you should be running this script from the directory in which the data
% files reside

% Written 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

% ------------- Step I what are the file names?
% processing file for 728PC002.adp
adpFile = '728PC002.adp'; % the raw binary data file
outFileRoot = '728pc'; % the name for the output files

% ------------- Step II what is the metadata?
%
% here is the metadata structure for this deployment
% if you use an empty char string denoted by ''
% for an external sensor serial number, then that data will not be included
% calibration info will be assigned to variable attributes dynamically
% based on the field names of items under the struct field cals.
% Do not nest more structs under cals, or adr2cdf will not record the data
% no two data loggers can have the same four digit mooring number
metadata.mooring_number = '7282'; %four digit NNNL, NNN mooring + L logger
metadata.deployment_start = '28-Oct-2003 13:27:00'; % the in water time
metadata.deployment_end = '21-jan-2004 21:09:00'; % the out of water time
metadata.lon = -78.77912; % decimal degrees, West = negative
metadata.lat = 33.72247; % decimal degrees, South = negative
metadata.declination = -8.22; %degrees west is negative
metadata.water_depth = 10.6; % m (DO NOT OMIT!)
metadata.system.serial = 'G136'; % serial number of the Hydra electronics
metadata.system.synch_to = 'nothing'; % G### | nothing
metadata.pcadp.serial = 'H30';
metadata.pcadp.height = 1.10; % height of probe above bottom, m
                                % used to reference anything coming from the PCADP probe
% for attached sensors, use '' for the serial number if there is no sensor, 
% and variable will not be written
metadata.extpress.serial = 'T60005'; % use the druck file name without the extension
metadata.extpress.height = 1.63; % m, same as adv if in probe
% pressure sensor calibration data will be read from the .drk file
% now add the external sensors, in this example, the OBS was on channel 1
% and the trans on channel 2
% anything under the metadat.ext#.cals structure gets written as an attibute
% in the netCDF file
metadata.ext1.serial = 'OBS2181'; % OBS### | TRANS### | XXX####
metadata.ext1.height = 0.6; % m
metadata.ext1.cals.gain = 0.33; % for OBS gain setting, mg/l per volt
% if coef is present, dohydraoptic will apply the coefficients to OBS data
% use [] if there is no calibration for concentration
metadata.ext1.cals.coef = []; % coefficients to a calibration eqation
% more than three coefficients will be accommodated as higher powers.
% these are not used but are informative
metadata.ext1.cals.units = 'kg/m-3'; % or g/l
metadata.ext1.cals.equation = 'Conc [kg/m-3] = coef(1)+coef(2)*V+coef(3)*V^2';
% anything under the metadat.ext#.cals structure gets written as an attibute
% in the netCDF file
metadata.ext2.serial = 'TRANS182'; % OBS### | TRANS### | XXX####
metadata.ext2.height = 1.20; % m
% the following are required for trans data
metadata.ext2.cals.path_length = 0.25; % focal length of the trans, 0.5 | 0.25 | 1 m, required!
metadata.ext2.cals.pre = [0 0]; % for trans, air reading: [unblocked blocked]
metadata.ext2.cals.post = [0 0]; % for trans, air reading: [unblocked blocked]
metadata.CTD.serial = 'SBE1779'; % 
metadata.CTD.height = 0.96; % m
% SBE37's record data in real units and do not need applied cals
metadata.metafile = mfilename('fullpath');
metadata.metafile_date = datestr(now);
metadata.metafile_version =  '1.0';
metadata.metafile_author = 'MM'
% if the following fields are missing, UNKOWN will be used
metadata.origin = 'USGS/WHFC'; % with collaborator's data, could be USC, etc.
metadata.experiment = 'Boston Longterm'; 
metadata.project = 'WHFC'; % might also use OFA funding agency, such as MWRA, EPA, WCMG

% it is useful to omit bursts at the beginning and the end that were not recorded
% when the instrument was in the water
% good bursts are 29 to 2423, or say 1 to Inf for all data
% metadata.goodbursts = [1 Inf];
metadata.goodbursts = [29 to 2423];

% these programs take a long time to run, so keep a diary
diary(sprintf('run%s',datestr(now,30)))

% ---------- Step III convert from vendor's raw binary to netCDF & get quality info
if 0,
    % known good in the water bursts are ??
    % it is useful to omit bursts at the beginning and the end that were not recorded
    % when the instrument was in the water
    diagnostics.verbose = 0;
    diagnostics.pcount = 0;
    diagnostics.fixpfreq = 0;
    diagnostics.rotationSTDs = 0; % output per profile data from profile header
    diagnostics.temperatureSTDs = 0; % output per profile data from profile header
    adp2cdf(adpFile,outFileRoot,metadata.goodbursts,metadata.goodbursts,diagnostics,metadata);
end

if 0,
    % make a copy of the files, fixbadadv works directly on them
    % use the already canged *s1.cdf to be changed further
    disp(sprintf('copying %s to %s',[outFileRoot 'b.cdf'], [outFileRoot 'b1.cdf']))
    copyfile([outFileRoot 'b.cdf'], [outFileRoot 'b1.cdf'])
    disp(sprintf('copying %s to %s',[outFileRoot 's.cdf'], [outFileRoot 's1.cdf']))
    copyfile([outFileRoot 's.cdf'], [outFileRoot 's1.cdf'])
end

% ---------- Step IV semi-automated cleaning


if 1,
    % function cleanhydra(cdfbFile, cdfsFile, cdfqFile, operation, ...
    %       'variables', {vnames}, 'settings', settings.n, ...
    %       'burstrange', [1:Inf]) 
    % Using the two interations below is reasonably robust over multiple data sets
    % 1) run thumbfinger with std of 10:
    %     thumbfinger operates on vars: x,v,z,heading, pitch, roll, and temperature
      seta.nsd=10;
       cleanhydra([outFileRoot 'b1.cdf'], [outFileRoot 's1.cdf'], ...
        [outFileRoot 'q1.cdf'], 'thumbfinger','settings',seta);
    % 2) run just deglitch1vector with nsd=2.8:
    %      deglitch1vector only operates on vars x,y,z
      seta.nsd=2.8;
       cleanhydra([outFileRoot 'b1.cdf'], [outFileRoot 's1.cdf'], ...
        [outFileRoot 'q1.cdf'], 'deglitch1vector','settings', seta);
 
    % here are some other things you can do with cleanhydra:
    % to modify the default settings, and just do bursts 1361-1365, use this syntax:
    %  seta.nsd=3; seta.ndt=10;
    %   cleanhydra([outFileRoot 'b1.cdf'], [outFileRoot 's1.cdf'], [outFileRoot 'q1.cdf'], ...
    %    'deglitch1vector','settings',seta,'burstrange', [1361:1365]);
    %
    %  or replace the first 5 values of pressure with the mean of the rest of the burst
    %  rm_badsamp can operate on ANY variable in the file, extsensors etc
    %    seta.var='pres'; seta.rvalue='mean'; seta.samps=[1:5];
    %   cleanhydra([outFileRoot 'b1.cdf'], [outFileRoot 's1.cdf'], [outFileRoot 'q1.cdf'], ...
    %    'rm_badsamp','settings',seta);
        
end

% ------------ Step V interactive flagging
% the PCADP does not have an equivalent tool to flagbadadv
% use browse hydraburst and ncbrowse to assess the data

% ---------- Step VI convert to scientific units and apply rotations
if 0,
    switches.ecorr_showplots = 0;
    switches.resolve_ambiguity = 1; % 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 = [];
    % so far, have not had pressure spike problems with the PCADP, this
    % override is not available for PCADP: switches.override_pmask = []; 

    % another time consuming step, so keep a diary
    diary(sprintf('run%s',datestr(now,30)))
    Pcadp = pcadp2nc([outFileRoot 'b.cdf'], [outFileRoot 's.cdf'], outFileRoot, switches, [1 Inf]);

end

diary off