%AQTYRAW2CDF convert Aquatec Aquatalk200TY for 2XX TY ascii output to netcdf
% function aqtyraw2cdf(csvfile, cdffile, settings);
% csvfile = raw data exported from Aquatalk reader program
% cdffile = name ofthe output file
% settings:
%     settings.MOORING = 'none'; %four digit NNNL, NNN mooring + L logger
%     settings.Deployment_date = '11 jun 2008'; % the in water time
%     settings.Recovery_date = '7 jul 2008'; % the out of water time
%     settings.longitude = 0; % decimal degrees, West = negative
%     settings.latitude = 0; % decimal degrees, South = negative
%     settings.magnetic_variation = 0; % degrees west is negative
%     settings.WATER_DEPTH = 3; % m (DO NOT OMIT!)
%     settings.WATER_DEPTH_NOTE = 'guessed at low tide'; 
%     settings.DATA_ORIGIN = 'WHSC'; % with collaborator's data, could be USC, etc.
%     settings.EXPERIMENT = 'dock test of wiper'; 
%     settings.DESCRIPTION = ''; 
%     settings.PROJECT = 'none'; % might also use OFA funding agency, such as MWRA, EPA, WCMG
%     settings.Data_Cmnt = ''; 
%     settings.SciPi = 'Martini'; 
%     settings.INST_TYPE = 'Aquatec Aqualogger 200TY';
%     settings.sensor_height = 3; % m of bed
%     settings.OBS_serial_number = '';

% How to export the data this file reads:
%   Download data from an Aquatec TY logger using Aquatalk 200TY

% Written by Marinna Martini 8 July 2008
% for the USGS Woods Hole Science Center

%%% START USGS BOILERPLATE -------------%
% Use of this program is described in:
%
% Program written in Matlab v7.0.0
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%
%%% END USGS BOILERPLATE --------------

function aqtyraw2cdf(csvfile, cdffile, theMeta)

% get the current SVN version- the value is automatically obtained in svn
% is the file's svn.keywords which is set to "Revision"
rev_info = 'SVN $Revision: $';
disp(sprintf('%s %s running',mfilename,rev_info))

% -------- verify inputs
if ~exist('csvfile','var'),
    [csvfile, path] = uigetfile('*.Asc','Open aanderaa 5059 ASCII export of RCM raw data');
    if isempty(csvfile), return; end
    csvfile= [path csvfile];
end
if ~exist('cdffile','var'), % if not specified, make it up
    [path,name] = fileparts(csvfile);
    cdffile = fullfile(path, [name '.cdf']);
end
if ~exist('theMeta','var'), % if not specified, abort
    disp('metadata missing')
    return
end

% global attributes which change
gattnames = {'MOORING', 'Deployment_date', 'Recovery_date', 'longitude', 'latitude',...
    'magnetic_variation', 'WATER_DEPTH', 'DATA_ORIGIN', 'EXPERIMENT','DESCRIPTION',...
    'PROJECT', 'WATER_DEPTH_NOTE','INST_TYPE'};
for ifield = 1:length(gattnames),
    if ~isfield(theMeta,gattnames{ifield}),
        disp(sprintf('Missing metadata for %s',gattnames{ifield}))
        eval(sprintf('theMeta.%s = [];',gattnames{ifield}))
    end
end
% fixed globals
theMeta.Conventions = 'PMEL-EPIC'; % ## variable conventions, ie: EPIC'
theMeta.DATA_SUBTYPE = 'MOORED'; % ## description of DATA_TYPE
theMeta.COORD_SYSTEM = 'GEOGRAPHIC'; % ## coordinate system of the data   '
theMeta.WATER_MASS = '?'; % ## water mass flag used for EPIC contouring programs
theMeta.POS_CONST = 0; % ## consistent position flag (=1 not consistent)
theMeta.DEPTH_CONST = 0; % ## consistent depth flag (=1 not consistent)
theMeta.DRIFTER = 0; % ## drifter flag (=1 if drifter)
theMeta.VAR_FILL = 1.E35; % ## missing or bad data value identifier
theMeta.FILL_FLAG = 0; % ## data fill flag (=1 if data has fill values)
theMeta.COMPOSITE = 0; % ## number of pieces in composite series
theMeta.OBS_manufacturer = 'Seapoint';

% get all the data at once, these are not huge files
fid = fopen(csvfile);
disp('Reading the csv file...')
tic

nheaderlines = 9;

