http://mathling.com/geometric/hyperbolic  library module

Imports

import module namespace path="http://mathling.com/geometric/path"
       at "../geo/path.xqy"

import module namespace edge="http://mathling.com/geometric/edge"
       at "../geo/edge.xqy"

import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy"

import module namespace ellipse="http://mathling.com/geometric/ellipse"
       at "../geo/ellipse.xqy"

import module namespace box="http://mathling.com/geometric/rectangle"
       at "../geo/rectangle.xqy"

import module namespace geom="http://mathling.com/geometric"
       at "../geo/euclidean.xqy"

import module namespace config="http://mathling.com/core/config"
       at "../core/config.xqy"

import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy"

import module namespace point="http://mathling.com/geometric/point"
       at "../geo/point.xqy"

Variables

Functions

declare function this:circle-through(
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)?
{
  let $ax := point:px($a)
  let $ay := point:py($a)
  let $bx := point:px($b)
  let $by := point:py($b)
  let $axby-aybx := $ax * $by - $ay * $bx
  where abs($axby-aybx) >= $this:TOLERANCE
  return (
    let $cx :=
      -0.5 *
        ($ay * $bx * $bx + $ay * $by * $by -
          ($ax * $ax + $ay * $ay + 1) * $by + $ay
        )
        div $axby-aybx
    let $cy :=
      0.5 *
        ($ax * $bx * $bx + $ax * $by * $by -
          ($ax * $ax + $ay * $ay + 1) * $bx + $ax
        )
        div $axby-aybx
    let $radius :=
      abs(
        -0.5 * math:sqrt(
          $ax * $ax * $ax * $ax * $bx * $bx +     
          $ax * $ax * $ax * $ax * $by * $by -     
          2 * $ax * $ax * $ax * $bx * $bx * $bx - 
          2 * $ax * $ax * $ax * $bx * $by * $by + 
          2 * $ax * $ax * $ay * $ay * $bx * $bx + 
          2 * $ax * $ax * $ay * $ay * $by * $by - 
          2 * $ax * $ax * $ay * $bx * $bx * $by - 
          2 * $ax * $ax * $ay * $by * $by * $by + 
          $ax * $ax * $bx * $bx * $bx * $bx +     
          2 * $ax * $ax * $bx * $bx * $by * $by + 
          $ax * $ax * $by * $by * $by * $by -     
          2 * $ax * $ay * $ay * $bx * $bx * $bx - 
          2 * $ax * $ay * $ay * $bx * $by * $by +
          $ay * $ay * $ay * $ay * $bx * $bx +
          $ay * $ay * $ay * $ay * $by * $by -
          2 * $ay * $ay * $ay * $bx * $bx * $by -
          2 * $ay * $ay * $ay * $by * $by * $by +
          $ay * $ay * $bx * $bx * $bx * $bx +
          2 * $ay * $ay * $bx * $bx * $by * $by +
          $ay * $ay * $by * $by * $by * $by -
          2 * $ax * $ax * $ax * $bx -
          2 * $ax * $ax * $ay * $by +
          4 * $ax * $ax * $bx * $bx -
          2 * $ax * $ay * $ay * $bx +
          8 * $ax * $ay * $bx * $by -
          2 * $ax * $bx * $bx * $bx -
          2 * $ax * $bx * $by * $by -
          2 * $ay * $ay * $ay * $by +
          4 * $ay * $ay * $by * $by -
          2 * $ay * $bx * $bx * $by -
          2 * $ay * $by * $by * $by +
          $ax * $ax - 2 * $ax * $bx +
          $ay * $ay - 2 * $ay * $by +
          $bx * $bx + $by * $by
        )
        div $axby-aybx
      )
    return ellipse:circle(point:point($cx, $cy), $radius)
  )
}

declare function this:arc-through(
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $ax := point:px($a)
  let $ay := point:py($a)
  let $bx := point:px($b)
  let $by := point:py($b)
  let $circle := this:circle-through($a, $b)
  return (
    if (empty($circle)) then ( (: near-aligned: straight line :)
      edge:edge($a, $b)
    ) else (
      let $c := ellipse:center($circle)
      let $r := ellipse:radius($circle)
      let $θ := -math:atan2(-$ay + point:py($c), $ax - point:px($c))
      let $θ2 := -math:atan2(-$by + point:py($c), $bx - point:px($c))
      (: Recenter the angles to [0, 2π] :)
      let $θ := util:remap-radians($θ)      
      let $θ2 := util:remap-radians($θ2)      
      let $delta := abs($θ - $θ2)
      return (
        if (($θ > $θ2 and $delta > math:pi()) or ($θ < $θ2 and $delta < math:pi()))
        then edge:arc($c, $r, $a, $b, false())
        else edge:arc($c, $r, $a, $b, true())
      )
    )
  )
}

