function [mdata, ncdata] = load_netcdf(filepath,varargin)
%LOAD_GETNC  Load CODAS netCDF file, mimicing load_getmat.m
%
%   Load the data from a short or long form netCDF file generated by
%   CODAS. Make a data structure as similar as possible to that from
%   the old load_getmat.m function.
%
% SYNTAX:
%    [mdata, ncdata] = load_netcdf(filepath)  % to get everything
%
%    mdata = load_netcdf(filepath)  % to just get a structure like
%                                  % the old load_getmat.m
%
%    [mdata, ncdata] = load_netcdf(filepath,'all') % to calculate some
%                      % (possibly useful) additional parameters
%
%    ncdata = load_netcdf(filepath,'nocodas') % do not throw
%                      % a warning if the file is not from CODAS;
%                      % do not create a fake load_getmat structure
%
%    The last syntax is actually a general-purpose netCDF importer
%    for simple format netCDFs without groups other than root:
%
%
% INPUT:
%    filepath    path and filename of the netcdf file
%
% OUTPUTS:
%    mdata       A data structure resembling that produced by load_getmat
%                See below for the meanings of the fields.
%    ncdata      A structure containing the entire netCDF dataset including
%                global attributes, variables, and variable attributes.
%                Data from each variable in the netCDF file is put in field
%                with the name of the variable (possibly cleaned up to fit
%                Matlab rules). Variable attributes are put in structures
%                with the data field names appended by '_att'.
%
%
% NOTES:
%    The short netCDF files do not contain all of the data that are in the
%    collected .mat files. Missing are
%        iblkprf
%        swcor
%        umeas_bt, vmeas_bt, and depth_bt
%        umean, vmean
%        w, wmean
%        e
%        head_misalign, scale factor
%        last_temp
%    In addition, edited data have already been removed from the netCDF
%    files, so u is the equivalent of u.*nanmask. Same for v, umeas, and
%    vmeas
%
%    The long form netCDF files have the above quantities, aside from
%    iblkfrof and umean, vmean, and wmean. Those last 3 are calculated,
%    but do not agree precisely with the load_getmat results.
%
%    Neither the short nor long netCDF files have the single-ping
%    statistics that can be returned by load_getmat(file,'all'). Those
%    are:
%        amp1, amp2, amp3, amp4   (individual beam signal return amplitude or rssi)
%        rvar_uu   variance of singleping residual (u)
%        rvar_uv   covariance of singleping residual (u,v)
%        rvar_vv   variance of singleping residual (v)
%        rvar_ff   variance of singleping residual (fwd)
%        rvar_fp   covariance of singleping residual (fwd,port)
%        rvar_pp   variance of singleping residual (port)
%
%        tsvar_uu   variance of singleping time series (u)
%        tsvar_uv   covariance of singleping time series (u,v)
%        tsvar_vv   variance of singleping time series (v)
%        tsvar_ff   variance of singleping time series (fwd)
%        tsvar_fp   covariance of singleping time series (fwd,port)
%        tsvar_pp   variance of singleping time series (port)
%        tsvar_n    number of points in singleping time series
%                   (when reading a long-form netCDF, tsvar_n
%                   is the number of pings in the ensemble rather
%                   than the number of pings actually averaged)
%
% Fields in the mdata output (from load_getmat.m)
% -------------------------------------------------------------------
%  fieldnames   size        meaning:
%
%    dday     1   x nprofs   zero-based decimal day (end of profile)
%    time     6   x nprofs   year,month,day,hour,minute,second (end of profile)
% iblkprf     2   x nprofs   database block and profile number
%     lon     1   x nprofs   longitude (from smoothed navigation; end of profile)
%     lat     1   x nprofs   latitude  (from smoothed navigation; end of profile)
% heading     1   x nprofs   average primary heading (usually gyro)
% tr_temp     1   x nprofs   average temperature
% last_temp   1   x nprofs   last temperature
% last_good_bin 1 x nprofs   last good bin (-1 is bad profile)
%     amp   nbins x nprofs   signal return amplitude or rssi (averaged)
%    amp1   nbins x nprofs   signal return amplitude or rssi (beam 1)
%    amp2   nbins x nprofs   signal return amplitude or rssi (beam 2)
%    amp3   nbins x nprofs   signal return amplitude or rssi (beam 3)
%    amp4   nbins x nprofs   signal return amplitude or rssi (beam 4)
%   swcor   nbins x nprofs   spectral width (nb150) or correlation (os or wh)
%   umeas   nbins x nprofs   measured  velocity (positive east, m/s)
%   umeas_bt  1   x nprofs   measured  velocity (positive east, m/s) from bottom track
%   uship     1   x nprofs   ship      velocity (positive east, m/s)
%       u   nbins x nprofs   ocean     velocity (positive east, m/s)
%   umean     1   x nprofs   ref layer velocity (positive east, m/s)
%   vmeas   nbins x nprofs   measured  velocity (positive north, m/s)
%   vmeas_bt  1   x nprofs   measured  velocity (positive north, m/s) from bottom track
%   vship     1   x nprofs   ship      velocity (positive north, m/s)
%       v   nbins x nprofs   ocean     velocity (positive north, m/s)
%   vmean     1   x nprofs   ref layer velocity (positive north, m/s)
%       w   nbins x nprofs   vertical  velocity (up   , m/s)
%   wmean     1   x nprofs   mean vertical velocity for each profile
%       e   nbins x nprofs   error velocity (m/s)
%   depth   nbins x nprofs   center depth for each bin and profile, positive down
%      pg   nbins x nprofs   percent good
%   pflag   nbins x nprofs   0 is unflagged, otherwise there is a code.  see www docs
% nanmask   nbins x nprofs   suitable for dot-multiply with velocity
%
%
% Steven Howell, 21 August 2024
% Tested with Matlab R2024a only.


