function aqd2nc(cdfFile,OutFile,metadata)
%   function aqd2nc('cdfFile','OutFile')
%   
%   Function to take raw aquadopp 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. Once data are in netCDF, you can plot to
%   determine when the first and last good ensembles are.
%                                                          
%   Required input:
%   'cdfFile' - Raw Netcdf data file base name (.cdf)
%   'OutFile' - base name for output data file - use mooring number
%   good_ens - vector of good ensembles [first_good last_good] based on in
%               and out of water times
if nargin<1
    help(mfilename)
    return
end
autonan_on
cdf = netcdf([cdfFile,'.cdf']);
% cdf = autonan(cdf,1);
%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);
start = datenum2julian(datenum(cdf.Deployment_date(:)));
finish = datenum2julian(datenum(cdf.Recovery_date(:)));
ind = find(jd>=start&jd<=finish);
gtime = gregorian(jd);
N = ind(end);
n = (find(diff(gtime(:,6))>0, 1 ) | find(diff(gtime(:,5))>0, 1 ));    
if n<N
    E = n;
    cropflag = 1;
else
    E = N;
    cropflag = 0;
end
S = ind(1);
names = ncnames(var(cdf));

%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);
jd = time + (time2/86400000);
[N,M] = size(vel1);

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
%convert analog input data to NTU's and 
if exist('AnaInp1','var')
    if isfield(AnaInp1,'cals')
        disp('User supplied optical calibration parameters for Channel 1 - performing optical conversions')
        opchan1 = doNORTEKoptic(metadata.AnalogInput1,AnaInp1);
    else 
        opchan1 = AnaInp1;
    end
end
if exist('AnaInp2','var')
    if isfield(AnaInp1,'cals')
        disp('User supplied optical calibration parameters for Channel 2 - performing optical conversions')
        opchan2 = doNORTEKoptic(metadata.AnalogInput1,AnaInp2);
    else 
        opchan2 = AnaInp2;
    end
end


theFillValue = 1e35;

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

%now do rotations from either Beam or XYZ into ENU
disp('Performing transformations to Earth Coordinates...')
%transformation matrix for S/N 
T = cdf{'TransMatrix'}(:);
T = T/4096;   % Scale the transformation matrix correctly to floating point numbers

VEL.U = zeros(N,1);VEL.V = zeros(N,1);VEL.W = zeros(N,1);
magvar = cdf.magnetic_variation*pi/180;

switch cdf.AQDCoordinateSystem(:)
    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'
        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;
        end
        %transform the data,rotate for magnetic variation - Negative for
        %Westerly variation
            vel = R*[vel1 vel2 vel3]';
            U = vel(1,:);V = vel(2,:);
            VEL.W = vel(3,:);
            VEL.U =  U*cos(magvar) + V*sin(magvar);
            VEL.V = -U*sin(magvar) + V*cos(magvar);
    case 'BEAM'
        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;
        end
        %transform the data,rotate for magnetic variation - Negative for
        %Westerly variation
            vel = R*[vel1 vel2 vel3]';
            U = vel(1,:);V = vel(2,:);
            VEL.W = vel(3,:);
            VEL.U =  U*cos(magvar) + V*sin(magvar);
            VEL.V = -U*sin(magvar) + V*cos(magvar);
end


%check amplitudes - Nortek lists the noise floor as 20-30 counts for the 
%2 Mhz profiler - 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.

disp('Screening data by predetermined Amplitude cutoff for this frequency system....')
if isfield(metadata,'amplitude_cutoff')
    cutoff = metadata.amplitude_cutoff;
else
    cutoff = 25;
end
    
ind = find(AMP1<cutoff | AMP2<cutoff | AMP3<cutoff);
VEL.U(ind) = theFillValue;
VEL.V(ind) = theFillValue;
VEL.W(ind) = theFillValue;

%average the amplitudes from the 3 beams
AGC = (AMP1 + AMP2 + AMP3)./3;

%calculate water level, listed as hght_18 in epic
WL = press + cdf.transducer_offset_from_bottom;

%now crop data based on water level data
%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.