declare function this:edge(
  $scale as xs:numeric,
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:arc-through($a, $b)=>geom:scale($scale)
}

declare function this:line(
  $scale as xs:numeric,
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:edge($scale, $start, $end)
}

declare function this:shape(
  $scale as xs:numeric,
  $points as map(xs:string,item()*)*
) as map(xs:string,item()*)
{
  this:polygon($scale, $points)
}

declare function this:box(
  $scale as xs:numeric,
  $box as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $points := (
    point:point(box:min-px($box),box:min-py($box)),
    point:point(box:min-px($box),box:max-py($box)),
    point:point(box:max-px($box),box:max-py($box)),
    point:point(box:max-px($box),box:min-py($box)),
    point:point(box:min-px($box),box:min-py($box))
  )
  return this:polygon($scale, $points)
}

declare function this:triangle(
  $scale as xs:numeric,
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*),
  $c as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:polygon( $scale, ($a, $b, $c) )
}

declare function this:supporting-line(
  $scale as xs:numeric,
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $ax := point:px($a)
  let $ay := point:py($a)
  let $bx := point:px($b)
  let $by := point:py($b)
  let $circle := this:circle-through($a, $b) (: Compute circle parameters :)
  return (
    if (empty($circle)) then (
      (: Project a and b points onto the circle :)
      let $θ := math:atan2($ay, $ax)
      let $θ2 := math:atan2($by, $bx)
      return (
        edge:edge(
          point:point($scale*math:cos($θ), $scale*math:sin($θ)),
          point:point($scale*math:cos($θ2), $scale*math:sin($θ2))
        )
      )
    ) else (
      ellipse:circle($point:ORIGIN, $scale*ellipse:radius($circle))
    )
  )
}

declare function this:all-support-lines(
  $scale as xs:numeric,
  $pts as map(xs:string,item()*)*
) as map(xs:string,item()*)*
{
  for $p at $i in $pts
  let $next := 
    if (empty($pts[$i + 1]))
    then $pts[1]
    else $pts[$i + 1]
  return
    this:supporting-line($scale, $p, $next)
}

declare function this:polygon(
  $scale as xs:numeric,
  $pts as map(xs:string,item()*)*
) as map(xs:string,item()*)
{
  let $edges := 
    for $p at $i in $pts
    let $next := 
      if (empty($pts[$i + 1]))
      then $pts[1]
      else $pts[$i + 1]
    where not(this:within-tolerance($p=>geom:scale($scale), $next=>geom:scale($scale), $this:TOLERANCE))
    return
      this:edge($scale, $p, $next)
  return (
    path:polygon($edges)
  )
}

declare function this:distance(
  $scale as xs:numeric,
  $p as map(xs:string,item()*),
  $q as map(xs:string,item()*)
) as xs:double
{
  let $pqd := geom:distance($p,$q)
  let $pd := geom:distance(point:point(0,0),$p)
  let $qd := geom:distance(point:point(0,0),$q)
  let $x :=  1 + 2*($pqd*$pqd)*$scale*$scale div (($scale*$scale - $pd*$pd)*($scale*$scale - $qd*$qd))
  (: acosh(x) = ln(x + sqrt(x^2 - 1)) :)
  return (
    math:log($x + math:sqrt($x*$x - 1))
  )
}

declare function this:tesselate(
  $scale as xs:numeric,
  $n as xs:integer,
  $k as xs:integer,
  $layer as xs:integer
) as map(xs:string,item()*)*
{
  if (($n - 2)*($k - 2) <= 4)
  then errors:error("GEOM-INVHYPER", ($n - 2) * ($k - 2))
  else ()
  ,
  this:remove-duplicates(
    $scale,
    let $center := this:center-polygon($n, $k)
    let $rule := 0
    return (
      $center,
      this:replicate-polygon($center, $rule, $n, $k, $layer - 1)
    )
  )
}

declare function this:sorted-hypers(
  $size as xs:numeric,
  $polys as map(xs:string,item()*)*
) as map(xs:string,item()*)*
{
  (: Sorted and turned into hyperpolygons :)
  let $origin := point:point(0,0)
  for $poly in $polys
  let $bounds := geom:bounding-box($poly=>geom:scale($size))
  let $center := geom:region-center($bounds)
  let $distance := round(this:distance($size, $origin, $center))
  let $angle := round(geom:path-angle($origin, $center))
  order by $distance ascending, $angle ascending
  return (
    this:polygon($size, geom:vertices($poly))=>
      map:put("layer", $distance)=>
      map:put("angle", $angle)
  )
}

