% 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

% modified to read global attributes from a common file 9/14/07 etm
%
% filename must contain a list of Global attributes
gatt=read_globalatts('..\globatt_837.txt');  % provides values to insert below

% ------------- Step I what are the file names?
% processing file for 728PC002.adp
adpFile = '837R002.adp'; % the raw binary data file

% make sure mooring is a string
if isnumeric(gatt.MOORING)
    gatt.MOORING=num2str(gatt.MOORING);
end
if length(gatt.MOORING) < 4
    gatt.MOORING=[gatt.MOORING outFileRoot(4)];
end
% ------------- 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 = gatt.MOORING; %four digit NNNL, NNN mooring + L logger
metadata.deployment_start = gatt.Deployment_date; % the in water time
metadata.deployment_end = gatt.Recovery_date; % the out of water time
metadata.lon = gatt.longitude; % decimal degrees, West = negative
metadata.lat = gatt.latitude; % decimal degrees, South = negative
metadata.declination = gatt.magnetic_variation; % degrees west is negative
metadata.water_depth = gatt.WATER_DEPTH; % m (DO NOT OMIT!)

metadata.system.serial = 'G357'; % serial number of the Hydra electronics
metadata.system.synch_to = 'Master to ABS'; % G### | nothing
metadata.pcadp.serial = 'H96';
metadata.pcadp.height = 0.985; % 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 = 'T52103'; % use the druck file name without the extension
metadata.extpress.height = 0.985; % 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
% anything under the metadat.ext#.cals structure gets written as an attibute
% in the netCDF file
metadata.ext1.serial = 'TRANS795'; % OBS### | TRANS### | XXX####
metadata.ext1.height = 2.24; % m
% the following are required for trans data
metadata.ext1.cals.path_length = 0.25; % focal length of the trans, 0.5 | 0.25 | 1 m, required!
metadata.ext1.cals.pre = [4.695 0.055]; % for trans, air reading: [unblocked blocked]
metadata.ext1.cals.post = [4.571 0.055]; % for trans, air reading: [unblocked blocked]
metadata.ext2.serial = 'OBS2277'; % OBS### | TRANS### | XXX####
metadata.ext2.height = 0.47; % m
metadata.ext2.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.ext2.cals.coef = [0.0255 134.5516]; % coefficients to a calibration eqation
% more than three coefficients will be accommodated as higher powers.
% these are not used but are informative
metadata.ext2.cals.units = 'NTU'; % or g/l
metadata.ext2.cals.equation = 'Conc [NTU] = coef(1)+coef(2)*V';
metadata.CTD.serial = ''; %
metadata.CTD.height = 0; % 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 =  '0.0';
metadata.metafile_author = 'EM'
% if the following fields are missing, UNKOWN will be used
metadata.origin = gatt.DATA_ORIGIN; % with collaborator's data, could be USC, etc.
metadata.experiment = gatt.EXPERIMENT;
metadata.description = gatt.DESCRIPTION;
metadata.conventions = gatt.Conventions;
metadata.project = gatt.PROJECT; % 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 19 to, 1527 or say 1 to Inf for all data
%metadata.goodbursts = [1 Inf];
metadata.goodbursts = [19 1527];


% set the root name for the output files
outFileRoot = [gatt.MOORING 'pca']; 

% 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'])
    disp(sprintf('copying %s to %s',[outFileRoot 'q.cdf'], [outFileRoot 'q1.cdf']))
    copyfile([outFileRoot 'q.cdf'], [outFileRoot 'q1.cdf'])
end

% ---------- Step IV semi-automated cleaning
if 0,
    % 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);
end
if 0
    % 2) run just deglitch1vector with nsd=2.8:
    %      deglitch1vector only operates on vars x,y,z
    disp(sprintf('copying %s to %s',[outFileRoot 'b1.cdf'], [outFileRoot 'b2.cdf']))
    copyfile([outFileRoot 'b1.cdf'], [outFileRoot 'b2.cdf'])
    disp(sprintf('copying %s to %s',[outFileRoot 's1.cdf'], [outFileRoot 's2.cdf']))
    copyfile([outFileRoot 's1.cdf'], [outFileRoot 's2.cdf'])
    disp(sprintf('copying %s to %s',[outFileRoot 'q1.cdf'], [outFileRoot 'q2.cdf']))
    copyfile([outFileRoot 'q1.cdf'], [outFileRoot 'q2.cdf'])
    % this data has a stormy section, so the "normal" threshold of 7 was
    % too small.
     seta.nsd=3.2; seta.threshold=18;  %(threshold from stdUres in stats file)
       cleanhydra([outFileRoot 'b2.cdf'], [outFileRoot 's2.cdf'], ...
        [outFileRoot 'q2.cdf'], 'deglitch1vector','settings', seta);

    % here are some other things you can do with cleanhydra:
    % to modify the default settings  (ndt expands the number of points
    % over which the filter is applied- default is 6), 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

if 0    %  USE CLEAN HYDRA TO CORRECT BRANGE
    disp(sprintf('copying %s to %s',[outFileRoot 'b2.cdf'], [outFileRoot 'b3.cdf']))
    copyfile([outFileRoot 'b2.cdf'], [outFileRoot 'b3.cdf'])
    disp(sprintf('copying %s to %s',[outFileRoot 's2.cdf'], [outFileRoot 's3.cdf']))
    copyfile([outFileRoot 's2.cdf'], [outFileRoot 's3.cdf'])
    disp(sprintf('copying %s to %s',[outFileRoot 'q2.cdf'], [outFileRoot 'q3.cdf']))
    copyfile([outFileRoot 'q2.cdf'], [outFileRoot 'q3.cdf'])
end
if 0
    seta.min=700; seta.max=950;
    seta.npts=25;
    seta.nstds=2.4;
    cleanhydra([outFileRoot 'b3.cdf'], [outFileRoot 's3.cdf'], ...
        [outFileRoot 'q3.cdf'], 'fix_vbrange','settings', seta, ...
        'variables',{'MeanRange1', 'MeanRange2', 'MeanRange3'});
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 'b3.cdf'], [outFileRoot 's3.cdf'], [outFileRoot 'A'], switches, [19 1527]);
end

diary off