FUNCTION enclose_pregion_in_pregion
(prgn : polar_complex_number_region, centre : complex_number_literal) : polar_complex_number_region;
FUNCTION angle(a : REAL) : REAL;
REPEAT WHILE a > PI; a := a - 2.0*PI; END_REPEAT;
REPEAT WHILE a <= -PI; a := a + 2.0*PI; END_REPEAT;
RETURN (a);
END_FUNCTION;
-- Find proper limits FOR direction interval
PROCEDURE angle_range(VAR amin, amax : REAL);
amin := angle(amin);
IF amin = PI THEN amin := -PI; END_IF;
amax := angle(amax);
IF amax <= amin THEN amax := amax + 2.0*PI; END_IF;
END_PROCEDURE;
-- Determine whether a direction is strictly within a direction interval
FUNCTION strictly_in(a : REAL;
aitv : finite_real_interval) : LOGICAL;
a := angle(a);
RETURN ({aitv.min < a < aitv.max} OR {aitv.min < a+2.0*PI < aitv.max});
END_FUNCTION;
-- Find min AND max AND related inclusion booleans among four candidates,
-- using a base direction chosen TO ensure the algebraic comparisons are valid.
PROCEDURE find_aminmax( ab,a0,a1,a2,a3 : REAL;
in0,in1,in2,in3 : BOOLEAN;
VAR amin,amax : REAL;
VAR amin_in,amax_in : BOOLEAN);
LOCAL
a : REAL;
END_LOCAL;
amin := angle(a0-ab); amin_in := in0;
amax := amin; amax_in := in0;
a := angle(a1-ab);
IF a = amin THEN amin_in := amin_in OR in1; END_IF;
IF a < amin THEN amin := a; amin_in := in1; END_IF;
IF a = amax THEN amax_in := amax_in OR in1; END_IF;
IF a > amax THEN amax := a; amax_in := in1; END_IF;
a := angle(a2-ab);
IF a = amin THEN amin_in := amin_in OR in2; END_IF;
IF a < amin THEN amin := a; amin_in := in2; END_IF;
IF a = amax THEN amax_in := amax_in OR in2; END_IF;
IF a > amax THEN amax := a; amax_in := in2; END_IF;
a := angle(a3-ab);
IF a = amin THEN amin_in := amin_in OR in3; END_IF;
IF a < amin THEN amin := a; amin_in := in3; END_IF;
IF a = amax THEN amax_in := amax_in OR in3; END_IF;
IF a > amax THEN amax := a; amax_in := in3; END_IF;
amin := amin+ab;
amax := amax+ab;
angle_range(amin,amax);
END_PROCEDURE;
LOCAL
ritp, ritv : real_interval;
aitp, aitv : finite_real_interval;
xp, yp, xc, yc, rmax, rmin, amin, amax, rc, acp, apc : REAL := 0.0;
rmax_in, rmin_in, amin_in, amax_in : BOOLEAN := FALSE;
rmxp, rmnp, x, y, r, a, ab, r0, a0, r1, a1, r2, a2, r3, a3 : REAL := 0.0;
in0, in1, in2, in3, inn : BOOLEAN := FALSE;
minclo, maxclo : open_closed := open;
END_LOCAL;
-- Extract elementary input information
IF NOT EXISTS (prgn) OR NOT EXISTS (centre) THEN RETURN (?); END_IF;
xp := prgn.centre.real_part;
yp := prgn.centre.imag_part;
ritp := prgn.distance_constraint;
aitp := prgn.direction_constraint;
xc := centre.real_part;
yc := centre.imag_part;
IF (xc = xp) AND (yc = yp) THEN RETURN (prgn); END_IF;
rc := SQRT((xp-xc)**2 + (yp-yc)**2);
acp := atan2(yp-yc,xp-xc);
apc := atan2(yc-yp,xc-xp);
rmnp := real_min(ritp);
-- Analyse cases BY existence OF max distance AND direction limits
IF max_exists(ritp) THEN
rmxp := real_max(ritp);
IF aitp.max - aitp.min = 2.0*PI THEN
-- annulus OR disk, WITH OR without slot OR puncture
inn := NOT max_included(aitp); -- slot exists;
a := angle(aitp.min); -- slot direction
rmax := rc+rmxp; rmax_in := max_included(ritp);
IF inn AND (acp = a) THEN rmax_in := FALSE; END_IF;
IF rc > rmxp THEN
a0 := ASIN(rmxp/rc);
amin := angle(acp-a0); amin_in := max_included(ritp);
IF amin = PI THEN amin := -PI; END_IF;
amax := angle(acp+a0); amax_in := amin_in;
IF amax < amin THEN amax := amax + 2.0*PI; END_IF;
rmin := rc-rmxp; rmin_in := amin_in;
IF inn THEN
-- slotted case
IF apc = a THEN rmin_in := FALSE; END_IF;
IF angle(amin+0.5*PI) = a THEN amin_in := FALSE; END_IF;
IF angle(amax-0.5*PI) = a THEN amax_in := FALSE; END_IF;
END_IF;
ELSE IF rc = rmxp THEN
amin := angle(acp-0.5*PI); amin_in := FALSE;
IF amin = PI THEN amin := -PI; END_IF;
amax := angle(acp+0.5*PI); amax_in := FALSE;
IF amax < amin THEN amax := amax + 2.0*PI; END_IF;
rmin := 0.0; rmin_in := max_included(ritp);
IF inn AND (apc = a) THEN rmin_in := FALSE; END_IF;
ELSE IF rc > rmnp THEN
IF inn AND (apc = a) THEN -- IN the slot
rmin := 0.0; rmin_in := FALSE;
amin := aitp.min; amin_in := FALSE;
amax := aitp.max; amax_in := FALSE;
ELSE
rmin := 0.0; rmin_in := TRUE;
amin := -PI; amin_in := FALSE;
amax := PI; amax_in := TRUE;
END_IF;
ELSE
rmin := rmnp-rc; rmin_in := min_included(ritp);
amin := -PI; amin_in := FALSE;
amax := PI; amax_in := TRUE;
IF inn THEN -- Special cases when aligned WITH slot
IF apc = a THEN
rmin_in := FALSE;
amin := aitp.min; amin_in := FALSE;
amax := aitp.max; amax_in := FALSE;
ELSE IF acp = a THEN
amin := aitp.min; amin_in := FALSE;
amax := aitp.max; amax_in := FALSE;
END_IF; END_IF;
END_IF;
END_IF; END_IF; END_IF;
ELSE -- direction range < 2*PI
-- Compute data FOR corners WITH respect TO xc,yc
x := xp + rmxp*cos(aitp.min) - xc;
y := yp + rmxp*sin(aitp.min) - yc;
r0 := SQRT(x**2 + y**2);
in0 := max_included(ritp) AND min_included(aitp);
IF r0 <> 0.0 THEN a0 := atan2(y,x); END_IF;
x := xp + rmxp*cos(aitp.max) - xc;
y := yp + rmxp*sin(aitp.max) - yc;
r1 := SQRT(x**2 + y**2);
in1 := max_included(ritp) AND max_included(aitp);
IF r1 <> 0.0 THEN a1 := atan2(y,x); END_IF;
x := xp + rmnp*cos(aitp.max) - xc;
y := yp + rmnp*sin(aitp.max) - yc;
r2 := SQRT(x**2 + y**2);
in2 := min_included(ritp) AND max_included(aitp);
IF r2 <> 0.0 THEN a2 := atan2(y,x); ELSE a2 := a1; in2 := in1; END_IF;
IF r1 = 0.0 THEN a1 := a2; in1 := in2; END_IF;
x := xp + rmnp*cos(aitp.min) - xc;
y := yp + rmnp*sin(aitp.min) - yc;
r3 := SQRT(x**2 + y**2);
in3 := min_included(ritp) AND min_included(aitp);
IF r3 <> 0.0 THEN a3 := atan2(y,x); ELSE a3 := a0; in3 := in0; END_IF;
IF r0 = 0.0 THEN a0 := a3; in0 := in3; END_IF;
IF rmnp = 0.0 THEN in2 := min_included(ritp); in3 := in2; END_IF;
IF (apc = angle(aitp.min)) OR (acp = angle(aitp.min)) THEN
in0 := min_included(aitp);
in3 := in0;
ELSE IF (apc = angle(aitp.max)) OR (acp = angle(aitp.max)) THEN
in1 := max_included(aitp);
in2 := in1;
END_IF; END_IF;
-- Find rmax
IF strictly_in(acp,aitp) THEN
rmax := rc+rmxp; rmax_in := max_included(ritp);
ELSE
rmax := r0; rmax_in := in0;
IF rmax = r1 THEN rmax_in := rmax_in OR in1; END_IF;
IF rmax < r1 THEN rmax := r1; rmax_in := in1; END_IF;
IF rmax = r2 THEN rmax_in := rmax_in OR in2; END_IF;
IF rmax < r2 THEN rmax := r2; rmax_in := in2; END_IF;
IF rmax = r3 THEN rmax_in := rmax_in OR in3; END_IF;
IF rmax < r3 THEN rmax := r3; rmax_in := in3; END_IF;
END_IF;
-- Find rmin
IF strictly_in(apc,aitp) THEN
IF rc >= rmxp THEN
rmin := rc-rmxp; rmin_in := max_included(ritp);
ELSE IF rc <= rmnp THEN
rmin := rmnp-rc; rmin_in := min_included(ritp);
ELSE
rmin := 0.0; rmin_in := TRUE;
END_IF; END_IF;
ELSE
rmin := r0; rmin_in := in0;
a := apc-aitp.min;
r := rc*COS(a);
IF {rmnp < r < rmxp} THEN -- USE nearest point on line segment
rmin := rc*SIN(ABS(a)); rmin_in := min_included(aitp);
END_IF;
a := apc-aitp.max;
r := rc*COS(a);
IF {rmnp < r < rmxp} THEN -- try nearest point on line segment
r := rc*SIN(ABS(a)); inn := max_included(aitp);
IF r = rmin THEN rmin_in := rmin_in OR inn; END_IF;
IF r < rmin THEN rmin := r; rmin_in := inn; END_IF;
END_IF;
IF r1 = rmin THEN rmin_in := rmin_in OR in1; END_IF;
IF r1 < rmin THEN rmin := r1; rmin_in := in1; END_IF;
IF r2 = rmin THEN rmin_in := rmin_in OR in2; END_IF;
IF r2 < rmin THEN rmin := r2; rmin_in := in2; END_IF;
IF r3 = rmin THEN rmin_in := rmin_in OR in3; END_IF;
IF r3 < rmin THEN rmin := r3; rmin_in := in3; END_IF;
END_IF;
-- Find amin AND amax, initially WITH respect TO base direction ab.
IF rc >= rmxp THEN -- outside outer circle
ab := acp;
find_aminmax(ab,a0,a1,a2,a3,in0,in1,in2,in3,amin,amax,amin_in,amax_in);
a := ACOS(rmxp/rc);
IF strictly_in(apc-a,aitp) THEN
amin := ab-ASIN(rmxp/rc); amin_in := max_included(ritp);
END_IF;
IF strictly_in(apc+a,aitp) THEN
amax := ab+ASIN(rmxp/rc); amax_in := max_included(ritp);
END_IF;
angle_range(amin,amax);
ELSE IF rc > rmnp THEN
ab := angle(0.5*(aitp.min+aitp.max)); -- REFERENCE direction
find_aminmax(ab,a0,a1,a2,a3,in0,in1,in2,in3,amin,amax,amin_in,amax_in);
ELSE
-- Using base direction midway IN prgn, compute all directions using
-- values which ensure a3 < a2 AND a0 < a1 algebraically.
ab := angle(0.5*(aitp.min+aitp.max)); -- REFERENCE direction
a0 := angle(a0-ab);
a1 := angle(a1-ab);
a2 := angle(a2-ab);
a3 := angle(a3-ab);
IF a3 > a2 THEN a2 := a2 + 2.0*PI; END_IF;
IF a0 > a1 THEN a0 := a0 + 2.0*PI; END_IF;
IF a3 < a0 THEN amin := a3; amin_in := in3;
ELSE amin := a0; amin_in := in0; END_IF;
IF a2 > a1 THEN amax := a2; amax_in := in2;
ELSE amax := a1; amax_in := in1; END_IF;
IF (amax - amin > 2.0*PI) OR
((amax - amin = 2.0*PI) AND (amin_in OR amax_in)) THEN
-- Cannot see out
amin := -PI; amin_in := FALSE;
amax := PI; amax_in := TRUE;
ELSE
amin := amin + ab;
amax := amax + ab;
angle_range(amin,amax);
END_IF;
END_IF; END_IF;
END_IF;
IF rmin_in THEN minclo := closed; END_IF;
IF rmax_in THEN maxclo := closed; END_IF;
ritv := make_finite_real_interval(rmin,minclo,rmax,maxclo);
ELSE -- NOT max_exists(ritp)
IF (rc > rmnp) AND strictly_in(apc,aitp) THEN
RETURN (?); -- No pregion exists. (Would require whole plane.)
END_IF;
IF aitp.max - aitp.min = 2.0*PI THEN
-- complement OF disk, WITH OR without slot
a := angle(aitp.min); -- slot direction
IF rc > rmnp THEN -- already excluded IF NOT aligned WITH slot
IF max_included(aitp) THEN
RETURN (?); -- No pregion exists. (Would require whole plane.)
END_IF;
rmin := 0.0; rmin_in := FALSE;
amin := aitp.min; amin_in := FALSE;
amax := aitp.max; amax_in := FALSE;
ELSE
rmin := rmnp-rc; rmin_in := min_included(ritp);
amin := -PI; amin_in := FALSE;
amax := PI; amax_in := TRUE;
IF NOT max_included(aitp) THEN -- Special cases when aligned WITH slot
IF apc = a THEN
rmin_in := FALSE;
amin := aitp.min; amin_in := FALSE;
amax := aitp.max; amax_in := FALSE;
ELSE IF acp = a THEN
amin := aitp.min; amin_in := FALSE;
amax := aitp.max; amax_in := FALSE;
END_IF; END_IF;
END_IF;
END_IF;
ELSE -- direction range < 2*PI
-- Compute data FOR corners WITH respect TO xc,yc (two at infinity)
a0 := angle(aitp.min);
in0 := FALSE;
a1 := angle(aitp.max);
in1 := FALSE;
x := xp + rmnp*cos(aitp.max) - xc;
y := yp + rmnp*sin(aitp.max) - yc;
r2 := SQRT(x**2 + y**2);
in2 := min_included(ritp) AND max_included(aitp);
IF r2 <> 0.0 THEN a2 := atan2(y,x); ELSE a2 := a1; in2 := in1; END_IF;
x := xp + rmnp*cos(aitp.min) - xc;
y := yp + rmnp*sin(aitp.min) - yc;
r3 := SQRT(x**2 + y**2);
in3 := min_included(ritp) AND min_included(aitp);
IF r3 <> 0.0 THEN a3 := atan2(y,x); ELSE a3 := a0; in3 := in0; END_IF;
IF rmnp = 0.0 THEN in2 := min_included(ritp); in3 := in2; END_IF;
IF (apc = angle(aitp.min)) OR (acp = angle(aitp.min)) THEN
in0 := min_included(aitp);
in3 := in0;
ELSE IF (apc = angle(aitp.max)) OR (acp = angle(aitp.max)) THEN
in1 := max_included(aitp);
in2 := in1;
END_IF; END_IF;
-- Find rmin
IF strictly_in(apc,aitp) THEN
rmin := rmnp-rc; rmin_in := min_included(ritp);
ELSE
rmin := r2; rmin_in := in2;
a := apc-aitp.min;
r := rc*COS(a);
IF rmnp < r THEN -- USE nearest point on aitp.min ray
rmin := rc*SIN(ABS(a)); rmin_in := min_included(aitp);
END_IF;
a := apc-aitp.max;
r := rc*COS(a);
IF rmnp < r THEN -- try nearest point on aitp.max ray
r := rc*SIN(ABS(a)); inn := max_included(aitp);
IF r = rmin THEN rmin_in := rmin_in OR inn; END_IF;
IF r < rmin THEN rmin := r; rmin_in := inn; END_IF;
END_IF;
IF r3 = rmin THEN rmin_in := rmin_in OR in3; END_IF;
IF r3 < rmin THEN rmin := r3; rmin_in := in3; END_IF;
END_IF;
-- Find amin AND amax
ab := angle(0.5*(aitp.min+aitp.max)); -- REFERENCE direction
IF rc > rmnp THEN
find_aminmax(ab,a0,a1,a2,a3,in0,in1,in2,in3,amin,amax,amin_in,amax_in);
ELSE
-- Using base direction midway IN prgn, compute all directions using
-- values which ensure a3 < a2 AND a0 < a1 algebraically.
a0 := angle(a0-ab);
a1 := angle(a1-ab);
a2 := angle(a2-ab);
a3 := angle(a3-ab);
IF a3 > a2 THEN a2 := a2 + 2.0*PI; END_IF;
IF a0 > a1 THEN a0 := a0 + 2.0*PI; END_IF;
IF a3 < a0 THEN amin := a3; amin_in := in3;
ELSE amin := a0; amin_in := in0; END_IF;
IF a2 > a1 THEN amax := a2; amax_in := in2;
ELSE amax := a1; amax_in := in1; END_IF;
IF (amax - amin > 2.0*PI) OR
((amax - amin = 2.0*PI) AND (amin_in OR amax_in)) THEN
-- Cannot see out
amin := -PI; amin_in := FALSE;
amax := PI; amax_in := TRUE;
IF (rmin = 0.0) AND rmin_in THEN
RETURN (?); -- No pregion exists. (Would require whole plane.)
END_IF;
ELSE
amin := amin + ab;
amax := amax + ab;
angle_range(amin,amax);
END_IF;
END_IF;
END_IF;
IF rmin_in THEN minclo := closed; END_IF;
ritv := make_real_interval_from_min(rmin,minclo);
END_IF;
minclo := open; maxclo := open;
IF amin_in THEN minclo := closed; END_IF;
IF amax_in THEN maxclo := closed; END_IF;
aitv := make_finite_real_interval(amin,minclo,amax,maxclo);
-- Construct polar region
RETURN (make_polar_complex_number_region(centre,ritv,aitv)); END_FUNCTION; -- enclose_pregion_in_pregion
|