declare function this:reflect(
  $pt as map(xs:string,item()*),
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  (: Find a unit vector $uv in the direction of the line :)
  (: $axby-bxay = ax*by - bx*ay :)
  let $axby-bxay := point:px($a) * point:py($b) - point:px($b) * point:py($a)
  return (
    if (false()) then (
      (: Euclidean reflection = 
         Ref(l,v) = 2 (v·l)/(l·l) - v = 2 Proj(l,v) - v
       :)
      (: distance between a and b √((ax - bx)² + (ax - by)²) :)
      let $delta :=
        math:sqrt(
          (point:px($a) - point:px($b))*(point:px($a) - point:px($b)) +
          (point:py($a) - point:py($b))*(point:py($a) - point:py($b))
        )
      let $uv :=
        point:point(
          (point:px($b) - point:px($a)) div $delta,
          (point:py($b) - point:py($a)) div $delta
        )
      (: 2 * ((px - ax)*ux + (py-ay)*uy) :)
      let $factor :=
        2.0 * (
          (point:px($pt) - point:px($a)) * point:px($uv) +
          (point:py($pt) - point:py($a)) * point:py($uv)
        )
      return
        (: 2 * $a + $factor * $uv - $p :)
        point:point(
          2.0 * point:px($a) + $factor * point:px($uv) - point:px($pt),
          2.0 * point:py($a) + $factor * point:py($uv) - point:py($pt)
        )
    ) else if (abs($axby-bxay) < $this:TOLERANCE) then ( (: straight :)
      geom:reflect($pt, $a, $b)
    ) else (
      (: Find the center of the circle through these points :)
      (:  
       :  s1 = (xa² + ya² + 1)/2
       :  s2 = (xb² + yb² + 1)/2
       :)
      let $s1 :=
        (1.0 + point:px($a) * point:px($a) + point:py($a) * point:py($a)) div 2.0
      let $s2 :=
        (1.0 + point:px($b) * point:px($b) + point:py($b) * point:py($b)) div 2.0
      let $center :=
        point:point(
          ($s1 * point:py($b) - $s2 * point:py($a)) div $axby-bxay,
          (point:px($a) * $s2 - point:px($b) * $s1) div $axby-bxay
        )
      (: √(x² + y² - 1) why -1? :)
      let $r :=
        math:sqrt(
          point:px($center)*point:px($center) +
          point:py($center)*point:py($center) - 1.0
        )
      (:
         r^2 / ((xp - xc)² + (yp - yc)²)
       :)
      let $factor :=
        $r * $r div
        (
          (point:px($pt)-point:px($center)) * (point:px($pt)-point:px($center)) +
          (point:py($pt)-point:py($center)) * (point:py($pt)-point:py($center))
        )
      return
        (: shift scale shift: $center + $factor * ($pt - $center) :)
        point:point(
          point:px($center) + $factor * (point:px($pt) - point:px($center)),
          point:py($center) + $factor * (point:py($pt) - point:py($center))
        )
    )
  )
}

declare function this:project(
  $scale as xs:numeric,
  $points as map(xs:string,item()*)*
) as map(xs:string,item()*)*
{
  for $pt in $points
  let $x := point:px($pt) div $scale
  let $y := point:py($pt) div $scale
  let $div := $x*$x + $y+$y
  return
    if ($div = 0) then (
      $pt
    ) else (
      geom:scale(
        point:point($x div $div, $y div $div),
        $scale
      )
    )
}

Original Source Code

