function Astruct = tributarymatrix( hdf_input, y_frame, attach_pts, varargin )

% [ A,  ind_ts ] = tributarymatrix( hdf_filename , y_frame, attach_pts, varargin ) 
% [ A,  ind_ts ] = tributarymatrix( hdf_struct , y_frame, attach_pts, varargin ) 
% [ A,  ind_ts ] = tributarymatrix( [] , y_frame, attach_pts, varargin ) --> prompts user to select HDF file

plot_on = 0; % (0 ==> no plots, 1 ==> some plots, 2 ==> all plots) default can be overridden by optional user input

n_frames = size(y_frame,1);

if ischar( hdf_input )
    hdf_filename = hdf_input;
    hdf_struct = load_hdf_coords( hdf_filename );
elseif isstruct( hdf_input )
    hdf_struct = hdf_input;
else
    [filename pathname]= uigetfile({'*.hdf','HDF-files (*.hdf)'},'Select file for plotting taps');
    hdf_full_filename = [pathname filename]
    hdf_struct = load_hdf_coords( hdf_full_filename );
    assignin('base','hdf_full_filename', hdf_full_filename); % make filename available in base workspace for convenience
end

% Extract model dimensions (full-scale) from HDF data structure:
L = hdf_struct.Building_Length_Ft;
W = hdf_struct.Building_Width_Ft;
H = hdf_struct.Building_Height_Ft; % eave height
R = hdf_struct.Roof_Slope_in_12*W/2; % slope is given as rise (in inches) per foot of horizontal distance

% Initialize dimensions of structure with dimensions of model (default):
L0 = L; W0 = W; H0 = H; R0 = R;

% Process any optional inputs: assign alternative dimensions of structure, activate plotting
if nargin>3
    if mod(nargin-3,2)~=0
        error('Optional input arguments must occur in pairs: ...''PropertyName'',PropertyValue... ')
    else
        for i = 1:2:nargin-4
            if ischar(varargin{i}) & max(strcmp(varargin{i},{'plot_on','L0','W0','H0','R0'}))
                if isnumeric(varargin{i+1})
                    eval([varargin{i} '=' num2str(varargin{i+1}) ';']);
                else
                    error(['Non-numeric value specified for optional input >> ' varargin{i} ' <<']);
                end
            else
                error(['The specified optional input >> ' num2str(varargin{i}) ' << is an invalid variable name.' ]);
            end
        end
    end
end

if plot_on, plot3Dtaps( hdf_struct ); h1 = gcf; pause; end;

tap_c0_3d = hdf_struct.Tap_Coordinates_3D; % tap #, face #, x-coord, y-coord, z-coord

% Face numbering and coordinate system in UWO HDF files:
%   Face 1 is the leeward face for zero wind direction, Face 3 is windward
%   (z-axis is upwards, walls are folded out here for display)
%   Internal pressure taps are labeled Face 0
%                                          
%                                          |x-axis
%                                          |
%                __________________________|
%               |                          |
%               |            4             |
%          _____|__________________________|_____
%         |     |                          |     |
%         |     |            6             |     |
%  0 deg  |     |                          |     |
%   ----> |  3  |----------ridge-----------|  1  |
%         |     |                          |     |
%         |     |            5             |     |
%  _______|_____|__________________________|_____|
%  y-axis       |                          |Origin
%               |            2             |
%               |__________________________|
%
%   Shift origin to opposite corner (on windward face for zero wind dir),
%   eliminate internal taps and end walls, and renumber faces as follows: 
%                __________________________
%               |                          |
%               |            1             |
%               |__________________________|_________________
%               |Origin                    |          y-axis
%               |            2             |     
%     0 deg     |                          |     
%       ---->   |--- ridge: face = 0 ------|    
%               |                          |     
%               |            3             |     
%               |__________________________|
%               |                          |
%               |            4             |
%               |__________________________|
%               |
%               |
%               |x-axis

