function [VEL] = AWAC2cdf(baseFile,cdfFile,metadata)
% function [VEL] = AWAC2cdf(basefile,outFile,metadata,timecorr)
%
%   Function to read in exported ascii current data and wave parameters
%   from Nortek AWAC. Data need to be converted to ascii from raw data file
%   using Nortek data conversion utility, and wave parameters need to be
%   calculated using Nortek Quickwave. This function will give a netCDF
%   datafile with all raw data, in raw coordinate system. Conversion to
%   earth coordinates from BEAM or XYZ coordinates and correction for
%   magnetic variation takes place in AWAC2nc.
%   It will also give output structures of current data and instrument
%   metadata. Depends upon all files output from the data conversion
%   utility, including the header file (.hdr), sensor file (.sen) and
%   amplitude (.a*) and current (.v*) data.
%
%   Inputs:    
%       'basefile' - base filename used for instrument deployment, upon
%       which
%       all data files are named
%       
%       'outFile' - your output filename 
%
% Program written in Matlab R2007
% Program updated in Matlab R2010a
% Program ran on Apple MacIntosh OS 10.6.7
%
% "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."

if nargin<1
    help(mfilename)
    return
end

% get the metadata about setup and log from deployment
files = dir;
names = cell(length(files),1);
for i = 1:length(files)
    names{i} = files(i).name;
end
if ~isempty(strmatch([baseFile '.hdr'],names))
    instmeta = readAWACHeader(baseFile);
else
    disp('Need the aquadopp header file (.hdr) for setup information')
    return
end
if ~isempty(strmatch([baseFile '.log'],names))
    logmeta = readAWAClog(baseFile);
else
end
%now load sensor data file for times and sensor data
senfile = strcat(baseFile,'.sen');
SEN = load(senfile);
VEL.battery = SEN(:,9);VEL.soundspeed = SEN(:,10);VEL.heading = SEN(:,11);
VEL.pitch = SEN(:,12);VEL.roll = SEN(:,13);VEL.pressure = SEN(:,14);
VEL.temperature = SEN(:,15);VEL.error = SEN(:,7);VEL.status = SEN(:,8);
% now get the external sensor data if it was logged, and convert to volts
if isfield(metadata,'AnalogInput1'),VEL.AnaInp1 = SEN(:,16)*5/65535;end
if isfield(metadata,'AnalogInput2'),VEL.AnaInp2 = SEN(:,17)*5/65535;end


%now get velocity and amplitude data from ascii output
disp('Reading in ASCII data files')
%now get velocity and amplitude data from ascii output
for i = 1:3
    afile = [baseFile,'.a',num2str(i)];aname = ['Vamp',num2str(i)];
    VEL.(aname) = load(afile);
    vfile = [baseFile,'.v',num2str(i)];vname = ['Vvel',num2str(i)];
    VEL.(vname) = load(vfile);
end
% a1file = strcat(baseFile,'.a1');a2file = strcat(baseFile,'.a2');a3file = strcat(baseFile,'.a3');
% VEL.Vamp1 = load(a1file);VEL.Vamp2 = load(a2file);VEL.Vamp3 = load(a3file);
% v1file = strcat(baseFile,'.v1');v2file = strcat(baseFile,'.v2');v3file = strcat(baseFile,'.v3');
% VEL.Vvel1 = load(v1file);VEL.Vvel2 = load(v2file);VEL.Vvel3 = load(v3file);

[N,M] = size(VEL.Vamp1);            
%create depth vector
blank = instmeta.BlankingDistance; 
blank2 = instmeta.BlankingDistance + metadata.transducer_offset_from_bottom; 
bin = instmeta.CellSize/100;DEP = metadata.WATER_DEPTH;
%Nortek lists the distance to the center of the first bin as the blanking
%distance plus one cell size
VEL.bindist = blank+(bin):bin:(bin*(M-1))+blank+(bin);	%depth vector for plotting distance from bottom
VEL.Depths = fliplr(DEP-((bin*(M-1))+blank2+(bin)):bin:DEP-(blank2+(bin)));	%depth vector for plotting by depth

%generally there is no clock drift in the data file with this instrument 
%(that's not to say the clock doesn't drift), so a straight
%translation to Julian Day is OK.  However, when the battery gets low, the
%instrument will shut off and restart at uneven intervals, so you need to
%check the .cdf file with the times to see if you get a complete record.
%Trimming the file based on start-stop cycle takes place in the conversion
%to .nc file.
jd = julian([SEN(:,3) SEN(:,1) SEN(:,2) SEN(:,4) SEN(:,5) SEN(:,6)]);
%look for small gaps in the record 
hic = find(round(diff(jd)*86400)~=instmeta.ProfileInterval);
if ~isempty(hic)
    disp(['There are ',num2str(length(hic)),' of gaps in the record'])