xquery version "3.1";
(:~
 : Module with functions providing some handy hyperbolic geometric operations.
 : All operations are performed on the unit circle centered on origin, then
 : scaled by the hyperbolic scale.
 :
 : Copyright© Mary Holstege 2020-2023
 : CC-BY (https://creativecommons.org/licenses/by/4.0/)
 : @since March 2023
 : @custom:Status Stable
 :)
module namespace this="http://mathling.com/geometric/hyperbolic"; 

import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy";
import module namespace config="http://mathling.com/core/config"
       at "../core/config.xqy";
import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy";
import module namespace geom="http://mathling.com/geometric"
       at "../geo/euclidean.xqy";
import module namespace point="http://mathling.com/geometric/point"
       at "../geo/point.xqy";
import module namespace ellipse="http://mathling.com/geometric/ellipse"
       at "../geo/ellipse.xqy";
import module namespace box="http://mathling.com/geometric/rectangle"
       at "../geo/rectangle.xqy";
import module namespace edge="http://mathling.com/geometric/edge"
       at "../geo/edge.xqy";
import module namespace path="http://mathling.com/geometric/path"
       at "../geo/path.xqy";

declare namespace map="http://www.w3.org/2005/xpath-functions/map";
declare namespace array="http://www.w3.org/2005/xpath-functions/array";
declare namespace math="http://www.w3.org/2005/xpath-functions/math";

declare %private variable $this:TOLERANCE as xs:double := 0.00001;


(:~
 : circle-through()
 : Compute the circle containing points a and b, and
 : orthogonal to the unit circle (a.k.a shortest path between
 : a and b in the hyperbolic sense).
 :
 : Returns empty if the shortest path is a straight line (diameter
 : of the unit circle).
 :)
declare function this:circle-through(
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)?
{
  let $ax := point:px($a)
  let $ay := point:py($a)
  let $bx := point:px($b)
  let $by := point:py($b)
  let $axby-aybx := $ax * $by - $ay * $bx
  where abs($axby-aybx) >= $this:TOLERANCE
  return (
    let $cx :=
      -0.5 *
        ($ay * $bx * $bx + $ay * $by * $by -
          ($ax * $ax + $ay * $ay + 1) * $by + $ay
        )
        div $axby-aybx
    let $cy :=
      0.5 *
        ($ax * $bx * $bx + $ax * $by * $by -
          ($ax * $ax + $ay * $ay + 1) * $bx + $ax
        )
        div $axby-aybx
    let $radius :=
      abs(
        -0.5 * math:sqrt(
          $ax * $ax * $ax * $ax * $bx * $bx +     
          $ax * $ax * $ax * $ax * $by * $by -     
          2 * $ax * $ax * $ax * $bx * $bx * $bx - 
          2 * $ax * $ax * $ax * $bx * $by * $by + 
          2 * $ax * $ax * $ay * $ay * $bx * $bx + 
          2 * $ax * $ax * $ay * $ay * $by * $by - 
          2 * $ax * $ax * $ay * $bx * $bx * $by - 
          2 * $ax * $ax * $ay * $by * $by * $by + 
          $ax * $ax * $bx * $bx * $bx * $bx +     
          2 * $ax * $ax * $bx * $bx * $by * $by + 
          $ax * $ax * $by * $by * $by * $by -     
          2 * $ax * $ay * $ay * $bx * $bx * $bx - 
          2 * $ax * $ay * $ay * $bx * $by * $by +
          $ay * $ay * $ay * $ay * $bx * $bx +
          $ay * $ay * $ay * $ay * $by * $by -
          2 * $ay * $ay * $ay * $bx * $bx * $by -
          2 * $ay * $ay * $ay * $by * $by * $by +
          $ay * $ay * $bx * $bx * $bx * $bx +
          2 * $ay * $ay * $bx * $bx * $by * $by +
          $ay * $ay * $by * $by * $by * $by -
          2 * $ax * $ax * $ax * $bx -
          2 * $ax * $ax * $ay * $by +
          4 * $ax * $ax * $bx * $bx -
          2 * $ax * $ay * $ay * $bx +
          8 * $ax * $ay * $bx * $by -
          2 * $ax * $bx * $bx * $bx -
          2 * $ax * $bx * $by * $by -
          2 * $ay * $ay * $ay * $by +
          4 * $ay * $ay * $by * $by -
          2 * $ay * $bx * $bx * $by -
          2 * $ay * $by * $by * $by +
          $ax * $ax - 2 * $ax * $bx +
          $ay * $ay - 2 * $ay * $by +
          $bx * $bx + $by * $by
        )
        div $axby-aybx
      )
    return ellipse:circle(point:point($cx, $cy), $radius)
  )
};

(:~
 : arc-through()
 : Compute the hyperbolic arc through the two points.
 : If the points are near-aligned it returns a straight edge, otherwise the
 : circular arc between the points.
 :
 : @param $a: one point
 : @param $b: another point
 : @return the hyperbolic arc from a to b
 :)
declare function this:arc-through(
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $ax := point:px($a)
  let $ay := point:py($a)
  let $bx := point:px($b)
  let $by := point:py($b)
  let $circle := this:circle-through($a, $b)
  return (
    if (empty($circle)) then ( (: near-aligned: straight line :)
      edge:edge($a, $b)
    ) else (
      let $c := ellipse:center($circle)
      let $r := ellipse:radius($circle)
      let $θ := -math:atan2(-$ay + point:py($c), $ax - point:px($c))
      let $θ2 := -math:atan2(-$by + point:py($c), $bx - point:px($c))
      (: Recenter the angles to [0, 2π] :)
      let $θ := util:remap-radians($θ)      
      let $θ2 := util:remap-radians($θ2)      
      let $delta := abs($θ - $θ2)
      return (
        if (($θ > $θ2 and $delta > math:pi()) or ($θ < $θ2 and $delta < math:pi()))
        then edge:arc($c, $r, $a, $b, false())
        else edge:arc($c, $r, $a, $b, true())
      )
    )
  )
};

(:~
 : edge()
 : Compute the scaled hyperbolic edge between a and b.
 : @param $scale: hyperbolic scale
 : @param $a: one point
 : @param $b: another point
 : @return the scaled hyperbolic arc from a to b
 :)
declare function this:edge(
  $scale as xs:numeric,
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:arc-through($a, $b)=>geom:scale($scale)
};

(:~
 : line()
 : Construct a hyperbolic straight line (edge), which is to say, the arc
 : between the points.
 :
 : @param: $scale: hyperbolic scale
 : @param: $start: start of edge
 : @param: $end: end of edge
 : @return scaled hyperbolic line segment (an arc)
 :)
declare function this:line(
  $scale as xs:numeric,
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:edge($scale, $start, $end)
};

(:~
 : polygon()
 : Construct a scaled hyperbolic polygon from the points.
 :
 : @param: $scale: hyperbolic scale
 : @param: $points: set of points
 : @return scaled hyperbolic polygon (arc edges)
 :)
declare function this:shape(
  $scale as xs:numeric,
  $points as map(xs:string,item()*)*
) as map(xs:string,item()*)
{
  this:polygon($scale, $points)
};

(:~
 : box()
 : Construct a scaled hyperbolic box from the axis-aligned box.
 :
 : @param: $scale: hyperbolic scale
 : @param $box: axis-aligned box the extreme points
 : @return a scaled hyperbolic polygon from the box's points
 :)
declare function this:box(
  $scale as xs:numeric,
  $box as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $points := (
    point:point(box:min-px($box),box:min-py($box)),
    point:point(box:min-px($box),box:max-py($box)),
    point:point(box:max-px($box),box:max-py($box)),
    point:point(box:max-px($box),box:min-py($box)),
    point:point(box:min-px($box),box:min-py($box))
  )
  return this:polygon($scale, $points)
};

(:~
 : triangle()
 : Construct a scaled hyperbolic triangle with vertices (a,b,c).
 : 
 : @param: $scale: hyperbolic scale
 : @param $a: one point
 : @param $b: second point
 : @param $c: third point
 : @return a scaled hyperbolic polygon from the three points
 :)
declare function this:triangle(
  $scale as xs:numeric,
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*),
  $c as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:polygon( $scale, ($a, $b, $c) )
};

