function [adcp,config,units]=adcpRDItoMat(ifile,meta)
%---------------------------------------------------------------------
% General:  Generate .mat files with instrument data given an rdi file file
%           in database
% Input:    ifile [String] - rdi file name
%           meta structure
% Output:   adcp
%           config
%           units
% Description:
%   Read ADCP data from rdi file "ifile" and 
%   1) save primary variables into a MATLAB structure
%   For SI & SoG, corrects for magnetic deviation, so U anv V are relative
%   to true north
% Modify Date/Author:  RKD 10/07 for RDI Workhorse
%       Reyna Jenkyns, June 5, 2009
%           modified for Long Ranger and details added
%       20100209, R.Jenkyns
%       - added rotations to Earth coords
%       20100427, R.Jenkyns
%       - remove magnetic u,v from structure, add backscatter for DMAS-3969
%---------------------------------------------------------------------
  
hdg=meta.deviceHeading;
if isnan(hdg), hdg=0; end    

% read binary .rdi into matlab, process and save in mat structure
[adcp,config,units]=vrdadcp(ifile,[1 inf]);% read entire binary file into structures

% determine whether adcp is fixed position
compassRange=iqr(adcp.compassHeading); % robust measure of spread, not sensitive to outliers
if compassRange>5, fixed=0;
else fixed=1;
end    

%    mtime=adcp.mtime; % the matlab times for each ensemble
if strcmp(config.coordSys,'Earth')
    if strcmp(config.sensorSrc_EZ(3),'1') && strcmp(config.sensorAvail(3),'1')% internal compass used to define true north
        if fixed
            % rotate to true Earth coords
            hdg=hdg*pi/180; % heading in radians
            ch=cos(hdg); sh=sin(hdg);
            M=[ch sh;-sh ch];
            mh=nanmean(adcp.compassHeading)*pi/180;
            mch=cos(mh); msh=sin(mh);
            Minv=[mch -msh;msh mch];
            [m,n,p]=size(adcp.velocity);
            v=reshape(adcp.velocity,4,n*p);
            X=M*Minv*v(1:2,:);
        else
            adcp.uMagnetic=squeeze(adcp.velocity(1,:,:)); % east velocity relative to magnetic North
            adcp.vMagnetic=squeeze(adcp.velocity(2,:,:)); % north velocity relative to magnetic north
            % rotate to true Earth coords
            % get magnetic deviation at this location
            [b,magdev]=igrf11syn(adcp.time(1),-meta.depth/1000,meta.lat,meta.lon);
            magdev=magdev*pi/180; % magnetic declination in radians
            ch=cos(magdev); sh=sin(magdev);
            M=[ch -sh;sh ch]; % for ccw rotation
            [m,n,p]=size(adcp.velocity);
            v=reshape(adcp.velocity,4,n*p);
            X=M*v(1:2,:);
        end    
        X=reshape(X,[2 n p]);
        adcp.u=squeeze(X(1,:,:));
        adcp.v=squeeze(X(2,:,:)); 
    elseif strcmp(config.sensorSrc_EZ(3),'0') % EH command used to get true ENU 
        adcp.u=squeeze(adcp.velocity(1,:,:)); 
        adcp.v=squeeze(adcp.velocity(2,:,:)); 
    end    
    adcp.w=squeeze(adcp.velocity(3,:,:));
    adcp.velocityError=squeeze(adcp.velocity(4,:,:));
elseif strcmp(config.coordSys,'Beam')    
    % transforming from beam to orthogonal
    if strcmp(config.beamPattern,'Convex')
        c=1;
    else
        c=-1;
    end
    if config.beamAngle==20
        a=1.4619; b=0.2660; d=1.0337;
    else %beamAngle=30
        a=1.0; b=0.2887; d=.7071;
    end    
    A=[c*a -c*a 0 0; 0 0 -c*a c*a; b b b b; d d -d -d];
    [m,n,p]=size(adcp.velocity);
    v=reshape(adcp.velocity,4,n*p);
    X=A*v;
    % three beam
    % HPR matrices
    hdg=hdg*pi/180; % heading in radians
    ch=cos(hdg);
    sh=sin(hdg);
    M1=[ch sh 0; -sh ch 0; 0 0 1];
    
    if strcmp(config.sensorAvail(4:5),'11')
        pmean=nanmean(adcp.pitch)*pi/180;
        rmean=nanmean(adcp.roll)*pi/180;  
        cp=cos(pmean); sp=sin(pmean);
        cr=cos(rmean); sr=sin(rmean);
        M2=[1 0 0;0 cp -sp;0 sp cp];
        M3=[cr 0 sr;0 1 0;-sr 0 cr];
    else
        M2=eye(3); M3=eye(3);
    end;    
    M=M1*M2*M3; 
    if strcmp(config.orient,'Up') %negate 1st & 3rd column
        M(1:3,1)=-M(1:3,1);
        M(1:3,3)=-M(1:3,3);
    end    
    Y=M*X(1:3,:);
    X=reshape([Y;X(4,:)],[m n p]); 
    adcp.u=squeeze(X(1,:,:));
    adcp.v=squeeze(X(2,:,:)); 
    adcp.w=squeeze(X(3,:,:));
    adcp.velocityError=squeeze(X(4,:,:));
end    
    
% BACKSCATTER    
%Ainslie M. A., McColm J. G., "A simplified formula for
% viscous and chemical absorption in sea water", Journal of the Acoustical
% Society of America, 103(3), 1671-1672, 1998.
T=nanmean(adcp.temperature);
S=nanmean(adcp.salinity);
D=nanmean(adcp.depth); %depth,metres
f=config.freq;
f2=42*exp(T/17);
alpha2=f;
for i=1:length(f)
    alpha2(i)=0.52*(1+T/43)*(S/35)*f2*f(i)^2/(f2^2+f(i)^2)*exp(-D/6000)+0.00049*f(i)^2*exp(-T/27-D/17000); %db/km
end
alpha2=alpha2/1000; % db/m 
BS=adcp.intens*0.45; %convert count to db, p142 manual
% 2 way beam spreading and absorption
lr=2*(20*log10((adcp.range)') + alpha2*adcp.range')*ones(1,size(BS,3));
[m,n]=size(lr);
LR=zeros(4,m,n);
for i=1:4, 
    LR(i,:,:)=lr; % make this a matrix (bins,1-4,time)
end
adcp.backscatter=BS+LR;  % back scatter beams 1-4 range attenuation correction = dB+20LogR
units.backscatter='dB';
config.soundAbsorptionCoefficient=alpha2;
units.soundAbsorptionCoefficient='dB/m';