% shift origin to opposite corner of model:
tap_c_3d = tap_c0_3d;
tap_c_3d(:,3) = W-tap_c0_3d(:,3);
tap_c_3d(:,4) = L-tap_c0_3d(:,4);
% eliminate internal pressure taps:
ext_ind = find(tap_c_3d(:,2)~=0);
tap_c_3d = tap_c_3d(ext_ind,:);
% eliminate taps on end walls:
ind_no_ends = intersect(find(tap_c_3d(:,4)~=0), find(tap_c_3d(:,4)~=L));
tap_c_3d = tap_c_3d(ind_no_ends,:);
% assign new face numbers:
tap_c_3d(find(tap_c_3d(:,3)==0),2) = 1;
tap_c_3d(find(tap_c_3d(:,3)>0 & tap_c_3d(:,3)<W/2 ),2) = 2;
tap_c_3d(find(tap_c_3d(:,3)==W/2 ),2) = 0;
tap_c_3d(find(tap_c_3d(:,3)>W/2 & tap_c_3d(:,3)<W),2) = 3;
tap_c_3d(find(tap_c_3d(:,3)==W),2) = 4;

if plot_on
    color_list = {'g','k','b','r','c'};
    figure(h1);
    for i = 0:4
        ind_face = find(tap_c_3d(:,2)==i);
        if length(ind_face)
            plot3(tap_c_3d(ind_face,3),tap_c_3d(ind_face,4),tap_c_3d(ind_face,5),[color_list{i+1} '.'])
            hold on;
            x_face = mean(tap_c_3d(ind_face,3));
            y_face = mean(tap_c_3d(ind_face,4));
            z_face = mean(tap_c_3d(ind_face,5));
            text(x_face, y_face, z_face, ['Face ' num2str(i)],'Color',color_list{i+1});
        end
    end
    xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis');
    title('Tap coordinates with origin shifted to opposite corner, end walls removed, and faces renumbered');
    axis equal
    axis vis3d
    hold off;
    pause;
end

% Scale tap coordinates to match dimensions of structure of interest:
scale = 0;
if H0~=H | W0~=W | L0~=L | R0~=R
    scale =1;
    % scale x-coordinates:
    tap_c_3d(:,3) = tap_c_3d(:,3)*W0/W;
    % scale y-coordinates:
    tap_c_3d(:,4) = tap_c_3d(:,4)*L0/L;
    % scale z-coordinates:
    roof_ind = find(tap_c_3d(:,5)>H);
    wall_ind = find(tap_c_3d(:,5)<=H);
    tap_c_3d(roof_ind,5) = H0+(tap_c_3d(roof_ind,5)-H)*R0/R;
    tap_c_3d(wall_ind,5) = tap_c_3d(wall_ind,5)*H0/H;
end

if plot_on && scale
    color_list = {'g','k','b','r','c'};
    figure(h1);
    for i = 0:4
        ind_face = find(tap_c_3d(:,2)==i);
        if length(ind_face)
            plot3(tap_c_3d(ind_face,3),tap_c_3d(ind_face,4),tap_c_3d(ind_face,5),[color_list{i+1} '.'])
            hold on;
            x_face = mean(tap_c_3d(ind_face,3));
            y_face = mean(tap_c_3d(ind_face,4));
            z_face = mean(tap_c_3d(ind_face,5));
            text(x_face, y_face, z_face, ['Face ' num2str(i)],'Color',color_list{i+1});
        end
    end
    xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis');
    title('Tap coordinates scaled to dimensions of structure of interest');
    axis equal
    axis vis3d
    hold off;
    pause;
end
    
% Define 2D local coordinates for each face:
% (y-axis coincides with global y-axis, w-axis explained below)
%               _____________________________________________
%               |Origin                    |          y-axis
%     0 deg     |                          |     
%       ---->   |                          |    
%               |__________________________|
%               |
%               |
%               |w-axis
%               
% On the side walls (faces 1 & 4) the local w-axis is in the vertical plane;
% On the roof (faces 2 and 3), the local w-axis is in the horizontal plane.
% The w-coordinate increases upwards on face 1 and downwards on face 4;
% The w-coordinate increases from left to right on faces 2 and 3.
%                                        .
%                                   .         .
%                              .    2           3  .
%                         .---------->   .------->      .
%                                                       |
%                         .                             |
%                       1 ^                             v 4
%                         |                             .
%                         |                             .
%                         .                             .
%          