% count the lines
buf = fgetl(fid);
nrec = 0;
while ischar(buf),
    nrec = nrec + 1;
    buf = fgetl(fid);
end
fseek(fid,0,'bof');
nrec = nrec - nheaderlines;
disp(sprintf('Found %d records',nrec))

% do header
fgetl(fid);
theMeta.readaqua_ver = fscanf(fid,'VERSION,%f');
% ,Serial Number, Log Start Date, Log Interval,Burst Length,Intersample, NumBursts,Description, Gain , Startup,Settling
fgetl(fid);
fgetl(fid);
buf = fgetl(fid);
% LOGGER,371-025,00:00:00 06/10/2008,10 Minutes,8 Samples,0 ms,1 Bursts,DockTest,AUTO, 1000 ms, 1000 ms
% parse it
[data, buf] = strtok(buf,','); % LOGGER
[theMeta.serial, buf] = strtok(buf,',');
[theMeta.start_time, buf] = strtok(buf,',');
[data, buf] = strtok(buf,','); 
theMeta.log_interval = sscanf(data,'%d Minutes');
[data, buf] = strtok(buf,','); 
theMeta.burst_length = sscanf(data,'%d Samples');
[data, buf] = strtok(buf,','); 
theMeta.intersample_ms = sscanf(data,'%d ms');
[data, buf] = strtok(buf,','); 
theMeta.numbursts = sscanf(data,'%d Bursts');
[theMeta.user_description, buf] = strtok(buf,','); 
[theMeta.gain_setting, buf] = strtok(buf,','); 
[data, buf] = strtok(buf,','); 
theMeta.startup_time_ms = sscanf(data,'%d ms');
data = strtok(buf,','); 
theMeta.settling_time_ms = sscanf(data,'%d ms');
fgetl(fid);
fgetl(fid);
fgetl(fid);
fgetl(fid);
fgetl(fid);
disp(theMeta)

irec = 1;
tj = ones(nrec,1).*NaN;
rawturb = tj;
turb = tj;
gain = tj;
recnum = tj;
%for i=1:40,
buf = fgetl(fid);
while ischar(buf),
    data = [];
    [data, count] = sscanf(buf,'DATA,%d/%d/%d %d:%d:%d,%d,%f,%d');
    if count ~= 9, % it's not what we are looking for
        % DATA,6/10/2008 1:10,0,0,100
        [data, count] = sscanf(buf,'DATA,%d/%d/%d %d:%d,%d,%f,%d');
        data = [data(1:5); 0; data(6:8)]';
    elseif count ~=8,
        disp('Data does not match expected format')
    end
    if ~isempty(data),
        % record length is wrong (misses time every midnight) etm mods below
        tj(irec) = julian([data(3) data(1) data(2) data(4) data(5) data(6)]);
        rawturb(irec) = data(7);
        turb(irec) = data(8);
        gain(irec) = data(9);
        recnum(irec) = irec;
    end
    irec = irec + 1;

    buf = fgetl(fid);
    if ~rem(irec,1000), disp(sprintf('Converted %d records',irec)); end
end
fclose(fid);

cdf = netcdf(cdffile, 'clobber');
if isempty(cdf), return, end

definencfile(cdf,theMeta);

%endef(cdf);

cdf{'time'}(1:length(tj)) = floor(tj);
cdf{'time2'}(1:length(tj)) = (tj-floor(tj)).*(1000*24*3600);
cdf.start_time = ncchar(datestr(datenum(gregoria(tj(1)))));
cdf.stop_time = ncchar(datestr(datenum(gregoria(tj(end)))));
%cdf.DELTA_T = ncchar(num2str(theMeta.log_interval));
cdf.DELTA_T = ncchar(num2str(mean(diff(tj)).*(24*3600)));
cdf{'depth'}(1) = theMeta.WATER_DEPTH;
cdf{'lon'}(1) = theMeta.longitude;
cdf{'lat'}(1) = theMeta.latitude;
cdf{'recnum'}(1:length(recnum)) = recnum;
cdf{'NEP_56'}(1:length(rawturb)) = rawturb;
cdf{'Sed_981'}(1:length(turb)) = turb;
cdf{'gain'}(1:length(gain)) = gain;

add_vardesc(cdf);
add_fillvalues(cdf,1E35);

cdf.history = ncchar(sprintf('Created by %s %s, declination of %f applied.',...
    mfilename, rev_info,theMeta.magnetic_variation));
add_minmaxvalues(cdf);

