function VEL = AWAC2nc(CDFFile,outFile,metadata,varargin)
%   function AWAC2nc('CDFFile','OutFile',metadata,[good_ens])
%   
%   Function to take raw AWAC profile data already in netCDF format, 
%   perform rotations, correct pressure data and write to  EPIC netCDF. 
%   Data can be in either beam, earth or XYZ coordinates. Must do aqd2cdf
%   first to get raw data from binary to Matlab, and then translate raw
%   data into a netCDF file. You must have deployment and recovery times
%   entered into the metadata file in order to crop the bad data before and
%   after deployment, otherwise all data will be kept.
%                                                          
%   Required input:
%   'cdfFile' - Raw Netcdf data file base name (.cdf) without extension
%   'OutFile' - base name for output data file - use mooring number
%   'metadata' - metadata from the run script
%   Optional Input
%       good_ens - a vector of the first and last good ensembles to trim the
%       data by. If left empty, the deployment date and recovery date will
%       be used.

% 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
if nargin>3
    good_ens = varargin;
end
autonan_on
cdf = netcdf([CDFFile,'.cdf'],'nowrite');

%first determine if timeseries is not conintuous, an indication in the
%aquadopp that the battery level was too low and it shut down and began 
%intermittent start/stop cycle, as battery power recharged. Although
%data beyond this point is theoretically good, we'll cut it off for the
%Best Basic Version to have a continuous timeseries

time = cdf{'time'}(:);
time2 = cdf{'time2'}(:);
jd = time + (time2/86400000);
deploy = datenum2julian(datenum(cdf.Deployment_date(:),0));
recover = datenum2julian(datenum(cdf.Recovery_date(:),0));
indx = find(jd>=deploy & jd<=recover);
gtime = gregorian(jd);
N = length(time);
n = min(find(diff(gtime(:,6))>0));    
if n<N
    E = n;S=1;
    cropflag = 1;
elseif isfield('metadata','good_ens')
    if isinf(good_ens(2)),good_ens(2) = N;end;
    E = good_ens(2);
    cropflag = 2;
    S = good_ens(1);
else 
    E = indx(end);
    S = indx(1);
    cropflag = 0;
end
names = ncnames(var(cdf));
Histcomment = '';
%read in good data
vel1 = cdf{'VEL1'}(S:E,:);
vel2 = cdf{'VEL2'}(S:E,:);
vel3 = cdf{'VEL3'}(S:E,:);
heading = cdf{'Heading'}(S:E,:);
pitch = cdf{'Pitch'}(S:E,:);
roll = cdf{'Roll'}(S:E,:);
press = cdf{'Pressure'}(S:E,:);
temp = cdf{'Temperature'}(S:E,:);
AMP1 = cdf{'AMP1'}(S:E,:);
AMP2 = cdf{'AMP2'}(S:E,:);
AMP3 = cdf{'AMP3'}(S:E,:);
time = cdf{'time'}(S:E);
time2 = cdf{'time2'}(S:E);
if ~isempty(strmatch('AnalogInput1',names));
    AnaInp1 = cdf{'AnalogInput1'}(S:E);
    nextsens = 1;
end
if ~isempty(strmatch('AnalogInput2',names));
    AnaInp2 = cdf{'AnalogInput2'}(S:E);
    nextsens = nextsens+1;
end

VEL.jd = time + (time2/86400000);
[N,M] = size(vel1);

theFillValue = 1e35;

%get actual water depth from mean pressure
Wdepth = nanmean(press) + cdf.transducer_offset_from_bottom;

% %correct bin depth estimates based on actual water depth
% blank = cdf.BlankingDistance;
% blank2 = cdf.BlankingDistance + cdf.transducer_offset_from_bottom;
% bin = cdf{'depth'}.bin_size(:);
% VEL.bindist = blank+(bin):bin:(bin*(M-1))+blank+(bin);	%depth vector for plotting distance from bottom
% VEL.Depths = fliplr(Wdepth-((bin*(M-1))+blank2+(bin)):bin:Wdepth-(blank2+(bin)));	%depth vector for plotting by depth
sensdepth = Wdepth-cdf.transducer_offset_from_bottom;

% update metadata from Aquadopp header to CMG standard so that various
% profilers have the same attribute wording.  Redundant, but necessary
metadata.bin_count = cdf.NumberOfCells(:);
metadata.bin_size = cdf.CellSize./100; % from cm to m
metadata.blanking_distance = cdf.BlankingDistance; % 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 = cdf.Salinity;
metadata.salinity_set_by_user_units = 'ppt';
metadata.frequency = cdf.Frequency;
metadata.beam_pattern = cdf.BeamPattern;
metadata.beam_angle = cdf.BeamAngle;

