function modifiedSigma = computeFarFieldBCCoefficients(self)
    % Function that computes the modification to the coefficients necessary
    % to be consistent with the inclusion of a PML with the given
    % parameters
    % NOTE: this code assumes the edges with far field conditions are
    % straight lines and faces are planes that are parallel to the
    % coordinate axes.
    % There is some basic checking to that effect, but it is not 
    % exhaustive.

    %   Copyright 2020-2021 The MathWorks, Inc.


    modifiedSigma = @(loc,state) 0;
    geo = self.Geometry;
    % get approximate center of geometry
    center = mean([min(geo.Vertices) ; max(geo.Vertices)]);
    tol = sqrt(eps( norm(max(geo.Vertices) - min(geo.Vertices)) ));
    [p,e,~] = meshToPet(self.Mesh);
    dimensions = ["x","y","z"];
    for bcIndex = 1:length(self.BoundaryConditions.BCAssignments)
        bcRec = self.BoundaryConditions.BCAssignments(bcIndex);
        if isempty(bcRec.FarField)
            continue;
        end
        thickness = bcRec.Thickness;
        scaling = bcRec.Scaling;
        exponent = bcRec.Exponent;
        % generic modification function
        step = @(v) (v > 0).*( scaling*abs( v ).^exponent );
        for rIDs = 1:length(bcRec.RegionID)
            facet = bcRec.RegionID(rIDs);
            validateNonInteriorFacet(self,facet);
            % determine face orientation
            if self.IsTwoD
                nodes = verticesAttachedToEdges(geo,p,e,facet);
                dim = find(abs(nodes(:,1) - nodes(:,2)) < tol);
            else
                nodes = verticesAttachedToFace(self.Mesh,facet);
                d1 = nodes(:,3) - nodes(:,1);
                d2 = nodes(:,3) - nodes(:,2);
                dim = find(abs(d1) + abs(d2) < tol);
            end
            % error if dim is not found or multiple dims found (this
            % indicates an invalid geometry)
            if isempty(dim) || length(dim) > 1
                if self.IsTwoD
                    error(message('pde:ElectromagneticModel:invalidAbsorbingBCGeo2D'))
                else
                    error(message('pde:ElectromagneticModel:invalidAbsorbingBCGeo3D'));
                end
            end
            orientation = sign(nodes(dim,1) - center(dim));
            % PML should be zero at dim == offset
            offset = nodes(dim,1) - orientation * thickness;
            % update sigma based on facet info
            coord = dimensions(dim);
            if orientation < 0
                modifiedSigma = @(loc,state) modifiedSigma(loc,state) + ...
                    step(-loc.(coord) + offset);
            else
                modifiedSigma = @(loc,state) modifiedSigma(loc,state) + ...
                    step(loc.(coord) - offset);
            end
        end
    end
end


% find vertices attached to edges
function vae=verticesAttachedToEdges(geometry,p,e,edgeID)

% Example 1: with function handle
% pdem=createpde;
% gm = pdem.geometryFromEdges(@lshapeg);
% vae=verticesAttachedToEdges(gm,1);
%
% Example 2: with decsg
% model = createpde;
% gm2 = geometryFromEdges(model,dg)
% figure; pdegplot(gm2,'VertexLabels','on','EdgeLabels','on','FaceLabels','on')
% vae=verticesAttachedToEdges(gm2,7)
%

    if isequal(class(geometry),'pde.AnalyticGeometry') ...
            && ~isa(geometry.geom,'function_handle')
        dsgm = geometry.geom;
        sx = dsgm(2,edgeID)';
        ex = dsgm(3,edgeID)';
        sy = dsgm(4,edgeID)';
        ey = dsgm(5,edgeID)';
        vae = [sx sy; ex ey]';
    else % discrete geometry
        ind = find(e(5,:)==edgeID);
        sx = p(1,e(1,ind(1))); sy = p(2,e(1,ind(1))); 
        ex = p(1,e(1,ind(end))); ey = p(2,e(1,ind(end))); 
        vae = [sx sy ; ex ey]';
    end
end

function vaf = verticesAttachedToFace(mesh,faceID)
    % get node coords on face
    ind = findNodes(mesh,'region','Face',faceID);
    vaf = mesh.Nodes(:,ind(1:3)); 
end

function validateNonInteriorFacet(self,facet)
    % check if an edge is an interior edge and error if so
    region = "edge";
    if ~self.IsTwoD
        region = "face";
    end
    % get nodes on the facet
    nodesOnFacet = findNodes(self.Mesh,'region',region,facet);
    % check if nodes are also on the boundary
    if any(~isBoundaryNode(self.Mesh,nodesOnFacet))
        error(message('pde:ElectromagneticModel:interiorEdge'));
    end
end