[N] = length(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 aqd2nc. Amplitude of ',num2str(cutoff),' counts used for data screening'];
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_source = ncchar('Water depth determined by mean pressure value');
f.WATER_DEPTH_datum = ncchar('MSL');
sensdepth = cdf{'depth'}(:);

copy(cdf{'lat'},f,1,1);
copy(cdf{'lon'},f,1,1);
copy(cdf{'depth'},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{'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{'u_1205'} = ncfloat('time');
f{'u_1205'}(1:N) = 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'}.serial_number = ncchar(cdf.AQDHeadSerialNumber);
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'); 
f{'v_1206'}(1:N) = 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'}.serial_number = ncchar(cdf.AQDHeadSerialNumber);
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'); 
f{'w_1204'}(1:N) = 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'}.serial_number = ncchar(cdf.AQDHeadSerialNumber);
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{'hght_18'} = ncfloat('time');
% f{'hght_18'}(1:N) = WL;
% f{'hght_18'}.name = ncchar('hght');
% f{'hght_18'}.long_name = ncchar('height of sea surface from pressure data');
% f{'hght_18'}.generic_name = ncchar('height');
% f{'hght_18'}.FORTRAN_format = ncchar('f10.2');
% f{'hght_18'}.units = ncchar('m');
% f{'hght_18'}.epic_code = nclong(18);
% f{'hght_18'}.sensor_depth = sensdepth;
% f{'hght_18'}.minimum = ncfloat(0);
% f{'hght_18'}.maximum = ncfloat(0);
% f{'hght_18'}.serial_number = ncchar(cdf.AQDHeadSerialNumber);
% f{'hght_18'}.valid_range = ncfloat([0 1000]);
% f{'hght_18'}.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'}.minimum = ncfloat(0);
f{'P_1'}.maximum = ncfloat(0);
f{'P_1'}.serial_number = ncchar(cdf.AQDHeadSerialNumber);
f{'P_1'}.valid_range = ncfloat([0 10000]);
f{'P_1'}.FillValue_ = theFillValue;
f{'P_1'}.comment = ncchar('Absolute pressure minus atmospheric pressure at time of sensor zeroing');

f{'AGC_1202'} = ncfloat('time');
f{'AGC_1202'}(1:N) = 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'}.serial_number = ncchar(cdf.AQDHeadSerialNumber);
f{'AGC_1202'}.minimum = ncfloat(0);
f{'AGC_1202'}.maximum = ncfloat(0);
f{'AGC_1202'}.FillValue_ = theFillValue;

% if ~isempty(cdf{'AnalogInput1'})
%     copy(cdf{'AnalogInput1'},f,1,1)
% end
% 
% if ~isempty(cdf{'AnalogInput2'})
%     copy(cdf{'AnalogInput2'},f,1,1)
% end

if exist('AnaInp1','var')
    varname = 'NEP1_56';f{varname} = ncfloat('time');Aobj = f{varname};
    Aobj(1:N) = opchan1(:,1);
    atts = att(cdf{'AnalogInput1'});
    for i = 1:length(atts),copy(atts{i},Aobj),end;
%     Aobj.long_name = ncchar('Nephylometer');
    if isfield(metadata.AnalogInput1,'cals')
        if ~isempty(metadata.AnalogInput1.cals.NTUcoef)
            Aobj.units = ncchar('NTU');
        else
            Aobj.units = ncchar('volts');
        end
    end
    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);
    end
end
if exist('AnaInp2','var')
    varname = 'NEP2_56';f{varname} = ncfloat('time');Aobj = f{varname};
    Aobj(1:N) = opchan2(:,1);
    atts = att(cdf{'AnalogInput2'});
%     Aobj.long_name = ncchar('Nephylometer');
    for i = 1:length(atts),copy(atts{i},Aobj),end;
    if isfield(metadata.AnalogInput1,'cals')
        if ~isempty(metadata.AnalogInput1.cals.NTUcoef)
            Aobj.units = ncchar('NTU');
        else
            Aobj.units = ncchar('volts');
        end
    end
    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);
    end
end
add_minmaxvalues(f)
autonan_off

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