%create depth vector
disp('calculating depths and bin distances')
% bindist is range to bin from head
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));
VEL.depths = metadata.WATER_DEPTH - metadata.transducer_offset_from_bottom - VEL.bindist;


%now do rotations from either Beam or XYZ into ENU
disp('Performing transformations to Earth Coordinates...')
%now do rotations from either Beam or XYZ into ENU
%transformation matrix for S/N 
T = cdf{'TransMatrix'}(:);
% it seems that some old aquadopp data has a different transmatrix
% old stuff is integer
if ~rem(T(1,1),floor(T(1,1))), % integer that needs scaling
    fprintf('The TransMatrix has been scaled by T=T/4096\n')
    T = T/4096;   % Scale the transformation matrix correctly to floating point numbers
end

VEL.U = zeros(N,M);VEL.V = zeros(N,M);VEL.W = zeros(N,M);
magvar = cdf.magnetic_variation*pi/180;
%correct heading for magnetic variation
heading = heading + cdf.magnetic_variation(:);

switch cdf.CoordinateSystem(:)
    case 'ENU'
        disp('Data are already in Earth Coordinates, applying magnetic variation')
        VEL.U =  vel1*cos(magvar) + vel2*sin(magvar);
        VEL.V = -vel1*sin(magvar) + vel2*cos(magvar);
        VEL.W = vel3;
    case 'XYZ'
        disp('Data are in XYZ Coordinates, applying transformation and magnetic variation')
        for i = 1:N
            hh = pi*(heading(i)-90)/180;
            pp = pi*pitch(i)/180;
            rr = pi*roll(i)/180;

            H = [cos(hh) sin(hh) 0; -sin(hh) cos(hh) 0; 0 0 1];

            % Make tilt matrix
            P = [cos(pp) -sin(pp)*sin(rr) -cos(rr)*sin(pp);...
                0             cos(rr)          -sin(rr);  ...
                sin(pp) sin(rr)*cos(pp)  cos(pp)*cos(rr)];
            % Make resulting transformation matrix depending on coord system
            R = H*P;
        %transform the data,rotate for magnetic variation - Negative for
        %Westerly variation
            for j = 1:M
                vel = R*[vel1(i,j) vel2(i,j) vel3(i,j)]';
                VEL.U(i,j) = vel(1,1);VEL.V(i,j) = vel(2,1);VEL.W(i,j) = vel(3,1);
            %             VEL.U =  U*cos(magvar) + V*sin(magvar);   old way of rotating
            %             VEL.V = -U*sin(magvar) + V*cos(magvar);
            end
        end
    case 'BEAM'
        disp('Data are in BEAM Coordinates, applying transformation and magnetic variation')
        for i = 1:N
            hh = pi*(heading(i)-90)/180;
            pp = pi*pitch(i)/180;
            rr = pi*roll(i)/180;

            H = [cos(hh) sin(hh) 0; -sin(hh) cos(hh) 0; 0 0 1];

            % Make tilt matrix
            P = [cos(pp) -sin(pp)*sin(rr) -cos(rr)*sin(pp);...
                0             cos(rr)          -sin(rr);  ...
                sin(pp) sin(rr)*cos(pp)  cos(pp)*cos(rr)];
            % Make resulting transformation matrix depending on coord system
            R = H*P*T;
        %transform the data,rotate for magnetic variation - Negative for
        %Westerly variation
            for j = 1:M
                vel = R*[vel1(i,j) vel2(i,j) vel3(i,j)]';
                VEL.U(i,j) = vel(1,1);VEL.V(i,j) = vel(2,1);VEL.W(i,j) = vel(3,1);
            %             VEL.U =  U*cos(magvar) + V*sin(magvar);   old way of rotating
            %             VEL.V = -U*sin(magvar) + V*cos(magvar);
            end
        end
end


minU = min(VEL.U);maxU = max(VEL.U);
minV = min(VEL.V);maxV = max(VEL.V);
minW = min(VEL.W);maxW = max(VEL.W);
%check amplitudes - Nortek lists the noise floor as 20-30 counts for the 
%2 Mhz profiler, 25-35 counts, for the 1Mhz and 30-40 counts for the 600kHz.
%Picking a value in the higher range seems to cut off a lot of data, so
%we'll choose the lower end and researchers can trim more if desired.
freq = num2str(cdf.Frequency(:));
if isfield(metadata,'cutoff_ampl')
    cutoff = metadata.cutoff_ampl;
    Histcomment = ['User selected amplitude cutoff of ',num2str(cutoff),' counts'];
    disp(['User selected amplitude cutoff of ',num2str(cutoff),' counts for screening bad data'])
