function VEL = aqdprf2cdf(baseFile,cdfFile,metadata)
%function aqd2cdf('baseFile','cdfFile')
%
%   Function to read in raw Nortek Aquadopp profiler data and create  a cdf file of
%   the raw data. Looks for output Matlab 
%
%   baseFile - base filename for raw Nortek Aquadopp Profiler data (.prf)
%
%   cdfFile - desired basefilename for output netcdf file
%   
%   metaFile - name of metadata file
%
%   this will look for the Nortek header file (.hdr) and log file from
%   deployment (.log) for metadata about the instrument and setup to put
%   into the final netcdf file, so these need to be in the current
%   directory
%
% 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
% first get the metadata in the metaAQD.m file from the user
% if ~exist('metadata','var')
%     run(metaFile);
% 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 = readAQDprfHeader(baseFile);
else
    disp('Need the aquadopp header file (.hdr) for setup information')
    return
end
if ~isempty(strmatch([baseFile '.log'],names));
    logmeta = readAQDlog(baseFile);
else
end

%now load sensor data file for times and sensor data
disp('Loading ASCII Data files')
senfile = strcat(baseFile,'.sen');
SEN = load(senfile);
% VEL.error = SEN(:,7);VEL.status = SEN(:,8);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);
jd = julian([SEN(:,3) SEN(:,1) SEN(:,2) SEN(:,4) SEN(:,5) SEN(:,6)]);
int = (jd(2)-jd(1));

clear SEN;
%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
[N,M] = size(VEL.Vvel1);
% 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 verify the global metadata for standard EPIC and cmg stuff
% everything in metadata and instmeta get written as global attributes
% these also get copied to the .nc file
% metadata = check_gmeta(metadata);
% could be Nan, trap this...
if ~isfield(metadata,'initial_instrument_height') || isnan(metadata.initial_instrument_height),
    metadata.initial_instrument_height = 0;
end
metadata.serial_number = instmeta.AQDHeadSerialNumber; % this is an override with info embedded in the instrument

% update metadata from Aquadopp header to CMG standard so that various
% profilers have the same attribute wording.  Redundant, but necessary
metadata.bin_count = instmeta.AQDNumberOfCells;
metadata.bin_size = instmeta.AQDCellSize./100; % from cm to m
metadata.blanking_distance = instmeta.AQDBlankingDistance; % already in m
%Nortek lists the distance to the center of the first bin as the blanking
%distance plus one cell size
metadata.center_first_bin = metadata.blanking_distance+metadata.bin_size; % in m
metadata.salinity_set_by_user = instmeta.AQDSalinity;
metadata.salinity_set_by_user_units = 'ppt';
metadata.frequency = instmeta.AQDFrequency;
metadata.beam_width = instmeta.AQDBeamWidth;
metadata.beam_pattern = instmeta.AQDBeamPattern;
metadata.beam_angle = instmeta.AQDBeamAngle;
instmeta.AQDHeadRotation = metadata.head_rotation;
metadata = rmfield(metadata,'head_rotation');

%find gaps in the record and fill with NaN's - this is slow, but it works
%and it doesn't change the timebase like interp1gap
%look for small gaps in the record 
hic = find(round(diff(jd)*86400)~=instmeta.AQDProfileInterval);
if length(hic)>1
    disp(['There are ',num2str(length(hic)),' of gaps in the record'])
else
    disp('There are no gaps in the record')
