function [BClhs,iremove] = get_BC_lhs(rhogrid,BCrho,tt)
% get left hand side for boundary conditions and index of equation that
% should be removed

% #codegen

coder.extrinsic('error');

BClhs = zeros(1,size(tt.Lam,2)); % init for size
BCrho_eq_rho = (BCrho == rhogrid);
if any(BCrho_eq_rho)
    % rho exactly at rhogrid point, take splines already stored
    BClhs(:) = tt.Lam(BCrho_eq_rho,:);
else
    % rho at intermediate point, must evaluate splines
    % derive spline order from rho points
    order = tt.nsp-numel(rhogrid)+2;
    BClhs(:) = tt.scal*splineval(BCrho,rhogrid,order,[1 0 0])';
end

iremove = tt.nsp; % must always remove last equation since we didn't 
% calculate the PDE boudnary terms from the partial integration!

end