else
    switch freq
        case '600'
            cutoff = 35;
        case '1000'
            cutoff = 30;
        case '2000'
            cutoff = 25;
    end
    disp(['Screening data using an Amplitude cutoff of ', num2str(cutoff),' counts for this frequency system....'])
end

for j = 1:M
    ind = find((AMP1(:,j)<cutoff | AMP2(:,j)<cutoff | AMP3(:,j)<cutoff)| (isnan(AMP1(:,j))|isnan(AMP2(:,j))|isnan(AMP3(:,j))));
    VEL.U(ind,j) = theFillValue;
    VEL.V(ind,j) = theFillValue;
    VEL.W(ind,j) = theFillValue;
end
%average the amplitudes from the 3 beams
VEL.AGC = (AMP1 + AMP2 + AMP3)./3;

%calculate water level
WL = press + cdf.transducer_offset_from_bottom;

%convert analog input data to NTU's and 
if exist('AnaInp1','var')&&~isempty(strmatch(metadata.AnalogInput1.sensor_type,'OBS3+'))
    opchan1 = doNORTEKoptic(metadata.AnalogInput1,AnaInp1);
    opchan2 = doNORTEKoptic(metadata.AnalogInput2,AnaInp2);
elseif exist('AnaInp1','var')
    opchan1 = doNORTEKoptic(metadata.AnalogInput1,AnaInp1);
elseif exist('AnaInp2','var')
    opchan2 = doNORTEKoptic(metadata.AnalogInput2,AnaInp2);
else
end

if isfield(metadata,'trim_method')
    
    %now fill in bad data with fill values
    %Nortek lists the distance to the center of the first bin as the blanking
    %distance plus the cell size, and to be conservative, we should compare the
    %distance to the beginning of the cell with the water depth.
    blank = cdf{'depth'}.blanking_distance + cdf.transducer_offset_from_bottom;bin = cdf{'depth'}.bin_size(:);
    Dist2 = blank+(bin/2):bin:(bin*(M-1))+blank+(bin/2);	%depth vector for plotting distance from bottom
    for i = 1:N
        lastbin = max(find(Dist2<=WL(i)));
        if lastbin<M
            VEL.U(i,lastbin+1:M) = theFillValue;
            VEL.V(i,lastbin+1:M) = theFillValue;
            VEL.W(i,lastbin+1:M) = theFillValue;
        else
        end
    end
    %now get rid of bins that exist exclusively above the surface
    lastbin = max(find(VEL.bindist<=max(WL)));
    VEL.U = VEL.U(:,1:lastbin);
    VEL.V = VEL.V(:,1:lastbin);
    VEL.W = VEL.W(:,1:lastbin);
    VEL.AGC = VEL.AGC(:,1:lastbin);
    VEL.depths = VEL.depths(1:lastbin);
    VEL.bindist = VEL.bindist(1:lastbin);
    Histcomment = ['Data trimmed at the surface using pressure data' Histcomment];
end

[N,M] = size(VEL.U);
%write to file

f = netcdf([outFile,'-cal.nc'],'clobber');
atts = att(cdf);
copy(atts,f);

%fix a couple of attributes
Histcomment = ['Data rotated into Earth Coordinates by AWAC2nc using a magnetic variation of ',num2str(magvar),...
    '. Amplitude of ',num2str(cutoff),' counts used for data screening' Histcomment];
if cropflag
    history = ['Deployment shortened due to battery failure, data cropped by intermittent time stamp' ':' Histcomment];
    f.history = history;
else
    f.history = Histcomment;
end
f.WATER_DEPTH = ncdouble(Wdepth);
f.WATER_DEPTH_NOTE = ncchar('Water depth determined by mean pressure value');

copy(cdf{'lat'},f,1,1);
copy(cdf{'lon'},f,1,1);

disp('Writing Data to NetCDF....')
%Time is the record dimension
f('time') = 0;
f('depth') = M;
f('lat') = 1;
f('lon') = 1;

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{'time'}.FillValue_ = theFillValue;
f{'time'}.DATUM = 'Time (UTC) in USGS convention: 2440000 = 0000 h on May 23, 1968';

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

