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;
  -- Determine whether a REAL  is strictly within a REAL  interval
  FUNCTION strictly_in(z    : REAL;
                       zitv : real_interval) : LOGICAL;
    RETURN  ((NOT min_exists(zitv) OR  (z > real_min(zitv))) AND
      (NOT max_exists(zitv) OR  (z < real_max(zitv))));
  END_FUNCTION;
  -- Include direction IN  minmax collection
  PROCEDURE  angle_minmax(    ab, a            : REAL;
                             a_in             : BOOLEAN;
                         VAR  amin, amax       : REAL;
                         VAR  amin_in, amax_in : BOOLEAN);
    a := angle(a - ab);
    IF  amin = a THEN   amin_in := amin_in OR  a_in;  END_IF;
    IF  amin > a THEN   amin := a;  amin_in := a_in;  END_IF;
    IF  amax = a THEN   amax_in := amax_in OR  a_in;  END_IF;
    IF  amax < a THEN   amax := a;  amax_in := a_in;  END_IF;
  END_PROCEDURE;
  -- Include distance IN  max collection
  PROCEDURE  range_max(    r       : REAL;
                          incl    : BOOLEAN;
                      VAR  rmax    : REAL;
                      VAR  rmax_in : BOOLEAN);
    IF  rmax = r THEN   rmax_in := rmax_in OR  incl;   END_IF;
    IF  rmax < r THEN   rmax := r;  rmax_in := incl;  END_IF;
  END_PROCEDURE;
  -- Include distance IN min collection
  PROCEDURE range_min(    r       : REAL;
                          incl    : BOOLEAN;
                      VAR  rmin    : REAL;
                      VAR rmin_in : BOOLEAN);
    IF  rmin = r THEN   rmin_in := rmin_in OR  incl;  END_IF;
    IF  (rmin < 0.0) OR  (rmin > r) THEN   rmin := r;  rmin_in := incl;  END_IF;
  END_PROCEDURE;
  LOCAL
    xitv, yitv : real_interval;
    is_xmin, is_xmax, is_ymin, is_ymax : BOOLEAN;
    xmin, xmax, ymin, ymax, xc, yc : REAL  := 0.0;
    xmin_in, xmax_in, ymin_in, ymax_in : BOOLEAN  := FALSE;
    rmin, rmax : REAL  := -1.0;
    amin : REAL  := 4.0;
    amax : REAL  := -4.0;
    rmax_exists, outside : BOOLEAN  := TRUE;
    rmin_in, rmax_in, amin_in, amax_in : BOOLEAN := FALSE;
    ab, a, r : REAL := 0.0;
    incl : BOOLEAN;
    ritv : real_interval;
    aitv : finite_real_interval;
    minclo, maxclo : open_closed := open;
  END_LOCAL;
  IF  NOT  EXISTS  (crgn) OR  NOT  EXISTS  (centre) THEN   RETURN  (?);  END_IF;
  -- Extract elementary input information
  xitv := crgn.real_constraint;
  yitv := crgn.imag_constraint;
  xc := centre.real_part;
  yc := centre.imag_part;
  is_xmin := min_exists(xitv);
  is_xmax := max_exists(xitv);
  is_ymin := min_exists(yitv);
  is_ymax := max_exists(yitv);
  IF  is_xmin THEN   xmin := real_min(xitv);  xmin_in := min_included(xitv);  END_IF;
  IF  is_xmax THEN   xmax := real_max(xitv);  xmax_in := max_included(xitv);  END_IF;
  IF  is_ymin THEN   ymin := real_min(yitv);  ymin_in := min_included(yitv);  END_IF;
  IF  is_ymax THEN   ymax := real_max(yitv);  ymax_in := max_included(yitv);  END_IF;
  rmax_exists := is_xmin AND  is_xmax AND  is_ymin AND  is_ymax;
  -- Identify base direction WITH respect TO  which all relevant directions lie
  -- within +/- 0.5*PI, OR that the centre lies properly inside crgn.
  IF        is_xmin AND  (xc <= xmin)  THEN   ab := 0.0;
  ELSE  IF   is_ymin AND  (yc <= ymin)  THEN   ab := 0.5*PI;
  ELSE  IF   is_ymax AND  (yc >= ymax)  THEN   ab := -0.5*PI;
  ELSE  IF   is_xmax AND  (xc >= xmax)  THEN   ab := PI;
  ELSE                                      outside := FALSE;
  END_IF;  END_IF;  END_IF;  END_IF;
  IF  NOT  outside AND  NOT  rmax_exists THEN
    RETURN  (?);  -- No enclosing polar region EXISTS (requires whole plane)
  END_IF;
  -- Identify any closest point on a side but NOT  a corner.
  IF        is_xmin AND  (xc <= xmin) AND  strictly_in(yc,yitv) THEN
    rmin := xmin - xc;  rmin_in := xmin_in;
  ELSE  IF   is_ymin AND  (yc <= ymin) AND  strictly_in(xc,xitv) THEN
    rmin := ymin - yc;  rmin_in := ymin_in;
  ELSE  IF   is_ymax AND  (yc >= ymax) AND  strictly_in(xc,xitv) THEN
    rmin := yc - ymax;  rmin_in := ymax_in;
  ELSE  IF   is_xmax AND  (xc >= xmax) AND  strictly_in(yc,yitv) THEN
    rmin := xc - xmax;  rmin_in := xmax_in;
  END_IF;  END_IF;  END_IF;  END_IF;
  IF  is_xmin THEN
    IF  is_ymin THEN   -- Consider lower left corner
      r := SQRT((xmin-xc)**2 + (ymin-yc)**2);
      incl := xmin_in AND  ymin_in;
      IF  rmax_exists THEN   range_max(r,incl,rmax,rmax_in);  END_IF;
      IF  outside THEN
        IF  r > 0.0 THEN
          range_min(r,incl,rmin,rmin_in);
          a := angle(atan2(ymin-yc,xmin-xc) - ab);
          IF  xc = xmin THEN   incl := xmin_in;  END_IF;
          IF  yc = ymin THEN   incl := ymin_in;  END_IF;
          angle_minmax(ab,a,incl,amin,amax,amin_in,amax_in);
        ELSE   -- Centre at lower left corner
          rmin := 0.0;                rmin_in := xmin_in AND  ymin_in;
          amin := angle(0.0-ab);      amin_in := ymin_in;
          amax := angle(0.5*PI-ab);   amax_in := xmin_in;
        END_IF;
      END_IF;
    ELSE  IF  xc <= xmin THEN   -- Consider points near (xmin, -infinity)
      angle_minmax(ab,-0.5*PI,(xc=xmin) AND  xmin_in,amin,amax,amin_in,amax_in);
    END_IF;  END_IF;
    IF  NOT  is_ymax AND  (xc <= xmin) THEN   -- Consider points near (xmin, +infinity)
      angle_minmax(ab,0.5*PI,(xc=xmin) AND  xmin_in,amin,amax,amin_in,amax_in);
    END_IF;
  END_IF;
  IF  is_ymin THEN
    IF  is_xmax THEN   -- Consider lower right corner
      r := SQRT((xmax-xc)**2 + (ymin-yc)**2);
      incl := xmax_in AND  ymin_in;
      IF  rmax_exists THEN   range_max(r,incl,rmax,rmax_in);  END_IF;
      IF  outside THEN
        IF  r > 0.0 THEN
          range_min(r,incl,rmin,rmin_in);
          a := angle(atan2(ymin-yc,xmax-xc) - ab);
          IF  xc = xmax THEN   incl := xmax_in;  END_IF;
          IF  yc = ymin THEN   incl := ymin_in;  END_IF;
          angle_minmax(ab,a,incl,amin,amax,amin_in,amax_in);
        ELSE   -- Centre at lower right corner
          rmin := 0.0;                rmin_in := xmax_in AND  ymin_in;
          amin := angle(0.5*PI-ab);   amin_in := ymin_in;
          amax := angle(PI-ab);       amax_in := xmax_in;
        END_IF;
      END_IF;
    ELSE  IF  yc <= ymin THEN   -- Consider points near (+infinity, ymin)
      angle_minmax(ab,0.0,(yc=ymin) AND  ymin_in,amin,amax,amin_in,amax_in);
    END_IF;  END_IF;
    IF  NOT  is_xmin AND  (yc <= ymin) THEN   -- Consider points near (-infinity, ymin)
      angle_minmax(ab,PI,(yc=ymin) AND  ymin_in,amin,amax,amin_in,amax_in);
    END_IF;
  END_IF;
  IF  is_xmax THEN
    IF  is_ymax THEN   -- Consider upper right corner
      r := SQRT((xmax-xc)**2 + (ymax-yc)**2);
      incl := xmax_in AND  ymax_in;
      IF  rmax_exists THEN   range_max(r,incl,rmax,rmax_in);  END_IF;
      IF  outside THEN
        IF  r > 0.0 THEN
          range_min(r,incl,rmin,rmin_in);
          a := angle(atan2(ymax-yc,xmax-xc) - ab);
          IF  xc = xmax THEN   incl := xmax_in;  END_IF;
          IF  yc = ymax THEN   incl := ymax_in;  END_IF;
          angle_minmax(ab,a,incl,amin,amax,amin_in,amax_in);
        ELSE   -- Centre at lower left corner
          rmin := 0.0;                rmin_in := xmax_in AND  ymax_in;
          amin := angle(-PI-ab);      amin_in := ymax_in;
          amax := angle(-0.5*PI-ab);  amax_in := xmax_in;
        END_IF;
      END_IF;
    ELSE  IF  xc >= xmax THEN   -- Consider points near (xmax, +infinity)
      angle_minmax(ab,0.5*PI,(xc=xmax) AND  xmax_in,amin,amax,amin_in,amax_in);
    END_IF;  END_IF;
    IF  NOT  is_ymin AND  (xc >= xmax) THEN   -- Consider points near (xmax, -infinity)
      angle_minmax(ab,-0.5*PI,(xc=xmax) AND  xmax_in,amin,amax,amin_in,amax_in);
    END_IF;
  END_IF;
  IF  is_ymax THEN
    IF  is_xmin THEN   -- Consider upper left corner
      r := SQRT((xmin-xc)**2 + (ymax-yc)**2);
      incl := xmin_in AND  ymax_in;
      IF  rmax_exists THEN   range_max(r,incl,rmax,rmax_in);  END_IF;
      IF  outside THEN
        IF  r > 0.0 THEN
          range_min(r,incl,rmin,rmin_in);
          a := angle(atan2(ymax-yc,xmin-xc) - ab);
          IF  xc = xmin THEN   incl := xmin_in;  END_IF;
          IF  yc = ymax THEN   incl := ymax_in;  END_IF;
          angle_minmax(ab,a,incl,amin,amax,amin_in,amax_in);
        ELSE   -- Centre at lower right corner
          rmin := 0.0;                rmin_in := xmin_in AND  ymax_in;
          amin := angle(0.5*PI-ab);   amin_in := ymax_in;
          amax := angle(PI-ab);       amax_in := xmin_in;
        END_IF;
      END_IF;
    ELSE IF  yc >= ymax THEN   -- Consider points near (-infinity, ymax)
      angle_minmax(ab,PI,(yc=ymax) AND  ymax_in,amin,amax,amin_in,amax_in);
    END_IF;  END_IF;
    IF  NOT is_xmax AND  (yc >= ymax) THEN   -- Consider points near (+infinity, ymax)
      angle_minmax(ab,0.0,(yc=ymax) AND ymax_in,amin,amax,amin_in,amax_in);
    END_IF;
  END_IF;
  IF  outside THEN   -- Change direction origin FROM ab back TO zero
    amin := angle(amin+ab);
    IF  amin = PI THEN   amin := -PI;  END_IF;
    amax := angle(amax+ab);
    IF  amax <= amin THEN   amax := amax + 2.0*PI;  END_IF;
  ELSE
    amin := -PI;  amin_in := FALSE;
    amax := PI;   amax_in := FALSE;
  END_IF;
  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);
  minclo := open;
  IF  rmin_in THEN   minclo := closed;  END_IF;
  IF  rmax_exists THEN
    maxclo := open;
    IF rmax_in THEN  maxclo := closed;  END_IF;
    ritv := make_finite_real_interval(rmin,minclo,rmax,maxclo);
  ELSE
    ritv := make_real_interval_from_min(rmin,minclo);
  END_IF;
  RETURN (make_polar_complex_number_region(centre,ritv,aitv));