doall = false;
iscodas = true;

if ~isempty(varargin)
    for ii=1:length(varargin)
        if ischar(varargin{ii}) || isstring(varargin{ii})
            if strcmpi(string(varargin{ii}),"all")
                doall = true;
            elseif strcmpi(string(varargin{ii}),"nocodas")
                iscodas = false;
                doall = false;
            end
        end
    end
end

% Check for new enough releases of Matlab
% Obvious issues for older versions are the datetime functions
% and the netCDF routines

if isMATLABReleaseOlderThan("R2014b")
    if isMATLABReleaseOlderThan("2008b")
        error("Matlab version predates netCDF routines! At least R2008b is required")
        % No netCDF capability!
    else
        error("Your version of Matlab is too old! (earlier than R2014b)")
        % But note that it wouldn't be too hard to rewrite around the
        % lack of the datetime function
    end
end

ncdata = loadnc(filepath);

idx = strcmp('software',{ncdata.globalatts(:).name});
if ~any(idx) && strcmp('pycurrents',ncdata.globalatts(idx).value)~=0 && iscodas
    warn("This does not appear to be a CODAS file!")
    iscodas = false;
    doall = false;
else
    iscodas = true;
end

if isfield(ncdata,'heading')
    isshortnc = true;
else
    isshortnc = false;
end