disp(sprintf('Wrote %d records from %s to %s', length(cdf{'time'}(:)),...
    cdf.start_time(:), cdf.stop_time(:)))

close(cdf);

return % from main

% ********************************** subfunctions ********************************

function definencfile(cdf,theMeta)

% separate some stuff we know we don't want in the global attributes
omitGatts = {'sensor_height','OBS_serial_number'};
for i = 1:length(omitGatts),
    if isfield(theMeta, omitGatts{i}),
        eval(sprintf('miscMeta.%s = theMeta.%s;',omitGatts{i},omitGatts{i}))
        theMeta = rmfield(theMeta, omitGatts{i});
    end
end
% write metadata
metaFields = fieldnames(theMeta);
for i = 1:length(metaFields)
    %theField = metaFields{i};
    %theFieldDef = getfield(theMeta,theField);
    theFieldDef = theMeta.(metaFields{i});
    nRows = size(theFieldDef,1);
    nCols = size(theFieldDef,2);
    temp = [];
    if ischar(theFieldDef)
        if nRows>1
            for ii = 1:nRows
                temp = [temp,' ',theFieldDef(ii,1:nCols)];
            end
            theFieldDef = temp;
        end
        eval(['cdf.',metaFields{i},' = ncchar(theFieldDef);'])
    else
        eval(['cdf.',metaFields{i},'= ncfloat(theFieldDef);'])
    end
end

% write additional metadata
cdfb.CREATION_DATE = ncchar(datestr(now));
cdfb.DESCRIPT = ncchar('Aquascat Aqualogger 2XXTY raw data file');
cdfb.DATA_TYPE = ncchar('Aqualogger');

%% Dimensions:
cdf('time') = 0; %% (record dimension)
cdf('depth') = 1;
cdf('lon') = 1;
cdf('lat') = 1;

%% Variables and attributes:
cdf{'time'} = nclong('time');
cdf{'time'}.FORTRAN_format = ncchar('F10.2');
cdf{'time'}.units = ncchar('True Julian Day');
cdf{'time'}.type = ncchar('UNEVEN');
cdf{'time'}.epic_code = nclong(624);

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

cdf{'depth'} = ncfloat('depth'); %% 1 element.
cdf{'depth'}.FORTRAN_format = ncchar('F10.2');
cdf{'depth'}.units = ncchar('m');
cdf{'depth'}.type = ncchar('EVEN');
cdf{'depth'}.epic_code = nclong(3);
cdf{'depth'}(:) = theMeta.WATER_DEPTH;

cdf{'lon'} = ncfloat('lon'); %% 1 element.
cdf{'lon'}.FORTRAN_format = ncchar('F10.2');
cdf{'lon'}.units = ncchar('degrees_east');
cdf{'lon'}.type = ncchar('EVEN');
cdf{'lon'}.epic_code = nclong(502);
cdf{'lon'}(:) = theMeta.longitude;

cdf{'lat'} = ncfloat('lat'); %% 1 element.
cdf{'lat'}.FORTRAN_format = ncchar('F10.2');
cdf{'lat'}.units = ncchar('degrees_north');
cdf{'lat'}.type = ncchar('EVEN');
cdf{'lat'}.epic_code = nclong(500);
cdf{'lat'}(:) = theMeta.latitude;

cdf{'recnum'} = ncfloat('time');
cdf{'recnum'}.name = ncchar('recnum');
cdf{'recnum'}.long_name = ncchar('record number          ');
cdf{'recnum'}.generic_name = ncchar('rec');
cdf{'recnum'}.units = ncchar('count');
cdf{'recnum'}.valid_range = ncfloat([0 262100]);
cdf{'recnum'}.serial_number = ncchar(theMeta.serial);

%  56:NEP:BACKSCATTER INTENSITY    :nephylometer:v:f10.6: added for Dave Pashinski
% TODO - this is not exactly correct as we have counts
ncvarname = 'NEP_56';
cdf{ncvarname} = ncfloat('time'); cdfobj = cdf{ncvarname};
cdfobj.name = ncchar('nep');
cdfobj.long_name = ncchar('BACKSCATTER INTENSITY');
cdfobj.generic_name = ncchar('nephylometer');
cdfobj.epic_code = nclong(56);
%cdfobj.units = ncchar('v');
cdfobj.units = ncchar('count');
cdfobj.sensor_type = ncchar('nephylometer');
cdfobj.sensor_depth = ncdouble(cdf.WATER_DEPTH(:)-miscMeta.sensor_height);
cdfobj.initial_sensor_height = ncdouble(miscMeta.sensor_height);
cdfobj.serial = ncchar(miscMeta.OBS_serial_number);
cdfobj.minimum = ncfloat(0);
cdfobj.maximum = ncfloat(0);
%cdfobj.valid_range = ncfloat([0 5]);
cdfobj.valid_range = ncfloat([0 2^16]);

