CREATE OR REPLACE FUNCTION ISRSrange_area(
    axis geometry,
    area geometry
) RETURNS geometry
AS $$
DECLARE
    area_subset geometry;
    result_geom geometry;
BEGIN
    -- In case area is a multipolygon, process the union of those
    -- polygons, which intersect with the axis. The union is to avoid
    -- problems with invalid/self-intersecting multipolygons
    SELECT ST_Union(a_dmp.geom)
        INTO STRICT area_subset
        FROM (SELECT ST_MakeValid(ST_Transform(geom, ST_SRID(axis)))
                FROM ST_Dump(area)) AS a_dmp (geom)
        WHERE ST_Intersects(a_dmp.geom, axis)
        HAVING ST_Union(a_dmp.geom) IS NOT NULL;

    WITH
        rotated_ends AS (
            SELECT ST_Collect(ST_Scale(
                    ST_Translate(e,
                        (ST_X(p1) - ST_X(p2)) / 2,
                        (ST_Y(p1) - ST_Y(p2)) / 2),
                    ST_Point(d, d), p1)) AS blade
                FROM (SELECT i, ST_PointN(axis, i) AS p1
                        FROM (VALUES (1), (-1)) AS idx (i)) AS ep,
                    ST_Rotate(ST_PointN(axis, i*2), pi()/2, p1) AS ep2 (p2),
                    ST_Makeline(p1, p2) AS e (e),
                    LATERAL (
                        SELECT (ST_MaxDistance(p1, area_subset) / ST_Length(e))
                            * 2) AS d (d)),
        range_area AS (
            -- Split area by orthogonal lines at the ends of the clipped axis
            SELECT (ST_Dump(ST_CollectionExtract(
                    ST_Split(area_subset, blade), 3))).geom
                FROM rotated_ends)
        -- From the polygons returned by the last CTE, select only those
        -- around the clipped axis
        SELECT ST_Multi(ST_CollectionExtract(ST_MakeValid(ST_Transform(
                ST_Union(range_area.geom), ST_SRID(area))), 3))
            INTO result_geom
            FROM range_area
            WHERE ST_Intersects(ST_Buffer(range_area.geom, -0.0001), axis);

    RETURN result_geom;
END;
    $$
    LANGUAGE plpgsql
    STABLE PARALLEL SAFE;