% Now make a version that looks like get_loadmat.m
mdata = struct();
if iscodas && isshortnc
    mdata.dday = ncdata.time';
    idx = find(strcmp('units',{ncdata.time_att(:).name}));
    reftime = datetime(extractAfter(ncdata.time_att(idx).value,'days since '));
    mdata.time = datevec(reftime+days(ncdata.time))';
    mdata.iblkprf = NaN([2,size(mdata.dday,2)]);
    mdata.lon = missing_to_nan(ncdata.lon, ncdata.lon_att)';
    mdata.lon = mod(mdata.lon+360,360); % .mat files seem to have 0<=lon<=360
    mdata.lat = missing_to_nan(ncdata.lat, ncdata.lat_att)';
    mdata.heading = missing_to_nan(ncdata.heading, ncdata.heading_att)';
    mdata.amp = missing_to_nan(ncdata.amp, ncdata.amp_att);
    mdata.swcor = NaN(size(mdata.amp));
    mdata.umeas = NaN(size(mdata.amp)); % Filled in later
    mdata.umeas_bt = NaN(size(mdata.lon));
    mdata.uship = missing_to_nan(ncdata.uship, ncdata.uship_att)';
    mdata.umean = NaN(size(mdata.uship));
    mdata.u = missing_to_nan(ncdata.u, ncdata.u_att);
    mdata.umeas = mdata.u - repmat(mdata.uship,[size(mdata.umeas,1),1]);
    mdata.vmeas = NaN(size(mdata.amp)); % Filled in later
    mdata.vmeas_bt = NaN(size(mdata.lon));
    mdata.vship = missing_to_nan(ncdata.vship, ncdata.vship_att)';
    mdata.vmean = NaN(size(mdata.vship));
    mdata.v = missing_to_nan(ncdata.v, ncdata.v_att);
    mdata.vmeas = mdata.v - repmat(mdata.vship,[size(mdata.vmeas,1),1]);
    mdata.w = NaN(size(mdata.v));
    mdata.wmean = NaN(size(mdata.dday));
    mdata.e = NaN(size(mdata.amp));
    mdata.depth = missing_to_nan(ncdata.depth, ncdata.depth_att);
    mdata.pg = missing_to_nan(ncdata.pg, ncdata.pg_att);
    mdata.pflag = missing_to_nan(ncdata.pflag, ncdata.pflag_att);
    mdata.nanmask = zeros(size(mdata.u)).*mdata.u + 1;
    mdata.head_misalign = NaN(size(mdata.dday));
    mdata.scale_factor = NaN(size(mdata.dday));
    mdata.last_good_bin = sum(mdata.pflag<4,1); % Seems to work
    mdata.tr_temp = missing_to_nan(ncdata.tr_temp, ncdata.tr_temp_att)';
    mdata.last_temp = NaN(size(mdata.dday));
    mdata.depth_bt = NaN(size(mdata.dday));
elseif iscodas && ~isshortnc
    mdata.dday = ncdata.time';
    idx = find(strcmp('units',{ncdata.time_att(:).name}));
    reftime = datetime(extractAfter(ncdata.time_att(idx).value,'days since '));
    mdata.time = datevec(reftime+days(ncdata.time))';
    mdata.iblkprf = NaN([2,size(mdata.dday,2)]);
    mdata.lon = missing_to_nan(ncdata.lon, ncdata.lon_att)';
    mdata.lon = mod(mdata.lon+360,360); % .mat files seem to have 0<=lon<=360
    mdata.lat = missing_to_nan(ncdata.lat, ncdata.lat_att)';
    mdata.heading = missing_to_nan(ncdata.ANCIL1_mn_heading, ncdata.ANCIL1_mn_heading_att)';
    mdata.amp = missing_to_nan(ncdata.amp_sound_scat, ncdata.amp_sound_scat_att);
    mdata.swcor = missing_to_nan(ncdata.spectral_width, ncdata.spectral_width_att);
    mdata.umeas = missing_to_nan(ncdata.u, ncdata.u_att);
    mdata.umeas_bt = missing_to_nan(ncdata.BT_u, ncdata.BT_u_att)';
    mdata.uship = missing_to_nan(ncdata.ACCESS_U_ship_absolute, ncdata.ACCESS_U_ship_absolute_att)';
    mdata.umean = NaN(size(mdata.uship));
    mdata.u = mdata.umeas + repmat(mdata.uship,[size(mdata.umeas,1),1]);
    mdata.vmeas = missing_to_nan(ncdata.v, ncdata.v_att);
    mdata.vmeas_bt = missing_to_nan(ncdata.BT_v, ncdata.BT_v_att)';
    mdata.vship = missing_to_nan(ncdata.ACCESS_V_ship_absolute, ncdata.ACCESS_V_ship_absolute_att)';
    mdata.vmean = NaN(size(mdata.vship));
    mdata.v = mdata.vmeas + repmat(mdata.vship,[size(mdata.vmeas,1),1]);
    mdata.w = missing_to_nan(ncdata.w, ncdata.w_att);
    mdata.wmean = NaN(size(mdata.dday));
    mdata.e = missing_to_nan(ncdata.error_vel, ncdata.error_vel_att);
    mdata.depth = missing_to_nan(ncdata.depth, ncdata.depth_att);
    mdata.pg = missing_to_nan(ncdata.percent_good, ncdata.percent_good_att);
    mdata.pflag = missing_to_nan(ncdata.profile_flags, ncdata.profile_flags_att);
    mdata.nanmask = double(mdata.pflag==0);
    mdata.nanmask(mdata.nanmask==0) = NaN;
    mdata.head_misalign = missing_to_nan(ncdata.ANCIL2_watrk_hd_misalign, ncdata.ANCIL2_watrk_hd_misalign_att)';
    mdata.scale_factor = missing_to_nan(ncdata.ANCIL2_watrk_scale_factor, ncdata.ANCIL2_watrk_scale_factor_att)';
    mdata.last_good_bin = missing_to_nan(ncdata.ACCESS_last_good_bin, ncdata.ACCESS_last_good_bin_att)';
    mdata.tr_temp = missing_to_nan(ncdata.ANCIL1_tr_temp, ncdata.ANCIL1_tr_temp_att)';
    mdata.last_temp = missing_to_nan(ncdata.ANCIL2_last_temp, ncdata.ANCIL2_last_temp_att)';
    mdata.depth_bt = missing_to_nan(ncdata.BT_depth, ncdata.BT_depth_att)';
    meanmask = NaN(size(mdata.nanmask));
    for jj=1:size(meanmask,2)
        meanmask(ncdata.CONFIG1_top_ref_bin(jj):ncdata.CONFIG1_bot_ref_bin(jj),jj)=1;
    end
    mdata.umean = mean(mdata.u.*meanmask.*mdata.nanmask,1,"omitnan");
    mdata.vmean = mean(mdata.v.*meanmask.*mdata.nanmask,1,"omitnan");
    mdata.wmean = mean(mdata.w.*meanmask.*mdata.nanmask,1,"omitnan");