else
    disp('There are no gaps in the record')
end

if ~isempty(hic) && metadata.fixgaps 
    disp('User instructed to fill gaps with NaNs - now filling')
    names = fieldnames(VEL);
    for i = 1:length(names)
        if length(VEL.(names{i}))==N            %i.e., only work on timeseries variables
            [newjd,newVEL.(names{i})] = interp1gap(jd,VEL.(names{i}));
        else
            newVEL.(names{i}) = VEL.(names{i}); %copy anything else that isn't a timeseries
        end
    end
    jd = newjd;VEL = newVEL;N = length(jd);
end
% won't actually use crop_method until the conversion to .nc
if isfield('trim_method',metadata),s = rmfield(metadata,'trim_method');end

%move time to center of ensemble
jd = jd + (instmeta.AverageInterval/86400/2);
time = floor(jd); time2 = (jd - floor(jd))*86400000;

%perform time correction if specified
if isfield(metadata,'TimeCorrection')
    newjd = jd + metadata.TimeCorrection/24;
    time = floor(newjd); time2 = (newjd - floor(newjd))*86400000;
    VEL.jd = newjd;
else
    VEL.jd = jd;
end

%write to a netCDF file
% autonan_on
f = netcdf([cdfFile '-raw.cdf'],'clobber');
theFillValue = 1e35;

%put global attributes from User metadata into file
Uatts = fieldnames(metadata);
for i = 1:length(Uatts)
    f.(Uatts{i}) = metadata.(Uatts{i});
end
%now put global attributes from instrument metadata into file
Iatts = fieldnames(instmeta);
for i = 1:length(Iatts)
    f.(Iatts{i}) = instmeta.(Iatts{i});
end
if exist('logmeta','var')
    atts3 = fieldnames(logmeta);
    for i = 1:length(atts3)
        f.(atts3{i}) = logmeta.(atts3{i});
    end
else
end
% MM fix up a couple of little metadata items
f.initial_instrument_height = metadata.transducer_offset_from_bottom;

f.StartTime = instmeta.DeploymentTime;
%fill in data from Nortek Matlab File
%define dimensions
%Time is the record dimension
f('time') = 0;
f('depth') = M;
f('lat') = 1;
f('lon') = 1;
f('matrix') = 3;

%convert from m/s to cm/s
VEL.V1 = VEL.Vvel1*100;
VEL.V2 = VEL.Vvel2*100;
VEL.V3 = VEL.Vvel3*100;
VEL.UserMetadata = metadata;
VEL.InstrumentMetadata = instmeta;
VEL = rmfield(VEL,'Vvel1');VEL = rmfield(VEL,'Vvel2');VEL = rmfield(VEL,'Vvel3');

disp('Writing to netCDF file...')
%fill variable values and define variable attributes

f{'time'} = nclong('time');
f{'time'}(1:N) = time;
f{'time'}.FORTRAN_format = ncchar('F10.2');
f{'time'}.units = ncchar('True Julian Day');
f{'time'}.type = ncchar('UNEVEN');
f{'time'}.epic_code = nclong(624);

f{'time2'} = nclong('time') ;
f{'time2'}(1:N) = time2;
f{'time2'}.FORTRAN_format = ncchar('F10.2');
f{'time2'}.units = ncchar('msec since 0:00 fMT');
f{'time2'}.type = ncchar('UNEVEN');
f{'time2'}.epic_code = nclong(624);

f{'lat'} = nclong('lat'); %% 1 element.
f{'lat'}(:) = metadata.latitude;
f{'lat'}.FORTRAN_format = ncchar('F10.2');
f{'lat'}.units = ncchar('Degrees North');
f{'lat'}.type = ncchar('EVEN');
f{'lat'}.epic_code = nclong(500);
f{'lat'}.name = ncchar('LAT');
f{'lat'}.long_name = ncchar('LATITUDE');
f{'lat'}.generic_name = ncchar('lat');
f{'lat'}.DATUM = ncchar(metadata.LatLonDatum);

f{'lon'} = nclong('lon'); %% 1 element.
f{'lon'}(:) = metadata.longitude;
f{'lon'}.FORTRAN_format = ncchar('f10.4');
f{'lon'}.units = ncchar('Degrees East');
f{'lon'}.type = ncchar('EVEN');
f{'lon'}.epic_code = nclong(502);
f{'lon'}.name = ncchar('LON');
f{'lon'}.long_name = ncchar('LONGITUDE');
f{'lon'}.generic_name = ncchar('lon');
f{'lon'}.DATUM = ncchar(metadata.LatLonDatum);