(:~ 
 : supporting-line()
 : Supporting lines (circles) through a and b; for debugging.
 :
 : @param $scale: hyperbolic scaling parameter
 : @param $a: one point
 : @param $b: the other point
 : @return either a circle or a straight edge
 :)
declare function this:supporting-line(
  $scale as xs:numeric,
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $ax := point:px($a)
  let $ay := point:py($a)
  let $bx := point:px($b)
  let $by := point:py($b)
  let $circle := this:circle-through($a, $b) (: Compute circle parameters :)
  return (
    if (empty($circle)) then (
      (: Project a and b points onto the circle :)
      let $θ := math:atan2($ay, $ax)
      let $θ2 := math:atan2($by, $bx)
      return (
        edge:edge(
          point:point($scale*math:cos($θ), $scale*math:sin($θ)),
          point:point($scale*math:cos($θ2), $scale*math:sin($θ2))
        )
      )
    ) else (
      ellipse:circle($point:ORIGIN, $scale*ellipse:radius($circle))
    )
  )
};

(:~
 : all-support-lines()
 : Get all the supporting lines from one point to the next; for debugging.
 :
 : @param $scale: hyperbolic scale
 : @param $pts: sequence of points
 : @return sequence of supporting hyperbolic edges
 :)
declare function this:all-support-lines(
  $scale as xs:numeric,
  $pts as map(xs:string,item()*)*
) as map(xs:string,item()*)*
{
  for $p at $i in $pts
  let $next := 
    if (empty($pts[$i + 1]))
    then $pts[1]
    else $pts[$i + 1]
  return
    this:supporting-line($scale, $p, $next)
};

(:~
 : polygon()
 : Construct a hyperbolic polygon with the provided vertices.
 : In general you want points in the unit circle, then scaled by scale
 : to size of your circle.
 :
 : @param: $scale: hyperbolic scale
 : @param: $pts: sequence of vertices
 : @return scaled hyperbolic polygon (arc edges)
 :)