elseif ~iscodas
    mdata = ncdata;
end
if doall
    mdata.mps = abs(mdata.uship + 1i*mdata.vship);
    %% rotate into forward and port (useful)
    % measured
    rads = (90-mdata.heading(:)')*pi/180;
    onesdown = ones(size(mdata.umeas(:,1)));

    fp = (mdata.umeas + 1i*mdata.vmeas).*exp(-1i*onesdown*rads);
    mdata.fmeas = real(fp);
    mdata.pmeas = imag(fp);
    % ocean
    fp =    (mdata.u + 1i*mdata.v).*exp(-1i*onesdown*rads);
    mdata.fvel  = real(fp);
    mdata.pvel  = imag(fp);
    if isshortnc
        mdata.tsvar_n = Nan(size(mdata.dday));
    else
        mdata.tsvar_n = missing_to_nan(ncdata.ANCIL1_pgs_sample, ncdata.ANCIL1_pgs_sample_att)';
    end
end

end

function ncdata = loadnc(filepath)
    % A simple netCDF loader. Gets all global attributes, variables in the
    % root goup, and the attributes of those variables.

    ncdata = struct();
    ncid = netcdf.open(filepath);

    % Check whether there are child groups.
    childGrps = netcdf.inqGrps(ncid);
    if ~isempty(childGrps)
        warn("This function only loads variables in the root group.")
    end

    [~, ~, numglobalatts, ~] = netcdf.inq(ncid);
    ncdata.globalatts = struct();
    for ii=0:(numglobalatts-1)
        ncdata.globalatts(ii+1).name = netcdf.inqAttName(ncid,netcdf.getConstant("NC_GLOBAL"),ii);
        ncdata.globalatts(ii+1).value = netcdf.getAtt(ncid,netcdf.getConstant("NC_GLOBAL"),ncdata.globalatts(ii+1).name);
    end

    varids = netcdf.inqVarIDs(ncid); %in the root group only
    for ii=varids
        [varname, ~, ~, numatts] = netcdf.inqVar(ncid,ii);
        vname = matlab.lang.makeValidName(string(varname));
        ncdata.(vname) = netcdf.getVar(ncid,ii);
        aname = vname + "_att";
        ncdata.(aname) = struct();
        for jj=0:(numatts-1)
            ncdata.(aname)(jj+1).name = netcdf.inqAttName(ncid,ii,jj);
            ncdata.(aname)(jj+1).value = netcdf.getAtt(ncid,ii,ncdata.(aname)(jj+1).name);
        end
    end
    netcdf.close(ncid)
end

function data = missing_to_nan(source,attributes)
    data = double(source);
    idx = strcmp('missing_value',{attributes(:).name});
    if any(idx)
        missing_value = attributes(idx).value;
        data(source==missing_value) = NaN;
    end
end