tap_lc = zeros(size(tap_c_3d,1),4);
tap_lc(:,1:2) = tap_c_3d(:,1:2); % tap #, face #
% y-coordinates:
tap_lc(:,4) = tap_c_3d(:,4); 
% s-coordinates:
tap_lc(find(tap_lc(:,2)==1),3)= tap_c_3d(find(tap_lc(:,2)==1),5);
tap_lc(find(tap_lc(:,2)==4),3)= tap_c_3d(find(tap_lc(:,2)==4),5);
hyp0 = sqrt((W0/2)^2+R0^2); % hypotenuse of triangle formed by W/2 and R
tap_lc(find(tap_lc(:,2)==2),3)= tap_c_3d(find(tap_lc(:,2)==2),3)*2*hyp0/W0;
tap_lc(find(tap_lc(:,2)==3),3)= (W0-tap_c_3d(find(tap_lc(:,2)==3),3))*2*hyp0/W0;

% duplicate the taps along the ridge line to be included in both faces 2 and 3:
tap_lc(find(tap_lc(:,2)==0),3)= W0/2;
tap_ridge_dup = tap_lc(find(tap_lc(:,2)==0),:);
tap_ridge_dup(:,2) = 3;
tap_lc(find(tap_lc(:,2)==0),2)= 2;
tap_lc = [tap_lc; tap_ridge_dup];

% define bounds of each face in s-y plane:
lc_bnds = [0 H0 0 L0; 0 hyp0 0 L0; 0 hyp0 0 L0; 0 H0 0 L0];

% Loop through faces, evaluate bounds of tributary area for each tap:
tap_trib_bnds = zeros(size(tap_lc,1),6); % initialize matrix for storing results: 
                                         %  tap #, face #, s_min, s_max, y_min, y_max

if plot_on, h2 = figure; end;
for i = 1:4
    ind_i = find(tap_lc(:,2)==i);
    n_tap_i = length(ind_i);
    tap_lc_i = tap_lc(ind_i,:);
    lc_bnds_i = lc_bnds(i,:);

    tap_trib_bnds_i = taptribbnds( tap_lc_i, lc_bnds_i );
    
    if plot_on
        figure(h2);
        plot(tap_lc_i(:,3), tap_lc_i(:,4), '+');
        hold on;
        plot([lc_bnds_i(1:2) lc_bnds_i(2:-1:1) lc_bnds_i(1)],...
            [lc_bnds_i(3) lc_bnds_i(3) lc_bnds_i(4) lc_bnds_i(4) lc_bnds_i(3)],'-');
        hold off;
        title(['Pressure tap layout on face ' num2str(i)]);
        pause;
        for j = 1:n_tap_i
            figure(h2);
            hold on;
            plot(tap_lc_i(j,3), tap_lc_i(j,4),'o');
            plot([tap_trib_bnds_i(j,3:4) tap_trib_bnds_i(j,4:-1:3) tap_trib_bnds_i(j,3)],...
                [tap_trib_bnds_i(j,5) tap_trib_bnds_i(j,5) tap_trib_bnds_i(j,6) tap_trib_bnds_i(j,6) tap_trib_bnds_i(j,5)], '-');
            hold off;
        end
        title(['Pressure tap layout on face ' num2str(i) ' with tributary areas indicated']);
        pause;
    end
    tap_trib_bnds(ind_i,:) = tap_trib_bnds_i;
end