declare function this:polygon(
  $scale as xs:numeric,
  $pts as map(xs:string,item()*)*
) as map(xs:string,item()*)
{
  let $edges := 
    for $p at $i in $pts
    let $next := 
      if (empty($pts[$i + 1]))
      then $pts[1]
      else $pts[$i + 1]
    where not(this:within-tolerance($p=>geom:scale($scale), $next=>geom:scale($scale), $this:TOLERANCE))
    return
      this:edge($scale, $p, $next)
  return (
    path:polygon($edges)
  )
};

(:~
 : distance()
 : Scaled hyperbolic distance between two points.
 :
 : @param $scale: hyperbolic scale
 : @param $p: one point
 : @param $q: the other point
 : @return hyperbolic distance between p and q, scaled by scale
 :)
declare function this:distance(
  $scale as xs:numeric,
  $p as map(xs:string,item()*),
  $q as map(xs:string,item()*)
) as xs:double
{
  let $pqd := geom:distance($p,$q)
  let $pd := geom:distance(point:point(0,0),$p)
  let $qd := geom:distance(point:point(0,0),$q)
  let $x :=  1 + 2*($pqd*$pqd)*$scale*$scale div (($scale*$scale - $pd*$pd)*($scale*$scale - $qd*$qd))
  (: acosh(x) = ln(x + sqrt(x^2 - 1)) :)
  return (
    math:log($x + math:sqrt($x*$x - 1))
  )
};

(: === tiling === :)

(:~
 : center-polygon()
 : Center polygon in n-k regular tiling.
 :
 : @param $n number of edges of tiles
 : @param $k number of tiles adjacent at each vertex
 : @return central polygon of a regular n-k tiling with straight edges
 :)
declare %private function this:center-polygon(
  $n as xs:integer,
  $k as xs:integer
) as map(xs:string,item()*)
{
  (: 
   : ABC is triangle in regular n,k tiling where
   :   A is center (also center of circle)
   :   B is a vertex of n-gon
   :   C is midpoint of side of n-gon adjacent to B
   :)

   let $angleA := math:pi() div $n
   let $angleB := math:pi() div $k
   let $angleC := math:pi() div 2
 
   (:
    : Compute distance s from A to B to make regular tiling
    :)
 
   let $sinA := math:sin($angleA)
   let $sinB := math:sin($angleB)
   let $s :=
     math:sin($angleC - $angleB - $angleA) div
     math:sqrt(1 - $sinB * $sinB - $sinA * $sinA)

   (: 
    : Determine the coordinates of the n vertices of the n-gon.
    : They're all at distance s from center of the Poincare disk
    :)
   let $points :=
     for $i in 1 to $n
     return (
       point:point(
         $s * math:cos($angleA * (3 + 2*$i)),
         $s * math:sin($angleA * (3 + 2*$i))
       )
     )
   return path:polygon(geom:to-edges($points))
};

(:~
 : replicate-polygon()
 : Create new polygon in tesselation from seed polygon.
 :
 : Rule codes
 :   0:  initial polygon.  Needs neighbors on all n sides
 :   1:  polygon already has 2 neighbors, but one less
 :       around corner needed
 :   2:  polygon already has 2 neighbors
 :   3:  polygon already has 3 neighbors
 :   4:  polygon already has 4 neighbors
 :)
declare %private function this:replicate-polygon(
  $polygon as map(xs:string,item()*),
  $rule as xs:integer,
  $n as xs:integer,
  $k as xs:integer,
  $layer as xs:integer
) as map(xs:string,item()*)*
{
  if ($layer = 0) then ()
  else (
    let $special := ($rule = 1)
    let $r := if ($special) then 2 else $rule
    let $start := if ($r = 4) then 3 else 2
    let $quantity := if ($k = 3 and $rule != 0) then $n - $r - 1 else $n - $r
    for $s in $start to $start + $quantity
    return (
      let $newp := this:reflect-polygon($polygon, 1 + ($s mod $n), $n, $k)
      let $newr := if ($k = 3 and $s = $start and $r != 0) then 4 else 3
      return (
        $newp,
        this:replicate-polygon($newp, $newr, $n, $k, $layer - 1)
      )
    )
  )
};

(:~
 : within-tolerance()
 : Points are two close to count as usefully distinct for polygon construction.
 :
 : @param $p: one point
 : @param $q: one point
 : @param $tolerance: how close is too close
 : @return true if x and y coordinates are too close
 :)