% 981:Sed:Sediment concentration   :sed:g/l:f10.2:sediment concentration
ncvarname = 'Sed_981';
cdf{ncvarname} = ncfloat('time'); cdfobj = cdf{ncvarname};
cdfobj.name = ncchar('sed');
cdfobj.long_name = ncchar('Sediment concentration');
cdfobj.generic_name = ncchar('sed');
cdfobj.epic_code = nclong(981);
%cdfobj.units = ncchar('g/l');
cdfobj.units = ncchar('NTU');
cdfobj.sensor_type = ncchar('nephylometer');
cdfobj.sensor_depth = ncdouble(cdf.WATER_DEPTH(:)-miscMeta.sensor_height);
cdfobj.initial_sensor_height = ncdouble(miscMeta.sensor_height);
cdfobj.serial = ncchar(miscMeta.OBS_serial_number);
cdfobj.minimum = ncfloat(0);
cdfobj.maximum = ncfloat(0);
cdfobj.valid_range = ncfloat([0 10000]);

ncvarname = 'gain';
cdf{ncvarname} = ncfloat('time'); cdfobj = cdf{ncvarname};
cdfobj.name = ncchar('gain');
cdfobj.long_name = ncchar('sensor gain');
cdfobj.generic_name = ncchar('gain');
cdfobj.epic_code = nclong(0);
%cdfobj.units = ncchar('g/l');
cdfobj.units = ncchar('percent');
cdfobj.sensor_type = ncchar('nephylometer');
cdfobj.sensor_depth = ncdouble(cdf.WATER_DEPTH(:)-miscMeta.sensor_height);
cdfobj.initial_sensor_height = ncdouble(miscMeta.sensor_height);
cdfobj.serial = ncchar(miscMeta.OBS_serial_number);
cdfobj.minimum = ncfloat(0);
cdfobj.maximum = ncfloat(0);
cdfobj.valid_range = ncfloat([0 100]);

return

function add_vardesc(cdf)
% construct the EPIC attribute var_desc on the fly
theVars = var(cdf);
theDims = dim(cdf);
var_desc = '';
for ivar = 1:length(theVars),
    adim = 0;
    for idim = 1:length(theDims),
        if strcmp(ncnames(theVars{ivar}),ncnames(theDims{idim}))...
                || strcmp(ncnames(theVars{ivar}),'time2'),
            adim = 1;
        end
    end
    if ~adim,
        var_desc = sprintf('%s:%s',var_desc,theVars{ivar}.name(:));
    end
end
var_desc = var_desc(2:end);
cdf.VAR_DESC = ncchar(var_desc);
disp(sprintf('Writing variables %s',var_desc));
return

function add_fillvalues(cdf,fill)
theVars = var(cdf);
% don't add _FillValue for coordinate vars (time, time2, lat, lon, depth)
for i = 6:length(theVars),
    switch datatype(theVars{i}),
        case 'float'
            theVars{i}.FillValue_ = ncfloat(fill);
        case 'double'
            theVars{i}.FillValue_ = ncdouble(fill);
        case 'long'
            theVars{i}.FillValue_ = nclong(fill);
    end
end
return

% ---------------------------------------------
function add_minmaxvalues(cdf)
theVars = var(cdf);
for i = 1:length(theVars),
    if ~strcmp(ncnames(theVars{i}),'time') && ~strcmp(ncnames(theVars{i}),'time2'),
        data = theVars{i}(:);
        [row, col] = size(data);
        if col == 1,
            theVars{i}.minimum = ncfloat(min(data));
            theVars{i}.maximum = ncfloat(max(data));
        elseif col == 2,
            theVars{i}.minimum = ncfloat(min(min(data)));
            theVars{i}.maximum = ncfloat(max(max(data)));
        elseif col > 2, % then these are amp and cor
            for icol = 1:col,
                mins(icol) = min(data(:,icol));
                maxs(icol) = max(data(:,icol));
            end
            theVars{i}.minimum = ncfloat(mins);
            theVars{i}.maximum = ncfloat(maxs);
        end
    end
end
return