f{'depth'} = ncdouble('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 (m)');
f{'depth'}.blanking_distance = cdf.BlankingDistance(:);
f{'depth'}.bin_size = cdf{'depth'}.bin_size(:);
f{'depth'}.transducer_offset_from_bottom = cdf{'depth'}.transducer_offset_from_bottom(:);
f{'depth'}.FillValue_ = theFillValue;
f{'depth'}.NOTE = ncchar('depths referenced to MSL from deployment pressure data');
f{'depth'}.DATUM = 'MSL';

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(cdf.BlankingDistance(:));
f{'bindist'}.bin_size = cdf{'depth'}.bin_size(:);
f{'bindist'}.transducer_offset_from_bottom = cdf{'depth'}.transducer_offset_from_bottom(:);

f{'u_1205'} = ncfloat('time','depth');
f{'u_1205'}(1:N,1:M) = VEL.U;
f{'u_1205'}.name = ncchar('u');
f{'u_1205'}.long_name = ncchar('Eastward Velocity');
f{'u_1205'}.feneric_name = ncchar('u');
f{'u_1205'}.FORTRAN_format = ncchar(' ');
f{'u_1205'}.units = ncchar('cm/s');
f{'u_1205'}.epic_code = nclong(1205);
f{'u_1205'}.sensor_type = cdf.instrument_type(:);
f{'u_1205'}.sensor_depth = sensdepth;
f{'u_1205'}.initial_sensor_height = cdf.initial_instrument_height;
f{'u_1205'}.serial_number = ncchar(cdf.instrument_serial_number);
f{'u_1205'}.minimum = ncfloat(0);
f{'u_1205'}.maximum = ncfloat(0);
f{'u_1205'}.valid_range = ncfloat([-1000 1000]);
f{'u_1205'}.FillValue_ = theFillValue;

f{'v_1206'} = ncfloat('time','depth'); 
f{'v_1206'}(1:N,1:M) = VEL.V;
f{'v_1206'}.name = ncchar('v');
f{'v_1206'}.lonf_name = ncchar('Northward Velocity');
f{'v_1206'}.feneric_name = ncchar('v');
f{'v_1206'}.FORTRAN_format = ncchar(' ');
f{'v_1206'}.units = ncchar('cm/s');
f{'v_1206'}.epic_code = nclong(1206);
f{'v_1206'}.sensor_type = cdf.instrument_type(:);
f{'v_1206'}.sensor_depth = sensdepth;
f{'v_1206'}.initial_sensor_height = cdf.initial_instrument_height;
f{'v_1206'}.serial_number = ncchar(cdf.instrument_serial_number);
f{'v_1206'}.minimum = ncfloat(0);
f{'v_1206'}.maximum = ncfloat(0);
f{'v_1206'}.valid_range = ncfloat([-1000 1000]);
f{'v_1206'}.FillValue_ = theFillValue;

f{'w_1204'} = ncfloat('time','depth'); 
f{'w_1204'}(1:N,1:M) = VEL.W;
f{'w_1204'}.name = ncchar('w');
f{'w_1204'}.lonf_name = ncchar('Vertical Velocity');
f{'w_1204'}.feneric_name = ncchar('w');
f{'w_1204'}.FORTRAN_format = ncchar(' ');
f{'w_1204'}.units = ncchar('cm/s');
f{'w_1204'}.epic_code = nclong(1204);
f{'w_1204'}.sensor_type = cdf.instrument_type(:);
f{'w_1204'}.sensor_depth = sensdepth;
f{'w_1204'}.initial_sensor_height = cdf.initial_instrument_height;
f{'w_1204'}.serial_number = ncchar(cdf.instrument_serial_number);
f{'w_1204'}.minimum = ncfloat(0);
f{'w_1204'}.maximum = ncfloat(0);
f{'w_1204'}.valid_range = ncfloat([-1000 1000]);
f{'w_1204'}.FillValue_ = theFillValue;

f{'P_1'} = ncfloat('time'); 
f{'P_1'}.name = ncchar('P');
f{'P_1'}(1:N) = press;
f{'P_1'}.long_name = ncchar('PRESSURE (DB)          ');
f{'P_1'}.generic_name = ncchar('depth');
f{'P_1'}.FORTRAN_format = ncchar('f10.1');
f{'P_1'}.units = ncchar('dbar');
f{'P_1'}.epic_code = nclong(1);
f{'P_1'}.sensor_depth = sensdepth;
f{'P_1'}.initial_sensor_height = cdf.initial_instrument_height;
f{'P_1'}.minimum = ncfloat(0);
f{'P_1'}.maximum = ncfloat(0);
f{'P_1'}.serial_number = ncchar(cdf.instrument_serial_number);
f{'P_1'}.valid_range = ncfloat([0 10000]);
f{'P_1'}.FillValue_ = 1.0e35;
f{'P_1'}.comment = ncchar('Absolute pressure minus atmospheric pressure at time of sensor zeroing');

