function [sx, sy] = hit_and_run_general(vx, vy, N, nBurn)
%HIT_AND_RUN_GENERAL  Hit-and-run for an arbitrary simple polygon (convex or
%   non-convex). Along each random line it finds *every* in-region segment and
%   samples uniformly across their union, so concavities are handled correctly
%   (unlike the convex-only chord in hit_and_run.m). Starts at the origin —
%   make_region's non-convex regions are star-shaped about it.
%
%   This runs in the numbl interpreter (sort + point-in-polygon tests don't
%   JIT), so it's used only for the non-convex demo at modest N.
nv = numel(vx);
px = 0;
py = 0;
total = nBurn + N;
sx = zeros(N, 1);
sy = zeros(N, 1);

for step = 1:total
    th = 2 * pi * rand;
    dx = cos(th);
    dy = sin(th);

    % All parameters t where the line p + t*d crosses the polygon boundary.
    ts = zeros(1, nv);
    m = 0;
    for i = 1:nv
        j = mod(i, nv) + 1;
        ex = vx(j) - vx(i);
        ey = vy(j) - vy(i);
        denom = dy * ex - dx * ey;
        if abs(denom) < 1e-12
            continue
        end
        wx = vx(i) - px;
        wy = vy(i) - py;
        sParam = (dx * wy - dy * wx) / denom;   % position along the edge
        if sParam >= 0 && sParam < 1
            m = m + 1;
            ts(m) = (wy * ex - wx * ey) / denom;  % position along the line
        end
    end
    if m < 2
        continue
    end
    ts = sort(ts(1:m));

    % In-region intervals are consecutive crossings whose midpoint is inside.
    totalLen = 0;
    for k = 1:m - 1
        tm = (ts(k) + ts(k + 1)) / 2;
        if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
            totalLen = totalLen + (ts(k + 1) - ts(k));
        end
    end
    if totalLen <= 0
        continue
    end

    % Pick a point uniformly across the union of in-region intervals.
    u = totalLen * rand;
    tpick = 0;
    for k = 1:m - 1
        tm = (ts(k) + ts(k + 1)) / 2;
        if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
            len = ts(k + 1) - ts(k);
            if u <= len
                tpick = ts(k) + u;
                break
            end
            u = u - len;
        end
    end
    px = px + tpick * dx;
    py = py + tpick * dy;

    if step > nBurn
        sx(step - nBurn) = px;
        sy(step - nBurn) = py;
    end
end
end

function inside = point_in_poly(x, y, vx, vy)
%POINT_IN_POLY  Ray-casting test for a point against polygon (VX, VY).
n = numel(vx);
inside = false;
j = n;
for i = 1:n
    if ((vy(i) > y) ~= (vy(j) > y)) && ...
            (x < (vx(j) - vx(i)) * (y - vy(i)) / (vy(j) - vy(i)) + vx(i))
        inside = ~inside;
    end
    j = i;
end
end