function [normals] = dflow_his_unit_vectors(boundary_type,...
    elems_ordered)
%Calculates unit vector for channel-floodplain boundary
%   Detailed explanation goes here

[k,~] = size(elems_ordered);
elem_xy = elems_ordered(:,1:2);
normals = zeros(k-1,6);

% Unit vector calcs

for n=1:k-1;
    
    x = (elem_xy(n,1)+elem_xy(n+1,1))/2;
    y = (elem_xy(n,2)+elem_xy(n+1,2))/2;
    length = sqrt((elem_xy(n+1,1)-elem_xy(n,1))^2 +...
        ((elem_xy(n+1,2)-elem_xy(n,2))^2));
    if n == 1
        down_dist = length / 2;
    else 
        down_dist = normals(n-1,6) + length / 2;
    end
    m = ((elem_xy(n,2)-elem_xy(n+1,2)) /...
        (elem_xy(n,1)-elem_xy(n+1,1)));
    if (elem_xy(n,1) ~= elem_xy(n+1,1)) &&...
            (elem_xy(n,2) ~= elem_xy(n+1,2))
        mprime = -1/m;
        mag = (mprime^2 + 1)^0.5;
        u = -1/mag;
        v = -mprime/mag;
    elseif (elem_xy(n,1) == elem_xy(n+1,1)) && ...
            (elem_xy(n,2) > elem_xy(n+1,2))
        u = -1;
        v = 0;
    elseif (elem_xy(n,1) == elem_xy(n+1,1)) && ...
            (elem_xy(n,2) < elem_xy(n+1,2))
        u = 1;
        v = 0;
    elseif (elem_xy(n,1) > elem_xy(n+1,1)) && ...
            (elem_xy(n,2) == elem_xy(n+1,2))
        u = 0;
        v = 1;
    elseif (elem_xy(n,1) < elem_xy(n+1,1)) && ...
            (elem_xy(n,2) == elem_xy(n+1,2))
        u = 0;
        v = -1;
    end

    % assign 0 for right bank, 1 for left bank (looking downstream)
    % assign 2 for upstream or downstream boundary, 
    % upstream and downstream boundary should be drawn from 
    % right to left bank

    if boundary_type == 0

        if (elem_xy(n,1) < elem_xy(n+1,1)) &&...
                (elem_xy(n,2) > elem_xy(n+1,2))
            u = abs(u) * -1;
            v = abs(v) * -1;
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) &&...
                (elem_xy(n,2) > elem_xy(n+1,2))
            u = abs(u) * -1;
            v = abs(v);
        elseif (elem_xy(n,1) < elem_xy(n+1,1)) &&...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = abs(u);
            v = abs(v) * -1;
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) &&...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = abs(u);
            v = abs(v);
        elseif (elem_xy(n,1) == elem_xy(n+1,1)) && ...
                (elem_xy(n,2) > elem_xy(n+1,2))
            u = -1;
            v = 0;
        elseif (elem_xy(n,1) == elem_xy(n+1,1)) && ...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = 1;
            v = 0;
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) && ...
                (elem_xy(n,2) == elem_xy(n+1,2))
            u = 0;
            v = 1;
        elseif (elem_xy(n,1) < elem_xy(n+1,1)) && ...
                (elem_xy(n,2) == elem_xy(n+1,2))
            u = 0;
            v = -1;
        end

    elseif boundary_type == 1

        if (elem_xy(n,1) < elem_xy(n+1,1)) &&...
                (elem_xy(n,2) > elem_xy(n+1,2))
            u = abs(u);
            v = abs(v);
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) &&...
                (elem_xy(n,2) > elem_xy(n+1,2))
            u = abs(u);
            v = abs(v) * -1;
        elseif (elem_xy(n,1) < elem_xy(n+1,1)) &&...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = abs(u) * -1;
            v = abs(v);
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) &&...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = abs(u) * -1;
            v = abs(v) * -1;
        elseif (elem_xy(n,1) == elem_xy(n+1,1)) && ...
                (elem_xy(n,2) > elem_xy(n+1,2))
            u = 1;
            v = 0;
        elseif (elem_xy(n,1) == elem_xy(n+1,1)) && ...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = -1;
            v = 0;
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) && ...
                (elem_xy(n,2) == elem_xy(n+1,2))
            u = 0;
            v = -1;
        elseif (elem_xy(n,1) < elem_xy(n+1,1)) && ...
                (elem_xy(n,2) == elem_xy(n+1,2))
            u = 0;
            v = 1;
        end

    elseif boundary_type == 2

        if ((elem_xy(n,1) < elem_xy(n+1,1)) &&...
                (elem_xy(n,2) > elem_xy(n+1,2)))
            u = abs(u) * -1;
            v = abs(v) * -1;
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) &&...
                (elem_xy(n,2) > elem_xy(n+1,2))
            u = abs(u) * -1;
            v = abs(v);
        elseif (elem_xy(n,1) < elem_xy(n+1,1)) &&...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = abs(u);
            v = abs(v) * -1;
        elseif (elem_xy(n,1) > elem_xy(n+1,1)) &&...
                (elem_xy(n,2) < elem_xy(n+1,2))
            u = abs(u);
            v = abs(v);
        else
            u = u;
            v = v;
        end    

    end

    normals(n,:)=[x,y,u,v,length,down_dist];

end

end