declare %private function this:within-tolerance(
  $p as map(xs:string,item()*),
  $q as map(xs:string,item()*),
  $tolerance as xs:numeric
) as xs:boolean
{
  (abs(point:px($p) - point:px($q)) le $tolerance) and
  (abs(point:py($p) - point:py($q)) le $tolerance)
};

(:~
 : same-polygon()
 : The two polygons are the same (or near enough). For pruning.
 :
 : @param $poly: one polygon
 : @param $other: the other polygon
 : @return true polygons are pretty much the same
 :)
declare %private function this:same-polygon(
  $poly as map(xs:string,item()*),
  $other as map(xs:string,item()*)
) as xs:boolean
{  
  let $pv := geom:vertices($other)
  let $ov := geom:vertices($other)
  return (
    (
      every $v in $pv satisfies (
        geom:shortest-distance($v, $other) <= $this:TOLERANCE
      )
    )
    and
    (
      every $v in $ov satisfies (
        geom:shortest-distance($v, $poly) <= $this:TOLERANCE
      )
    )
  )
};

(:~
 : remove-duplicates()
 : Prune duplicate polygons.
 :
 : @param $scale: hyperbolic scale
 : @param $polygons: sequence of polygons to prune
 : @return sequence of polygons with duplicates removed
 :)
declare %private function this:remove-duplicates(
  $scale as xs:numeric,
  $polygons as map(xs:string,item()*)*
) as map(xs:string,item()*)*
{
  if (count($polygons)=1) then $polygons else
  let $scaled :=
    for $poly in $polygons return $poly=>geom:scale($scale)=>geom:snap()
  for $poly at $i in $scaled
  where every $other in $scaled[position() < $i]
        satisfies not(this:same-polygon($poly, $other))
  return $polygons[$i]
};

(:~
 : tesselate()
 : Construct an n-k hyperbolic tesselation. For a tesselation to be possible
 : (n-2)(k-2) <= 4. The result is an set of unscaled straight-edged polygons:
 : you'll need to call sorted-hypers to convert them to properly scaled
 : hyperbolic polygons and then translate them to the appropriate position.
 :
 : @param $scale: hyperbolic scale (radius of circle)
 : @param $n: number of vertices in tiles
 : @param $k: number of polygons adjacent at each vertex
 : @param $layer: number of layers or rings of polygons to add
 : @return sequence of unscaled polygons
 : @error GEOM-INVHYPER if no such tesselation is possible 
 :)
declare function this:tesselate(
  $scale as xs:numeric,
  $n as xs:integer,
  $k as xs:integer,
  $layer as xs:integer
) as map(xs:string,item()*)*
{
  if (($n - 2)*($k - 2) <= 4)
  then errors:error("GEOM-INVHYPER", ($n - 2) * ($k - 2))
  else ()
  ,
  this:remove-duplicates(
    $scale,
    let $center := this:center-polygon($n, $k)
    let $rule := 0
    return (
      $center,
      this:replicate-polygon($center, $rule, $n, $k, $layer - 1)
    )
  )
};

(:~ 
 : sorted-hypers()
 : Sort the polygons produced from tessalate() and turn them into
 : scaled hyperbolic polygons. The hyperbolic polygons will be ordered
 : by distance from the center and angle, i.e. a spiral out from the center.
 : They will have the additional properties "layer" (based on distance) to
 : indicate which generation they were and "angle" indicating their angle
 : from the center. These can be used to apply colouring schemes.
 :
 : @param $size: scaling of the polygons
 : @param $polys: sequence of unscaled polygons
 : @return sequence of scaled hyperbolic polygons
 :)
declare function this:sorted-hypers(
  $size as xs:numeric,
  $polys as map(xs:string,item()*)*
) as map(xs:string,item()*)*
{
  (: Sorted and turned into hyperpolygons :)
  let $origin := point:point(0,0)
  for $poly in $polys
  let $bounds := geom:bounding-box($poly=>geom:scale($size))
  let $center := geom:region-center($bounds)
  let $distance := round(this:distance($size, $origin, $center))
  let $angle := round(geom:path-angle($origin, $center))
  order by $distance ascending, $angle ascending
  return (
    this:polygon($size, geom:vertices($poly))=>
      map:put("layer", $distance)=>
      map:put("angle", $angle)
  )
};

(:~
 : reflect-polygon()
 : Reflect the polygon through the point or the side indicated by the side s
 : to produce the resulting polygon
 :)
