function aqd2cdf(aqdFile,cdfFile,metadata)
%function aqd2cdf('aqdFile','cdfFile','metaFile')
%
%   Function to read in raw Nortek Aquadopp Single point current meter 
%   data and create  a cdf file of the raw data. 
%
%   aqdFile - base filename for raw Nortek Aquadopp Profiler data (.aqd)
%
%   cdfFile - desired basefilename for output netcdf file
%
%       

if nargin<1
    help(mfilename)
    return
end

% get the metadata about setup and log from deployment
files = dir;
for i = 1:length(files)
    names{i} = files(i).name;
end
if ~isempty(strmatch([aqdFile '.hdr'],names))
    instmeta = readAQDHeader(aqdFile);
else
end
if ~isempty(strmatch([aqdFile '.log'],names))
    logmeta = readAQDlog(aqdFile);
else
end
%now load data file 
datfile = strcat(aqdFile,'.dat');
disp(['Loading ASCII data file ' datfile])
data = load(datfile);
VEL.Vvel(:,1) = data(:,9);VEL.Vvel(:,2) = data(:,10);VEL.Vvel(:,3) = data(:,11);
VEL.Vamp1 = data(:,12);VEL.Vamp2 = data(:,13);VEL.Vamp3 = data(:,14);
VEL.battery = data(:,15);VEL.soundspeed = data(:,16);VEL.heading = data(:,17);
VEL.pitch = data(:,18);VEL.roll = data(:,19);VEL.pressure = data(:,20);
VEL.temperature = data(:,21);VEL.error = data(:,7);VEL.status = data(:,8);
% now get the external datasor data if it was logged, and convert to volts
if isfield(metadata,'AnalogInput1'),VEL.AnaInp1 = data(:,22)*5/65535;end
if isfield(metadata,'AnalogInput2'),VEL.AnaInp2 = data(:,23)*5/65535;end
[N,M] = size(VEL.Vvel);            

jd = julian([data(:,3) data(:,1) data(:,2) data(:,4) data(:,5) data(:,6)]);

%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.
%look for small gaps in the record 
hic = find(round(diff(jd)*86400)~=instmeta.AQDMeasurementInterval);
if length(hic)>1
    disp(['There are ',num2str(length(hic)),' gaps in the record'])
else
    disp('There are no gaps in the record')
end

if 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
            [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

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

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


DEPTH = metadata.WATER_DEPTH - metadata.transducer_offset_from_bottom;
%write to file
autonan_on
f = netcdf([cdfFile,'-raw.cdf'],'clobber');
%put global attributes from metadata into file
% metadata = rmfield(metadata,'AnalogInput1');
atts = fieldnames(metadata);
% for i = 1:length(atts),f.(atts{i}) = metadata.(atts{i});end
for i = 1:length(atts),
    if isstruct(metadata.(atts{i}))
        continue
    elseif ischar(metadata.(atts{i}))
        f.(atts{i}) = ncchar(metadata.(atts{i}));
    elseif isnumeric(metadata.(atts{i}))
        f.(atts{i}) = ncdouble(metadata.(atts{i}));
    end
end

%put global attributes from instrument metadata into file
Iatts = fieldnames(instmeta);
for i = 1:length(Iatts),f.(Iatts{i}) = instmeta.(Iatts{i});end
%put global attributes from logfile into file
if exist('logmeta','var')
    Latts = fieldnames(logmeta);
    for i = 1:length(Latts),f.(atts{i}) = logmeta.(Latts{i});end
end
% f.StartTime = instmeta.clockdeploy;
% f.SystemFrequency = num2str(instmeta.frequency);
% f.AQDCoordinateSystem = instmeta.coordsystem;
% % f.COORD_SYSTEM = VEL.coordsystem;
% f.DiagnosticsInterval = instmeta.diaginterval;
% f.INST_SER_NUM = instmeta.serialno;
%fill in data from Nortek Matlab File
%define dimensions
%Time is the record dimension
f('time') = 0;
f('depth') = 1;
f('lat') = 1;
f('lon') = 1;
f('matrix') = 3;

%convert from m/s to cm/s
VEL.Vvel = VEL.Vvel*100;

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{'VEL1'} = ncfloat('time');
f{'VEL1'}(1:N) = VEL.Vvel(:,1);
f{'VEL1'}.FillValue = ncdouble(999);
f{'VEL1'}.units = ncchar('cm/s');
f{'VEL1'}.Type = ncchar('Scalar');

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

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

switch instmeta.AQDCoordinateSystem
    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');
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{'AMP2'} = ncfloat('time');
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{'AMP3'} = ncfloat('time');
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{'Temperature'} = ncfloat('time');
f{'Temperature'}(1:N) = VEL.temperature;
f{'Temperature'}.units = ncchar('Defress C');
f{'Temperature'}.long_name = ncchar('ADCP Transducer Temp');

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{'Battery'} = ncfloat('time');
f{'Battery'}(1:N) = VEL.battery;
f{'Battery'}.units = ncchar('Volts DC');
f{'Battery'}.long_name = ncchar('Instrument Battery Voltage');

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

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

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

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

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

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{'depth'} = ncfloat('depth');
f{'depth'}(1) = DEPTH;
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 = nclong(instmeta.AQDBlankingDistance);
f{'depth'}.transducer_offset_from_bottom = metadata.transducer_offset_from_bottom;

f{'lon'} = ncfloat('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('NAD83');

f{'lat'} = ncfloat('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('NAD83');

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

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
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
end

autonan_off
close(f)