f{'AGC_1202'} = ncfloat('time','depth');
f{'AGC_1202'}(1:N,1:M) = VEL.AGC;
f{'AGC_1202'}.name = ncchar('AGC');
f{'AGC_1202'}.long_name = ncchar('Average Echo Intensity (AGC)');
f{'AGC_1202'}.generic_name = ncchar('AGC');
f{'AGC_1202'}.FORTRAN_format = ncchar('F5.1');
f{'AGC_1202'}.units = ncchar('counts');
f{'AGC_1202'}.epic_code = nclong(1202);
f{'AGC_1202'}.sensor_type = cdf.instrument_type(:);
f{'AGC_1202'}.sensor_depth = sensdepth;
f{'AGC_1202'}.initial_sensor_height = cdf.initial_instrument_height;
f{'AGC_1202'}.serial_number = ncchar(cdf.instrument_serial_number);
f{'AGC_1202'}.minimum = ncfloat(0);
f{'AGC_1202'}.maximum = ncfloat(0);
f{'AGC_1202'}.FillValue_ = theFillValue;

f{'T_28'} = ncfloat('time'); 
f{'T_28'}(1:N) = temp;
f{'T_28'}.name = ncchar('T');
f{'T_28'}.long_name = ncchar('TEMPERATURE (C)          ');
f{'T_28'}.generic_name = ncchar('temp');
f{'T_28'}.FORTRAN_format = ncchar('f10.2');
f{'T_28'}.units = ncchar('C');
f{'T_28'}.epic_code = nclong(20);
f{'T_28'}.sensor_depth = ncdouble(sensdepth);
f{'T_28'}.initial_sensor_height = cdf.initial_instrument_height;
f{'T_28'}.minimum = ncfloat(0);
f{'T_28'}.maximum = ncfloat(0);
f{'T_28'}.serial_number = ncchar(cdf.instrument_serial_number);
f{'T_28'}.valid_range = ncdouble([0 100]);
f{'T_28'}.FillValue_ = 1.0e35;

if exist('AnaInp1','var')
    varname = 'NEP1_56';f{varname} = ncfloat('time');Aobj = f{varname};
    Aobj(1:N) = opchan1(:,1);
    atts = att(cdf{'AnalogInput1'});
    Aobj.long_name = ncchar('Nephylometer');
    for i = 1:length(atts),copy(atts{i},Aobj),end;
    if ~isempty(metadata.AnalogInput1.cals.NTUcoef)
        Aobj.units = ncchar('NTU');
    else
        Aobj.units = ncchar('volts');
    end
    Aobj.sensor_depth = Wdepth - Aobj.initial_sensor_height;
    if size(opchan1,2)>1                  %i.e., we have sed conc data as well
        varname = 'Sed1_981';f{varname} = ncfloat('time');Aobj = f{varname};
        Aobj(1:N) = opchan1(:,2);
        Aobj.long_name = ncchar('Sediment Concentration');
        for i = 1:length(atts),copy(atts{i},Aobj),end;
        Aobj.units = ncchar(metadata.AnalogInput1.cals.SEDcoefUnits);
        Aobj.sensor_depth = Wdepth - Aobj.initial_sensor_height;
    end
end
if exist('AnaInp2','var')
    varname = 'NEP2_56';f{varname} = ncfloat('time');Aobj = f{varname};
    Aobj(1:N) = opchan2(:,1);
    Aobj.long_name = ncchar('Nephylometer');
    atts = att(cdf{'AnalogInput2'});
    for i = 1:length(atts),copy(atts{i},Aobj),end;
    if ~isempty(metadata.AnalogInput2.cals.NTUcoef)
        Aobj.units = ncchar('NTU');
    else
        Aobj.units = ncchar('volts');
    end
    Aobj.sensor_depth = Wdepth - Aobj.initial_sensor_height;
    if size(opchan2,2)>1                  %i.e., we have sed conc data as well
        varname = 'Sed2_981';f{varname} = ncfloat('time');Aobj = f{varname};
        Aobj(1:N) = opchan2(:,2);
        Aobj.long_name = ncchar('Sediment Concentration');
        for i = 1:length(atts),copy(atts{i},Aobj),end;
        Aobj.units = ncchar(metadata.AnalogInput2.cals.SEDcoefUnits);
        Aobj.sensor_depth = Wdepth - Aobj.initial_sensor_height;
    end
end

add_fillvalues(f,theFillValue)
% add_minmaxvalues(f)
autonan_off

disp(['Data written to ' outFile '-cal.nc'])
close(cdf)
close(f)