function [lb_rb_cell,model_time] =...
    obs2fpboundaries(file_num,filepath,filename,obs_marker)
% Function will extract all information in regard to channel-floodplain
% boundaries

polygon_27km = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_model_boundary_km0_27'],' ');

for i = 1:file_num
    fnameloc(i) = {strcat(filepath, filename(i))};
end

% Import observation file data
for j = 1:file_num
    if j == 1
        all_x_coords = (ncread(char(fnameloc{j}),'station_x_coordinate'))';
        all_y_coords = (ncread(char(fnameloc{j}),'station_y_coordinate'))';
        all_obs_ids_raw = (ncread(char(fnameloc{j}),'station_id'))';
        all_Vx = (ncread(char(fnameloc{j}),'x_velocity'))';
        all_Vy = (ncread(char(fnameloc{j}),'y_velocity'))';
        all_depth = (ncread(char(fnameloc{j}),'Waterdepth'))';
        all_wse = (ncread(char(fnameloc{j}),'waterlevel'))';
        time_stamp = ncread(char(fnameloc{j}),'time');
        max_time = max(time_stamp);
    elseif j > 1
        Vx_j = ncread(char(fnameloc{j}),'x_velocity');
        Vy_j = ncread(char(fnameloc{j}),'y_velocity');
        depth_j = ncread(char(fnameloc{j}),'Waterdepth');
        wse_j = ncread(char(fnameloc{j}),'waterlevel');
        time_stamp_j = ncread(char(fnameloc{j}),'time');
        
        all_Vx = [all_Vx; (Vx_j(:,time_stamp_j>max_time))'];
        all_Vy = [all_Vy; (Vy_j(:,time_stamp_j>max_time))'];
        all_depth = [all_depth; (depth_j(:,time_stamp_j>max_time))'];
        all_wse = [all_wse; (wse_j(:,time_stamp_j>max_time))'];
        time_stamp =...
            [time_stamp; time_stamp_j(time_stamp_j>max_time,:)];
        max_time = max(time_stamp);
        clear Vx_j Vy_j depth_j wse_j time_stamp_j
    end
end

%% Reduce all boundary points to those within 27 km polygon

in_polygon_yn = inpolygon(all_x_coords,all_y_coords,...
    polygon_27km(:,1),polygon_27km(:,2));

all_x_coords = all_x_coords(:,in_polygon_yn(:)==1);
all_y_coords = all_y_coords(:,in_polygon_yn(:)==1);
all_depth = all_depth(:,in_polygon_yn(:)==1);
all_obs_ids_raw = all_obs_ids_raw((in_polygon_yn(:)==1),:);
all_Vx = all_Vx(:,in_polygon_yn(:)==1);
all_Vy = all_Vy(:,in_polygon_yn(:)==1);
all_wse = all_wse(:,in_polygon_yn(:)==1);

%% Split all points into left and right boundaries

[~,obs_pts] = size(all_x_coords);
left_elems = sum(all_obs_ids_raw(:,1) == 'l');
right_elems = sum(all_obs_ids_raw(:,1) == 'r');
[time_steps,~] = size(time_stamp);

l_i = 0;
r_i = 0;
left_elem_pt = zeros(left_elems,2);
left_elem_Vx = zeros(left_elems, time_steps);
left_elem_Vy = zeros(left_elems, time_steps);
left_elem_depth = zeros(left_elems, time_steps);

right_elem_pt = zeros(right_elems,2);
right_elem_Vx = zeros(right_elems, time_steps);
right_elem_Vy = zeros(right_elems, time_steps);
right_elem_depth = zeros(right_elems, time_steps);

for i = 1:obs_pts
    if all_obs_ids_raw(i,1) == obs_marker(1) % obs_marker(1) == 'l'
         l_i = l_i + 1;
         left_elem_pt(l_i,1) = all_x_coords(1,i);
         left_elem_pt(l_i,2) = all_y_coords(1,i);
         left_elem_Vx(l_i,:) = all_Vx(:,i)';
         left_elem_Vy(l_i,:) = all_Vy(:,i)';
         left_elem_depth(l_i,:) = all_depth(:,i)';
    elseif all_obs_ids_raw(i,1) == obs_marker(2) % obs_marker(2) == 'r'
         r_i = r_i + 1;
         right_elem_pt(r_i,1) = all_x_coords(1,i);
         right_elem_pt(r_i,2) = all_y_coords(1,i);
         right_elem_Vx(r_i,:) = all_Vx(:,i)';
         right_elem_Vy(r_i,:) = all_Vy(:,i)';
         right_elem_depth(r_i,:) = all_depth(:,i)';
    end
    
end

left_mid_point_xy = zeros(left_elems-1,2);
right_mid_point_xy = zeros(right_elems-1,2);

for p_i = 1:left_elems-1
    left_mid_point_xy(p_i,1) = (left_elem_pt(p_i,1) +...
        left_elem_pt(p_i+1,1)) / 2;
    left_mid_point_xy(p_i,2) = (left_elem_pt(p_i,2) +...
        left_elem_pt(p_i+1,2)) / 2;
end

for p_j = 1:right_elems-1
    right_mid_point_xy(p_j,1) = (right_elem_pt(p_j,1) +...
        right_elem_pt(p_j+1,1)) / 2;
    right_mid_point_xy(p_j,2) = (right_elem_pt(p_j,2) +...
        right_elem_pt(p_j+1,2)) / 2;
end


%% Create matrices for boundary lengths based on element points

% left_segs_pts
% x_coord | y_coord | u | v | length | downstream_distance

% Left bank boundary calcs, boundary type = 1
left_seg_pts = dflow_his_unit_vectors(1,left_elem_pt);
[left_Q, left_Mo_flux, left_Mo_trans, left_Mo_length, l_area, l_length]=...
    dflow_mass_momentum_his_calc(1,left_elem_Vx, left_elem_Vy,...
    left_elem_depth,time_steps,left_seg_pts);

% Right bank boundary calcs, boundary type = 0
right_seg_pts = dflow_his_unit_vectors(0,right_elem_pt);
[right_Q, right_Mo_flux, right_Mo_trans, right_Mo_length, r_area,...
    r_length] =...
    dflow_mass_momentum_his_calc(0,right_elem_Vx, right_elem_Vy,...
    right_elem_depth,time_steps,right_seg_pts);

%% Calculate distances upstream and downstream of each segment

left_distances = updownstream_distances(left_seg_pts,left_elems);
right_distances = updownstream_distances(right_seg_pts,right_elems);

%% Assign left and right information to cells
% Cell 1 - element coordinates (NAD83 m)
% Cell 2 - Vx (m/s)
% Cell 3 - Vy (m/s)
% Cell 4 - Depth (m)
% Cell 5 - midpoint of flux line coordinate
% Cell 6 - Discharge (m^3/s)
% Cell 7 - Momentum flux per area (N/m^2)
% Cell 8 - Momentum flux (N)
% Cell 9 - Momentum per unit length (N/m)
% Cell 10 - Upstream and downstream distances (m)
% Cell 11 - Cumulative downstream distance (m)
% Cell 12 - Area of flux (m^2)
% Cell 13 - Interface element length (m)

lb_cell = cell(12,1);
rb_cell = cell(12,1);

lb_cell(1) = {left_elem_pt};
lb_cell(2) = {left_elem_Vx};
lb_cell(3) = {left_elem_Vy};
lb_cell(4) = {left_elem_depth};
lb_cell(5) = {left_mid_point_xy};
lb_cell(6) = {left_Q};
lb_cell(7) = {left_Mo_flux};
lb_cell(8) = {left_Mo_trans};
lb_cell(9) = {left_Mo_length};
lb_cell(10) = {left_distances};
lb_cell(11) = {left_distances(:,1)};
lb_cell(12) = {l_area};
lb_cell(13) = {l_length};

rb_cell(1) = {right_elem_pt};
rb_cell(2) = {right_elem_Vx};
rb_cell(3) = {right_elem_Vy};
rb_cell(4) = {right_elem_depth};
rb_cell(5) = {right_mid_point_xy};
rb_cell(6) = {right_Q};
rb_cell(7) = {right_Mo_flux};
rb_cell(8) = {right_Mo_trans};
rb_cell(9) = {right_Mo_length};
rb_cell(10) = {right_distances};
rb_cell(11) = {right_distances(:,1)};
rb_cell(12) = {r_area};
rb_cell(13) = {r_length};

lb_rb_cell = cell(1,2);
lb_rb_cell(1,1) = {lb_cell};
lb_rb_cell(1,2) = {rb_cell};

model_time = time_stamp;

end