f{'depth'} = ncfloat('depth');
f{'depth'}(1:M) = VEL.Depths;
f{'depth'}.FORTRAN_format = ncchar('F10.2');
f{'depth'}.units = ncchar('m');
f{'depth'}.type = ncchar('EVEN');
f{'depth'}.epic_code = nclong(3);
f{'depth'}.long_name = ncchar('depth of bin relative to input Mean Water depth (m)');
f{'depth'}.blanking_distance = ncdouble(instmeta.BlankingDistance);
f{'depth'}.bin_size = ncdouble(bin);
f{'depth'}.transducer_offset_from_bottom = metadata.transducer_offset_from_bottom;

f{'bindist'} = ncfloat('depth');
f{'bindist'}(1:M) = VEL.bindist;
f{'bindist'}.long_name = ncchar('Distance of bin off instrument head');
f{'bindist'}.units = ncchar('cm/s');
f{'bindist'}.blanking_distance = ncdouble(instmeta.BlankingDistance);
f{'bindist'}.bin_size = ncdouble(bin);
f{'bindist'}.transducer_offset_from_bottom = metadata.transducer_offset_from_bottom;

f{'VEL1'} = ncfloat('time','depth');
f{'VEL1'}(1:N,:) = VEL.V1;
f{'VEL1'}.FillValue = ncdouble(999);
f{'VEL1'}.units = ncchar('cm/s');
f{'VEL1'}.Type = ncchar('Scalar');
f{'VEL1'}.FillValue_ = theFillValue;

f{'VEL2'} = ncfloat('time','depth');
f{'VEL2'}(1:N,:) = VEL.V2;
f{'VEL2'}.FillValue = ncdouble(999);
f{'VEL2'}.units = ncchar('cm/s');
f{'VEL2'}.Type = ncchar('Scalar');
f{'VEL2'}.FillValue_ = theFillValue;

f{'VEL3'} = ncfloat('time','depth');
f{'VEL3'}(1:N,:) = VEL.V3;
f{'VEL3'}.FillValue = ncdouble(999);
f{'VEL3'}.units = ncchar('cm/s');
f{'VEL3'}.Type = ncchar('Scalar');
f{'VEL3'}.FillValue_ = theFillValue;

switch instmeta.CoordinateSystem
    case 'ENU'
        f{'VEL1'}.long_name = ncchar('Eastward current velocity');
        f{'VEL2'}.long_name = ncchar('Northward current velocity');
        f{'VEL3'}.long_name = ncchar('Vertical current velocity');
    case 'XYZ'
        f{'VEL1'}.long_name = ncchar('Current velocity in X Direction');
        f{'VEL2'}.long_name = ncchar('Current velocity in Y Direction');
        f{'VEL3'}.long_name = ncchar('Current velocity in Z Direction');
    case 'BEAM'
        f{'VEL1'}.long_name = ncchar('Beam 1 current velocity');
        f{'VEL2'}.long_name = ncchar('Beam 2 current velocity');
        f{'VEL3'}.long_name = ncchar('Beam 3 current velocity');
end

f{'AMP1'} = ncfloat('time','depth');
f{'AMP1'}(1:N,:) = VEL.Vamp1;
f{'AMP1'}.FillValue = ncdouble(999);
f{'AMP1'}.units = ncchar('counts');
f{'AMP1'}.long_name = ncchar('Beam 1 echo amplitude');
f{'AMP1'}.Type = ncchar('Scalar');
f{'AMP1'}.FillValue_ = theFillValue;

f{'AMP2'} = ncfloat('time','depth');
f{'AMP2'}(1:N,:) = VEL.Vamp2;
f{'AMP2'}.FillValue = ncdouble(999);
f{'AMP2'}.units = ncchar('counts');
f{'AMP2'}.long_name = ncchar('Beam 2 echo amplitude');
f{'AMP2'}.Type = ncchar('Scalar');
f{'AMP2'}.FillValue_ = theFillValue;

f{'AMP3'} = ncfloat('time','depth');
f{'AMP3'}(1:N,:) = VEL.Vamp3;
f{'AMP3'}.FillValue = ncdouble(999);
f{'AMP3'}.units = ncchar('counts');
f{'AMP3'}.long_name = ncchar('Beam 3 echo amplitude');
f{'AMP3'}.Type = ncchar('Scalar');
f{'AMP3'}.FillValue_ = theFillValue;

f{'Temperature'} = ncfloat('time');
f{'Temperature'}(1:N) = VEL.temperature;
f{'Temperature'}.units = ncchar('Degress C');
f{'Temperature'}.long_name = ncchar('ADCP Transducer Temp');
f{'Temperature'}.FillValue_ = theFillValue;

f{'Pressure'} = ncfloat('time');
f{'Pressure'}(1:N) = VEL.pressure;
f{'Pressure'}.units = ncchar('Dbar');
f{'Pressure'}.long_name = ncchar('Pressure record at instrument head');
f{'Pressure'}.FillValue_ = theFillValue;