declare %private function this:reflect-polygon(
  $polygon as map(xs:string,item()*),
  $s as xs:integer,
  $n as xs:integer,
  $k as xs:integer
) as map(xs:string,item()*)?
{
  let $vertices := geom:vertices($polygon)
  let $start := $vertices[$s]
  let $end := $vertices[$s + 1] (: polygon vertices include closing point :)
  let $points := (
    for $i in 1 to $n
    let $j := 1 + ($n + $s - $i + 1) mod $n
    order by $j ascending
    return this:reflect($vertices[$i], $start, $end)
  )
  return (
    if (count($points) le 2)
    then util:log("too few points="||geom:quote($points))
    else path:polygon(geom:to-edges($points))
  )
};

(:~
 : reflect()
 : Reflect the point across the given line. Within hyperbolic unit circle.
 : @param $pt: the point
 : @param: $a: start of edge
 : @param: $b: end of edge
 : @return reflected point
 :)
declare function this:reflect(
  $pt as map(xs:string,item()*),
  $a as map(xs:string,item()*),
  $b as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  (: Find a unit vector $uv in the direction of the line :)
  (: $axby-bxay = ax*by - bx*ay :)
  let $axby-bxay := point:px($a) * point:py($b) - point:px($b) * point:py($a)
  return (
    if (false()) then (
      (: Euclidean reflection = 
         Ref(l,v) = 2 (v·l)/(l·l) - v = 2 Proj(l,v) - v
       :)
      (: distance between a and b √((ax - bx)² + (ax - by)²) :)
      let $delta :=
        math:sqrt(
          (point:px($a) - point:px($b))*(point:px($a) - point:px($b)) +
          (point:py($a) - point:py($b))*(point:py($a) - point:py($b))
        )
      let $uv :=
        point:point(
          (point:px($b) - point:px($a)) div $delta,
          (point:py($b) - point:py($a)) div $delta
        )
      (: 2 * ((px - ax)*ux + (py-ay)*uy) :)
      let $factor :=
        2.0 * (
          (point:px($pt) - point:px($a)) * point:px($uv) +
          (point:py($pt) - point:py($a)) * point:py($uv)
        )
      return
        (: 2 * $a + $factor * $uv - $p :)
        point:point(
          2.0 * point:px($a) + $factor * point:px($uv) - point:px($pt),
          2.0 * point:py($a) + $factor * point:py($uv) - point:py($pt)
        )
    ) else if (abs($axby-bxay) < $this:TOLERANCE) then ( (: straight :)
      geom:reflect($pt, $a, $b)
    ) else (
      (: Find the center of the circle through these points :)
      (:  
       :  s1 = (xa² + ya² + 1)/2
       :  s2 = (xb² + yb² + 1)/2
       :)
      let $s1 :=
        (1.0 + point:px($a) * point:px($a) + point:py($a) * point:py($a)) div 2.0
      let $s2 :=
        (1.0 + point:px($b) * point:px($b) + point:py($b) * point:py($b)) div 2.0
      let $center :=
        point:point(
          ($s1 * point:py($b) - $s2 * point:py($a)) div $axby-bxay,
          (point:px($a) * $s2 - point:px($b) * $s1) div $axby-bxay
        )
      (: √(x² + y² - 1) why -1? :)
      let $r :=
        math:sqrt(
          point:px($center)*point:px($center) +
          point:py($center)*point:py($center) - 1.0
        )
      (:
         r^2 / ((xp - xc)² + (yp - yc)²)
       :)
      let $factor :=
        $r * $r div
        (
          (point:px($pt)-point:px($center)) * (point:px($pt)-point:px($center)) +
          (point:py($pt)-point:py($center)) * (point:py($pt)-point:py($center))
        )
      return
        (: shift scale shift: $center + $factor * ($pt - $center) :)
        point:point(
          point:px($center) + $factor * (point:px($pt) - point:px($center)),
          point:py($center) + $factor * (point:py($pt) - point:py($center))
        )
    )
  )
};

(:~
 : project()
 : Inversion in unit circle.
 : Φ(x,y) = (u,v) = (x/(x² + y²),y/(x² + y²))
 :
 : http://homepages.gac.edu/~hvidsten/geom-text/web-chapters/hyper-transf.pdf:
 : "all other hyperbolic reflections about lines that are not diameters can be
 : expressed as inversion through the circle which defines the line"
 :)
declare function this:project(
  $scale as xs:numeric,
  $points as map(xs:string,item()*)*
) as map(xs:string,item()*)*
{
  for $pt in $points
  let $x := point:px($pt) div $scale
  let $y := point:py($pt) div $scale
  let $div := $x*$x + $y+$y
  return
    if ($div = 0) then (
      $pt
    ) else (
      geom:scale(
        point:point($x div $div, $y div $div),
        $scale
      )
    )
};