function lineout = read_NSIDC_OIB_flight_data(fullFileName,epsg)

% ******************************************************************

%  MATLAB FUNCTION read_NSIDC_OIB_flight_data.m

%  Read standard-format NSIDC Operation IceBridge MCoRDS L2
%  ice thickness data files and extract key output variables.

%  Operation IceBridge MCoRDS Level 2 ice thickness files are available
%  in csv format via https://doi.org/10.5067/GDQ0CUCVTE2Q.

%  Inputs
%  fullFileName: string of the full file name and directory
%  epsg: projected co-ordinate system used to calculate x/y co-ordinates
%  [default for Greenland is 3413]

%  Outputs
%  lineout: structure containing variables needed for subsequent analysis

%  This function requires MATLAB version 2020b or later.

% *******************************************************************

%% Set paths to external datasets

% Path to BedMachine geoid correction file
% Set this as the local path to the BedMachine NetCDF file, 
% available via https://doi.org/10.5067/GMEVBWFLWA7X
BM_path = 'BedMachineGreenland-v5.nc';

% Path to rebound correction file
% Set this as the local path to the BedMachine ice unloading correction file, 
% available via https://doi.org/10.18739/A2280509Z
rebound_path = 'Total-Isostatic-Adjustment-Greenland-BedMachinev5.nc';

%% Load data

line = importdata(fullFileName);
data = sortrows(line.data,3,'ascend');   % sort in ascending order of time

%% Separate channels

lat = data(:,1);          % latitude (WGS84)
lon = data(:,2);          % longitude (WGS84)
SoD = data(:,3);          % second of day (UTC)
thk = data(:,4);          % ice thickness (m)
elev = data(:,5);         % aircraft elevation (m WGS84)
frame = data(:,6);        % frame ID (YYYYMMDDSSFFF)
srf_range = data(:,7);    % surface range (m) (aircraft to surface)
bed_range = data(:,8);    % bed range (m) (aircraft to bed)
quality = data(:,9);      % quality flag (1 = good; 0 = bad)

%% Projection into x/y coordinates

p2 = projcrs(epsg);
[x,y] = projfwd(p2,lat,lon);

%% Calculate along-track distance

d = zeros(size(x));

for i = 2:length(x)
    d(i,1) = sqrt((x(i,1)-x(i-1,1))^2 + (y(i,1)-y(i-1,1))^2) + d(i-1,1);
end

%% Calculate along-track flight line azimuth

az = zeros(length(d),1);

for i = 1:length(d)-1
    az(i,1) = rad2deg(atan2((y(i+1,1)-y(i,1)),(x(i+1,1)-x(i,1))));
end

az(end,1) = az(end-1,1);

% Translate onto standard 0–>360 compass bearing
az_bearing = zeros(numel(az),1);
for i = 1:numel(az) 
    if az(i) < 0
        az_bearing(i) = 90 - az(i);
    elseif az(i) >= 0 && az(i) <= 90
        az_bearing(i) = 90 - az(i);
    elseif az(i) > 90 && az(i) <= 180
        az_bearing(i) = 450 - az(i);
    else
        disp('Error: azimuth outside expected range!')
    end
end

%% Calculate elevations relative to global mean sea level (EIGEN-6C4 geoid)

% Extract in situ elevation relative to WGS84 ellipsoid
srf_elev = elev - srf_range;
bed_elev = elev - bed_range;

% Shift onto geoid (global mean sea level)
xg = double(ncread(BM_path,'x'));
yg = double(ncread(BM_path,'y'));
G = ncread(BM_path,'geoid')';
[XG,YG] = meshgrid(xg,yg);
geoid = interp2(XG,YG,double(G),x,y);
elev_geoid = elev - geoid;
srf_elev_geoid = srf_elev - geoid;
bed_elev_geoid = bed_elev - geoid;

%% Adjust bed elevation for the isostatic response to ice unloading

xr = double(ncread(rebound_path,'x'));
yr = double(ncread(rebound_path,'y'));
R = ncread(rebound_path,'solid surface change')';
[XR,YR] = meshgrid(xr,yr);
rebound = interp2(XR,YR,double(R),x,y);
bed_elev_ice_free = bed_elev_geoid + rebound;

%% Create a mask for the modern grounding line

mask = double(ncread(BM_path,'mask')');
gl_mask = mask;
gl_mask(mask==0 | mask==3 | mask==4) = 0;   % outside the grounding line
gl_mask(mask==1 | mask==2) = 1;             % inside the grounding line
gl_mask_interp = interp2(XG,YG,gl_mask,x,y);

%% Remove any duplicate points

[d_unique,idx_unique,~] = unique(d);

%% Assemble output structure

lineout = struct;
lineout.x = x(idx_unique);
lineout.y = y(idx_unique);
lineout.d = d_unique;
lineout.azimuth = az_bearing(idx_unique);
lineout.SoD = SoD(idx_unique);
lineout.frame = frame(idx_unique);
lineout.elev = elev_geoid(idx_unique);
lineout.thk = thk(idx_unique);
lineout.srf_elev = srf_elev_geoid(idx_unique);
lineout.bed_elev = bed_elev_geoid(idx_unique);
lineout.bed_elev_ice_free = bed_elev_ice_free(idx_unique);
lineout.gl_mask = gl_mask_interp(idx_unique);
lineout.quality = quality(idx_unique);

end