end
% if length(hic)>1 && metadata.fixgaps
%     disp('User instructed to fill gaps with NaNs - now filling'),tic;
%     jdnew = (jd(1):int/1440:jd(end))';
%     ind = find(round(diff(jd)*1440)>round(int));
%     names = fieldnames(VEL);m = length(names);
%     DATA = VEL.(names{1});parsi{1} = size(VEL.(names{1}),2);
%     for i = 2:m,
%         DATA = [DATA VEL.(names{i})];parsi{i} = size(VEL.(names{i}),2);
% %             if length(VEL.(names{i}))==N,data(:,i) = VEL.(names{i});else,continue,end
%     end
%     data= DATA;[N,M] = size(DATA);
%     for i = 1:length(ind)
%         jdgb = jd(ind(i));jdge = jd(ind(i)+1);jdfill = (jdgb:int:jdge)';
%         ind2 = near(jdnew,jdgb,1)-1;
%         jdtemp = [jdnew(1:ind2); jdfill; jd(ind(i)+1:end)];
%         jdnew = jdtemp;
%         
%         fill = (ones(length(jdfill),M)*NaN);
%         dummy = [data(1:ind2,:); fill; DATA(ind(i)+1:end,:)];
%         data = dummy;
%     end
%     count = 1;
%     for i = m,
%         indxs = [count count+parsi{i}];
%         newVEL.(names{i}) = data(:,indxs(1):indxs(2));
% %             if length(VEL.(names{i}))==N,data(:,i) = VEL.(names{i});else,continue,end
%     end
%     jd = jdnew;N = length(jd);
%     ptime = toc/60;
%     disp([num2str(ptime) ' minutes'])
%     fillvals = zeros(length(data),1);fillvals(isnan(data(:,5))) = NaN;
%     oldVEL = VEL;VEL = newVEL;clear newVEL
% end
%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.

if length(hic)>1 && 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
[N,M] = size(VEL.Vamp1);            

% bindist is range to bin from head
disp(metadata.orientation)
fprintf('Center_first_bin = %f\n',metadata.center_first_bin);
fprintf('bin_size = %f\n',metadata.bin_size);
fprintf('bin_count = %f\n',metadata.bin_count);
metadata.bin_count
VEL.bindist = metadata.center_first_bin:metadata.bin_size:(metadata.center_first_bin+((metadata.bin_count-1)*metadata.bin_size));
switch metadata.orientation
    case 'UP'
        disp('User instructed that instrument was pointing UP')
        % depth, or distance below surface, is a positive number below the
        % surface, negative above the surface, for CMG purposes and consistency with ADCP
        VEL.Depths = metadata.WATER_DEPTH - metadata.transducer_offset_from_bottom - VEL.bindist;
        Depth_NOTE = 'uplooking bin depths = water_depth - transducer offset from bottom - bindist';
    case 'DOWN'
        disp('User instructed that instrument was pointing DOWN')
        VEL.Depths = metadata.WATER_DEPTH - metadata.transducer_offset_from_bottom + VEL.bindist;
        Depth_NOTE = 'downlooking bin depths = water_depth - transducer_offset_from_bottom + bindist';
end

%move time to center of ensemble
jd = jd + (instmeta.AQDAverageInterval/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 file
rawname = [cdfFile,'-raw.cdf'];
f = netcdf(rawname,'clobber');

%put global attributes from metadata into file

atts = fieldnames(metadata);
for i = 1:length(atts)
    if isstruct(metadata.(atts{i}))
        continue
    else
        f.(atts{i}) = metadata.(atts{i});
        f.trim_method = []; % won't actually use crop_method until the conversion to .nc
    end
end
Iatts = fieldnames(instmeta);
for i = 1:length(Iatts)
    f.(Iatts{i}) = instmeta.(Iatts{i});
end
if exist('logmeta','var')
    Latts = fieldnames(logmeta);
    for i = 1:length(Latts)
        f.(Latts{i}) = logmeta.(Latts{i});
    end
else
end

% f.SystemFrequency = num2str(instmeta.AQDFrequency);
% f.AQDCoordinateSystem = instmeta.AQDCoordinateSystem;
% f.INST_SER_NUM = instmeta.AQDInstrumentSerialNumber;
%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;

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'} = 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(metadata.LatLonDatum);

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(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 = nclong(instmeta.AQDBlankingDistance);
f{'depth'}.bin_size = ncdouble(metadata.bin_size);
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 from instrument head');
f{'bindist'}.blanking_distance = ncdouble(instmeta.AQDBlankingDistance);
f{'bindist'}.bin_size = ncdouble(metadata.bin_size);
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{'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{'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');

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','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{'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{'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{'Temperature'} = ncfloat('time');
f{'Temperature'}(1:N) = VEL.temperature;
f{'Temperature'}.units = ncchar('Degress 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{'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('Degrees');
f{'Heading'}.long_name = ncchar('Instrument Heading');

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



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

disp(['Finished writing raw data to ' rawname])
close(f)