function Q = charge_on_conductor(electrostatic_result, conductor_edges, show_geometry)
% CHARGE_ON_CONDUCTOR Find charge on the specified conductor.
%
% To use this function, firstly solve your problem geometry using the
% electrostatics solver in the Matlab PDE toolbox. Then you can use this
% function to find the charge on each conductor.
%
% This function calculates the enclosed charge by using Gauss's law. It
% evaluates the electric displacement directly on the mesh nodes along the
% boundary.
% 
% Q = charge_on_conductor(R, edges) calculates the charge enclosed within
% the region specified by the list of edges, using the electrostatics
% result R. The edges specified here must be one of the conductors for
% which you specified a voltage boundary condition when setting up the
% problem. You must specify a list of edges that forms a closed polygon, or
% else the result will be incorrect.
% 
% charge_on_conductor(..., true) produces a plot showing the geometry of
% the conductor and the direction of the normal vector (for checking).

if nargin < 3
    show_geometry = false;
end

if show_geometry
    figure();
end

Q = 0;
% For each edge that makes up the conductor
for edge_id = conductor_edges
    % find the mesh node IDs that line along this edge
    node_ids = findNodes(electrostatic_result.Mesh, "region", "edge", edge_id);

    % find their coordinates in (x,y) space
    node_coords = electrostatic_result.Mesh.Nodes(:, node_ids);

    % calculate the distance of each point from the start of the edge
    node_distances = sqrt(sum((node_coords - node_coords(:,1)).^2));

    % sort the points in order of distance
    [node_distances, sort_indices] = sort(node_distances);
    node_coords = node_coords(:, sort_indices);

    % calculate the unit vector that lies along the edge
    first = node_coords(:,1);
    last = node_coords(:,end);
    edge_unit = (last - first) ./ norm(last - first);

    % Calculate the direction of the normal to the edge. First rotate the
    % edge unit by 90 degrees, then check that stepping in that direction
    % goes outside the conductor. If it goes inside the conductor, flip the
    % sign.
    len = norm(last - first);
    centre = (first + last)/2;
    normal = [0 -1; 1 0] * edge_unit;
    test_point = centre + normal*len*1e-5;
    if isnan(electrostatic_result.interpolateElectricPotential(test_point(1), test_point(2)))
        % we are inside the conductor. try flipping the direction of the
        % normal
        normal = -normal;
        test_point = centre + normal*len*1e-5;
        if isnan(electrostatic_result.interpolateElectricPotential(test_point(1), test_point(2)))
            % still can't find a point outside of the conductor. This
            % isn't working
            error('Cannot find the normal to the surface. You must specify edge IDs that refer to conductors (not insulating regions).');
        end
    end

    % calculate the electric flux dotted with the normal
    Dx = electrostatic_result.ElectricFluxDensity.Dx(node_ids);
    Dy = electrostatic_result.ElectricFluxDensity.Dy(node_ids);
    D = normal(1)*Dx + normal(2)*Dy;

    % perform the line integral
    Q = Q + trapz(node_distances, D);

    if show_geometry
        plot([first(1) last(1)], [first(2) last(2)], 'b');
        hold on;
        plot(test_point(1), test_point(2), '*');
        quiver(centre(1), centre(2), normal(1)*len, normal(2)*len);
        text(centre(1), centre(2), sprintf("%i", edge_id));
    end
end

if show_geometry
    xlabel('x');
    ylabel('y');
end

end