% Reflect local tap coordinates and tributary bounds about axes of symmetry:
% (b) Reflect about x = W0/2 (ridge axis) by reflecting face numbers:
tap_lc_b = tap_lc;
tap_lc_b(find(tap_trib_bnds(:,2)==1),2)=4;
tap_lc_b(find(tap_trib_bnds(:,2)==2),2)=3;
tap_lc_b(find(tap_trib_bnds(:,2)==3),2)=2;
tap_lc_b(find(tap_trib_bnds(:,2)==4),2)=1;
tap_trib_bnds_b = tap_trib_bnds; 
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==1),2)=4;
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==2),2)=3;
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==3),2)=2;
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==4),2)=1;
% (c) Reflect about y = L0/2 by reflecting y-coordinates:
tap_lc_c = tap_lc;
tap_lc_c(:,4) = L0-tap_lc(:,4);
tap_trib_bnds_c = tap_trib_bnds;
tap_trib_bnds_c(:,5)=L0-tap_trib_bnds(:,6); % y_min becomes y_max after reflection
tap_trib_bnds_c(:,6)=L0-tap_trib_bnds(:,5); % y_max becomes y_min after reflection
% (d) Reflect about both axis by reflecting face numbers and y-coordinates:
tap_lc_d = tap_lc;
tap_lc_d(find(tap_trib_bnds(:,2)==1),2)=4;
tap_lc_d(find(tap_trib_bnds(:,2)==2),2)=3;
tap_lc_d(find(tap_trib_bnds(:,2)==3),2)=2;
tap_lc_d(find(tap_trib_bnds(:,2)==4),2)=1;
tap_lc_d(:,4) = L0-tap_lc(:,4);
tap_trib_bnds_d = tap_trib_bnds;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==1),2)=4;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==2),2)=3;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==3),2)=2;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==4),2)=1;
tap_trib_bnds_d(:,5)=L0-tap_trib_bnds(:,6); % y_min becomes y_max after reflection
tap_trib_bnds_d(:,6)=L0-tap_trib_bnds(:,5); % y_max becomes y_min after reflection

% Select the taps whose tributary area overlaps the tributary area of the specified frames:
ind_trib = [];
for f = 1:n_frames
    ind_trib_ab = intersect(find(tap_trib_bnds(:,6)>y_frame(f,1)),find(tap_trib_bnds(:,5)<y_frame(f,3)));
    ind_trib_cd = intersect(find(tap_trib_bnds_c(:,6)>y_frame(f,1)),find(tap_trib_bnds_c(:,5)<y_frame(f,3)));
    ind_trib = union( ind_trib_ab, ind_trib_cd );
%end

    tap_lc = tap_lc(ind_trib,:);
    tap_lc_b = tap_lc_b(ind_trib,:);
    tap_lc_c = tap_lc_c(ind_trib,:);
    tap_lc_d = tap_lc_d(ind_trib,:);
    tap_trib_bnds = tap_trib_bnds(ind_trib,:);
    tap_trib_bnds_b = tap_trib_bnds_b(ind_trib,:);
    tap_trib_bnds_c = tap_trib_bnds_c(ind_trib,:);
    tap_trib_bnds_d = tap_trib_bnds_d(ind_trib,:);
    tap_num = tap_lc(:,1);
    [tf, tap_ind] = ismember(tap_num, hdf_struct.Tap_Position_List);
    Astruct(f).tap_ind = tap_ind;

    % error checking (can be removed eventually)
    % min(tap_lc(:,1)==tap_trib_bnds(:,1));
    if ismember(0, tf)
        error('Tap number not found in "Tap_Position_List" field in HDF file.');
    elseif ~isempty(setxor( tap_num, hdf_struct.Tap_Position_List(tap_ind)))
        error('Error in matching tap number with "Tap_Position_List" in HDF file.');
    end

    if plot_on
        % Find indices in tap_c_3d of taps that are tributary to current frame:
        [ tf2, ind_3d ] = ismember( tap_num, tap_c_3d(:,1) );
        color_list = {'b','r','g','c','m','y'};

        figure(h1);
        plot3([0 W W W/2 0 0],[0 0 0 0 0 0],[0 0 H H+R H 0],'k-', ...
        [0 W W W/2 0 0],[L L L L L L],[0 0 H H+R H 0],'k-', ...
        [0 0],[0 L],[0 0],'k-',[0 0],[0 L],[H H],'k-',...
        [W W],[0 L],[0 0],'k-',[W W],[0 L],[H H],'k-',...
        [W/2 W/2],[0 L],[H+R H+R],'k-'); hold on;
        plot3(tap_c_3d(:,3),tap_c_3d(:,4),tap_c_3d(:,5),'k.');
        plot3(tap_c_3d(ind_3d,3),tap_c_3d(ind_3d,4),tap_c_3d(ind_3d,5),'ro');