f{'Battery'} = ncfloat('time');
f{'Battery'}(1:N) = VEL.battery;
f{'Battery'}.units = ncchar('Volts DC');
f{'Battery'}.long_name = ncchar('Instrument Battery Voltage');
f{'Battery'}.FillValue_ = theFillValue;

f{'Pitch'} = ncfloat('time');
f{'Pitch'}(1:N) = VEL.pitch;
f{'Pitch'}.units = ncchar('Degrees');
f{'Pitch'}.long_name = ncchar('Instrument Pitch');
f{'Pitch'}.FillValue_ = theFillValue;

f{'Roll'} = ncfloat('time');
f{'Roll'}(1:N) = VEL.roll;
f{'Roll'}.units = ncchar('Degrees');
f{'Roll'}.long_name = ncchar('Instrument Roll');
f{'Roll'}.FillValue_ = theFillValue;

f{'Heading'} = ncfloat('time');
f{'Heading'}(1:N) = VEL.heading;
f{'Heading'}.units = ncchar('Degrees');
f{'Heading'}.long_name = ncchar('Instrument Heading');
f{'Heading'}.FillValue_ = theFillValue;

f{'Status'} = ncfloat('time');
f{'Status'}(1:N) = VEL.status;
f{'Status'}.units = ncchar('');
f{'Status'}.long_name = ncchar('Nortek status code');
f{'Status'}.FillValue_ = theFillValue;

f{'Error'} = ncfloat('time');
f{'Error'}(1:N) = VEL.error;
f{'Error'}.units = ncchar('');
f{'Error'}.long_name = ncchar('Nortek error code');
f{'Error'}.FillValue_ = theFillValue;

f{'Soundspeed'} = ncfloat('time');
f{'Soundspeed'}(1:N) = VEL.soundspeed;
f{'Soundspeed'}.units = ncchar('m/s');
f{'Soundspeed'}.long_name = ncchar('Speed of sound used in velocity calculations');
f{'Soundspeed'}.FillValue_ = theFillValue;

if isfield(metadata,'AnalogInput1')
    varname = 'AnalogInput1';f{varname} = ncfloat('time');Aobj = f{varname};
    Aobj(1:N) = VEL.AnaInp1;
    Aobj.units = ncchar('Volts');
    Aobj.sensor_type = metadata.AnalogInput1.sensor_type;
    Aobj.sensor_manufacturer = metadata.AnalogInput1.manufacturer;
    Aobj.serial_number = metadata.AnalogInput1.serial_number;
    Aobj.initial_sensor_height = metadata.AnalogInput1.sensor_height;
    if isfield(metadata.AnalogInput1,'range'),Aobj.range = metadata.AnalogInput1.range;end
    if isfield(metadata.AnalogInput1,'cals'),Aobj.NTUCoefficients = ncdouble(metadata.AnalogInput1.cals.NTUcoef(:));end
    if isfield(metadata.AnalogInput1,'cals')
        if isfield(metadata.AnalogInput1.cals,'SEDcoef')
            f{'AnalogInput1'}.SEDCoefficients = ncdouble(metadata.AnalogInput1.cals.SEDcoef(:));
        end
    end
    Aobj.FillValue_ = theFillValue;
end

if isfield(metadata,'AnalogInput2')
    varname = 'AnalogInput2';f{varname} = ncfloat('time');Aobj = f{varname};
    Aobj(1:N) = VEL.AnaInp2;
    Aobj.units = ncchar('Volts');
    Aobj.sensor_type = metadata.AnalogInput1.sensor_type;
    Aobj.sensor_manufacturer = metadata.AnalogInput1.manufacturer;
    Aobj.serial_number = metadata.AnalogInput1.serial_number;
    Aobj.initial_sensor_height = metadata.AnalogInput1.sensor_height;
    if isfield(metadata.AnalogInput2,'range'),Aobj.range = metadata.AnalogInput1.range;end
    if isfield(metadata.AnalogInput2,'cals'),Aobj.NTUCoefficients = ncdouble(metadata.AnalogInput2.cals.NTUcoef(:));end
    if isfield(metadata.AnalogInput2,'cals')
        if isfield(metadata.AnalogInput2.cals,'SEDcoef')
            f{'AnalogInput1'}.SEDCoefficients = ncdouble(metadata.AnalogInput2.cals.SEDcoef(:));
        end
    end
    Aobj.FillValue_ = theFillValue;
end


f{'TransMatrix'} = ncfloat('matrix','matrix');
f{'TransMatrix'}(:) = instmeta.TransMatrix;
f{'TransMatrix'}.long_name = ncchar('Transformation Matrix for this AWAC');

% autonan_off
close(f)