%    for i = 1:n_frames
        f_clr = mod(f,6); f_clr(find(i_clr==0))=6;
        clr = color_list{f_clr};
        plot3([0 0 W0/2 W0 W0], y_frame(f,1)*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
        plot3([0 0 W0/2 W0 W0], y_frame(f,2)*ones(1,5),[0 H0 H0+R0 H0 0],[clr '-'],'LineWidth',2);
        plot3([0 0 W0/2 W0 W0], y_frame(f,3)*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
        plot3([0 0 W0/2 W0 W0], (L0-y_frame(f,1))*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
        plot3([0 0 W0/2 W0 W0], (L0-y_frame(f,2))*ones(1,5),[0 H0 H0+R0 H0 0],[clr '-'],'LineWidth',2);
        plot3([0 0 W0/2 W0 W0], (L0-y_frame(f,3))*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
%    end
        xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis');
        title('Tap coordinates plotted with specified frame and tributary taps circled in red.');
        axis equal
        axis vis3d
        rotate3d on;
        hold off;
        pause;
    end

%for f = 1:n_frames
    Astruct(f).Aa = trib_mx( y_frame(f,:), tap_trib_bnds, tap_lc );
    Astruct(f).Ab = trib_mx( y_frame(f,:), tap_trib_bnds_b, tap_lc_b );
    Astruct(f).Ac = trib_mx( y_frame(f,:), tap_trib_bnds_c, tap_lc_c );
    Astruct(f).Ad = trib_mx( y_frame(f,:), tap_trib_bnds_d, tap_lc_d );
end

function A_mx = trib_mx( y_f, tap_bnds, lc )

    % find indices of taps that overlap current frame:
    trib_ind = intersect(find(tap_bnds(:,6)>y_f(1)),find(tap_bnds(:,5)<y_f(3)));
    % Adjust tributary bounds to lie within the tributary area of the current frame:
    tap_bnds(find(tap_bnds(:,6)>y_f(3)),6) = y_f(3);
    tap_bnds(find(tap_bnds(:,5)<y_f(1)),5) = y_f(1);

    A_mx = zeros(size(attach_pts,1), length(tap_ind)); % initialize tributary area matrix
    if plot_on==2, h3 = figure; end;
    for i = 1:4
        % find taps on current face that overlap current frame:
        tap_ind_i = intersect( trib_ind, find(lc(:,2)==i));
        n_tap_i = length(tap_ind_i);
        lc_i = lc(tap_ind_i,:);
        lc_bnds_i = lc_bnds(i,:);
        tap_bnds_i = tap_bnds(tap_ind_i, :);

        if plot_on==2
            figure(h3);
            plot(lc_i(:,3), lc_i(:,4), '+');
            hold on;
            plot(lc_bnds_i(1:2),[y_f(1) y_f(1)],'k:',...
                lc_bnds_i(1:2),[y_f(2) y_f(2)],'k-',...
                lc_bnds_i(1:2),[y_f(3) y_f(3)],'k:');
            hold off;
            title(['Pressure taps on face ' num2str(i) ' that are tributary to specified frame.']);
            pause;
            for j = 1:n_tap_i
                figure(h3);
                hold on;
                plot(lc_i(j,3), lc_i(j,4),'o');
                plot([tap_bnds_i(j,3:4) tap_bnds_i(j,4:-1:3) tap_bnds_i(j,3)],...
                    [tap_bnds_i(j,5) tap_bnds_i(j,5) tap_bnds_i(j,6) tap_bnds_i(j,6) tap_bnds_i(j,5)], '-');
                hold off;
            end
            title(['Pressure tap layout on face ' num2str(i) ' with tributary areas indicated']);
            pause;
        end
        % find attachment points on current face:
        attach_ind_i = find(attach_pts(:,1)==i);
        attach_s_i = attach_pts(attach_ind_i,2); % s-coordinates

        for j = 0:length(attach_ind_i)

            if j == 0
                s_min = lc_bnds_i(1);
                s_max = attach_s_i(j+1);
            elseif j == length(attach_ind_i)
                s_min = attach_s_i(j);
                s_max = lc_bnds_i(2);
            else
                s_min = attach_s_i(j);
                s_max = attach_s_i(j+1);
            end
            dS = s_max - s_min;

            if dS<=0
                continue
            end
            % adjust tributary bounds to lie within current panel:
            tap_bnds_ij = tap_bnds_i;
            tap_bnds_ij(find(tap_bnds_ij(:,3)>s_max),3) = s_max;
            tap_bnds_ij(find(tap_bnds_ij(:,4)>s_max),4) = s_max;
            tap_bnds_ij(find(tap_bnds_ij(:,3)<s_min),3) = s_min;
            tap_bnds_ij(find(tap_bnds_ij(:,4)<s_min),4) = s_min;

            if y_f(1)~=y_f(2)
                % treat contributions from span A preceeding frame of interest:
                dY = y_f(2)-y_f(1);
                tap_bnds_ij_A = tap_bnds_ij;
                tap_bnds_ij_A(find(tap_bnds_ij_A(:,5)>y_f(2)),5) = y_f(2);
                tap_bnds_ij_A(find(tap_bnds_ij_A(:,6)>y_f(2)),6) = y_f(2);
                s_p = tap_bnds_ij_A(:,3)'-s_min; % "s prime" (lever arm from panel origin to corner of tributary area)
                ds = tap_bnds_ij_A(:,4)'-tap_bnds_ij_A(:,3)'; 
                y_p = y_f(2)-tap_bnds_ij_A(:,6)'; % "y prime" (lever arm from panel origin to corner of tributary area)
                dy = tap_bnds_ij_A(:,6)'-tap_bnds_ij_A(:,5)';

                if j == 0
                    % first purlin: all load transferred to attachment point at s_max:
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a2 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a1 = zeros(size(a2));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a2_1 = ds.*dy.*(y_p+.5*dy)/dY;
                    a1_1 = zeros(size(a2_1));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                elseif j == length(attach_ind_i)
                    % last purlin: all load transferred to attachment point at s_min:
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a1 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a2 = zeros(size(a1));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_1 = ds.*dy.*(y_p+.5*dy)/dY;
                    a2_1 = zeros(size(a1_1));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                else
                    % tributary area vectors for 2 attachment points at w_min (1) and w_max (2):
                    a1 = ds.*dy.*(dS-s_p-.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    a2 = ds.*dy.*(s_p+.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_1 = ds.*dy.*(dS-s_p-.5*ds).*(y_p+.5*dy)/(dS*dY);
                    a2_1 = ds.*dy.*(s_p+.5*ds).*(y_p+.5*dy)/(dS*dY);
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                end

                if plot_on==2
                    ind_nz = find(a1 | a2);
                    for k = ind_nz
                        figure(h3);
                        clf;
                        plot(lc_i(:,3), lc_i(:,4), '+'); % plot all taps on current face
                        hold on;
                        plot(lc_bnds_i(1:2),[y_f(1) y_f(1)],'k:',...
                            lc_bnds_i(1:2),[y_f(2) y_f(2)],'k-',...
                            lc_bnds_i(1:2),[y_f(3) y_f(3)],'k:'); % plot frames
                        plot([s_min s_min], [y_f(1) y_f(2)],'k-',...
                            [s_max s_max], [y_f(1) y_f(2)],'k-'); % plot girts/purlins
                        plot([s_min s_max],[y_f(2) y_f(2)],'kx'); % mark attachment points to current frame
                        plot([s_min s_max],[y_f(1) y_f(1)],'kx'); % mark attachment points to previous frame
                        plot(lc_i(k,3), lc_i(k,4),'o'); % circle current tap
                        % plot bounds of tributary area for current tap on current span:
                        plot(s_min+s_p(k)+[0 ds(k) ds(k) 0 0],y_f(2)-y_p(k)-[0 0 dy(k) dy(k) 0],'-');
                        % label total tributary area of current tap:
                        text(lc_i(k,3), lc_i(k,4),...
                            ['Area = ' num2str(ds(k)*dy(k)) ' (sum = ' num2str(a1(k)+a2(k)+a1_1(k)+a2_1(k)) ')']); 
                        % label portion of tributary area applied to each attachment point:
                        text(s_min,y_f(2),['Area = ' num2str(a1(k))]);
                        text(s_max,y_f(2),['Area = ' num2str(a2(k))]);
                        text(s_min,y_f(1),['Area = ' num2str(a1_1(k))]);
                        text(s_max,y_f(1),['Area = ' num2str(a2_1(k))]);
                        title(['Application of tributary area for tap # ' num2str(lc_i(k,1)) ' on face ' num2str(i)]);
                        hold off;
                        pause;
                    end
                end

            end
            if y_f(3)~=y_f(2)
                % treat contributions from span B, beyond frame of interest:
                dY = y_f(3)-y_f(2);
                tap_bnds_ij_B = tap_bnds_ij;
                tap_bnds_ij_B(find(tap_bnds_ij_B(:,5)<y_f(2)),5) = y_f(2);
                tap_bnds_ij_B(find(tap_bnds_ij_B(:,6)<y_f(2)),6) = y_f(2);
                s_p = tap_bnds_ij_B(:,3)'-s_min;
                ds = tap_bnds_ij_B(:,4)'-tap_bnds_ij_B(:,3)';
                y_p = tap_bnds_ij_B(:,5)'-y_f(2);
                dy = tap_bnds_ij_B(:,6)'-tap_bnds_ij_B(:,5)';

                if j == 0
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a2 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a1 = zeros(size(a2));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a2_3 = ds.*dy.*(y_p+.5*dy)/dY;
                    a1_3 = zeros(size(a2_3));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                elseif j == length(attach_ind_i)
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a1 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a2 = zeros(size(a1));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_3 = ds.*dy.*(y_p+.5*dy)/dY;
                    a2_3 = zeros(size(a1_3));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                else
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a1 = ds.*dy.*(dS-s_p-.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    a2 = ds.*dy.*(s_p+.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_3 = ds.*dy.*(dS-s_p-.5*ds).*(y_p+.5*dy)/(dS*dY);
                    a2_3 = ds.*dy.*(s_p+.5*ds).*(y_p+.5*dy)/(dS*dY);
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                end

                if plot_on==2
                    ind_nz = find(a1 | a2);
                    for k = ind_nz
                        figure(h3);
                        clf;
                        plot(lc_i(:,3), lc_i(:,4), '+'); % plot all taps on current face
                        hold on;
                        plot(lc_bnds_i(1:2),[y_f(1) y_f(1)],'k:',...
                            lc_bnds_i(1:2),[y_f(2) y_f(2)],'k-',...
                            lc_bnds_i(1:2),[y_f(3) y_f(3)],'k:'); % plot frames
                        plot([s_min s_min], [y_f(2) y_f(3)],'k-',...
                            [s_max s_max], [y_f(2) y_f(3)],'k-'); % plot girts/purlins
                        plot([s_min s_max],[y_f(2) y_f(2)],'kx'); % mark attachment points to current frame
                        plot([s_min s_max],[y_f(3) y_f(3)],'kx'); % mark attachment points to next frame
                        plot(lc_i(k,3), lc_i(k,4),'o'); % circle current tap
                        % plot bounds of tributary area for current tap on current span:
                        plot(s_min+s_p(k)+[0 ds(k) ds(k) 0 0],y_f(2)+y_p(k)+[0 0 dy(k) dy(k) 0],'-');
                        % label total tributary area of current tap:
                        text(lc_i(k,3), lc_i(k,4),...
                            ['Area = ' num2str(ds(k)*dy(k)) ' (sum = ' num2str(a1(k)+a2(k)+a1_3(k)+a2_3(k)) ')']); 
                        % label portion of tributary area applied to each attachment point:
                        text(s_min,y_f(2),['Area = ' num2str(a1(k))]);
                        text(s_max,y_f(2),['Area = ' num2str(a2(k))]);
                        text(s_min,y_f(3),['Area = ' num2str(a1_3(k))]);
                        text(s_max,y_f(3),['Area = ' num2str(a2_3(k))]);
                        title(['Application of tributary area for tap # ' num2str(lc_i(k,1)) ' on face ' num2str(i)]);
                        hold off;
                        pause;
                    end
                end

            end

        end

    end

end

function tap_trib_bnds_i = taptribbnds( tap_lc_i, lc_bnds_i )

    tap_row_y = unique(tap_lc_i(:,4));
    tap_col_w = unique(tap_lc_i(:,3));
    n_tap_row_y = zeros(size(tap_row_y));
    n_tap_col_w = zeros(size(tap_col_w));
    for I = 1:length(tap_row_y)
        n_tap_row_y(I) = length(find(tap_lc_i(:,4) == tap_row_y(I))); % number of taps per row
    end
    for I = 1:length(tap_col_w)
        n_tap_col_w(I) = length(find(tap_lc_i(:,3) == tap_col_w(I))); % number of taps per column
    end
    d_n_row = diff(n_tap_row_y); % change in number of taps per row
    d_n_col = diff(n_tap_col_w); % change in number of taps per column
    
    % if tap array is a simple grid: constant number of taps per row and per column
    if max(abs(d_n_row))==0 & max(abs(d_n_col))==0 

        tap_trib_bnds_i = simplegrid( tap_lc_i, lc_bnds_i );
    
    % if the number of taps per row changes only once on the face:
    elseif length(find(d_n_row))==1
        % slice the current face horizontally to yeild two simple grids:
        y_slice = mean([tap_row_y(find(d_n_row)) tap_row_y(find(d_n_row)+1)]);
        ind_A = find(tap_lc_i(:,4)<y_slice);
        ind_B = find(tap_lc_i(:,4)>y_slice);
        tap_lc_i_A = tap_lc_i(ind_A,:);
        tap_lc_i_B = tap_lc_i(ind_B,:);
        lc_bnds_i_A = lc_bnds_i; lc_bnds_i_A(4) = y_slice;
        lc_bnds_i_B = lc_bnds_i; lc_bnds_i_B(3) = y_slice;

        tap_trib_bnds_i_A = simplegrid( tap_lc_i_A, lc_bnds_i_A );
        tap_trib_bnds_i_B = simplegrid( tap_lc_i_B, lc_bnds_i_B );
        
        tap_trib_bnds_i = zeros(size(tap_lc_i,1),6);
        tap_trib_bnds_i(ind_A, :) = tap_trib_bnds_i_A;
        tap_trib_bnds_i(ind_B, :) = tap_trib_bnds_i_B;

    else
        error('Script can currently handle only one change in the number of taps per row on any face.')
    end

end

function tap_trib_bnds_i = simplegrid( tap_lc_i, lc_bnds_i )

    tap_trib_bnds_i = zeros(size(tap_lc_i,1),6);
    tap_trib_bnds_i(:,1:2) = tap_lc_i(:,1:2);
    tap_row_y = unique(tap_lc_i(:,4));
    tap_col_w = unique(tap_lc_i(:,3));

    % assign the left edge of the current face as the lower bound for taps in the first column:
    tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(1)),3) = lc_bnds_i(1);
    % assign the right edge of the current face as the upper bound for taps in the last column:
    tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(end)),4) = lc_bnds_i(2);
    % loop through intermediate tap rows and assign the midpoints as upper and lower bounds:
    for J = 1:length(tap_col_w)-1
        tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(J)),4) = tap_col_w(J)+0.5*(tap_col_w(J+1)-tap_col_w(J));
        tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(J+1)),3) = tap_col_w(J)+0.5*(tap_col_w(J+1)-tap_col_w(J));
    end
    % assign the lower edge of the current face as the lower bound for taps in the first row:
    tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(1)),5) = lc_bnds_i(3);
    % assign the upper edge of the current face as the uper bound for taps in the last row:
    tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(end)),6) = lc_bnds_i(4);
    % loop through intermediate tap rows and assign the midpoints as upper and lower bounds:
    for J = 1:length(tap_row_y)-1
        tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(J)),6) = tap_row_y(J)+0.5*(tap_row_y(J+1)-tap_row_y(J));
        tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(J+1)),5) = tap_row_y(J)+0.5*(tap_row_y(J+1)-tap_row_y(J));
    end
    
end

end