
   0 | module Geom.Angle
   1 | 
   2 | import Chem.Util
   3 | import Derive.Prelude
   4 | 
   5 | %language ElabReflection
   6 | %default total
   7 | 
   8 | public export
   9 | Epsilon : Double
  10 | Epsilon = 1.0e-12
  11 | 
  12 | --------------------------------------------------------------------------------
  13 | --          Angle
  14 | --------------------------------------------------------------------------------
  15 | 
  16 | export
  17 | TwoPi : Double
  18 | TwoPi = 2 * pi
  19 | 
  20 | ||| Convert a floating point number to an angle
  21 | ||| in the interval [0, 2 * pi).
  22 | export
  23 | normalize : Double -> Double
  24 | normalize x =
  25 |   if x < 0 then
  26 |     let f := ceiling (abs x / TwoPi) in x + f * TwoPi
  27 |   else if x >= TwoPi then
  28 |     let f := floor (x / TwoPi) in x - f * TwoPi
  29 |   else x
  30 | 
  31 | ||| A normalized angle in the range [0,2 * pi).
  32 | public export
  33 | record Angle where
  34 |   constructor A
  35 |   ||| The normalized value
  36 |   value : Double
  37 | 
  38 |   ||| The original floating point number from which the angle
  39 |   ||| was created.
  40 |   0 org : Double
  41 | 
  42 |   ||| Proof that we did not forget the normalization step.
  43 |   {auto 0 prf : value === normalize org}
  44 | 
  45 | %runElab derive "Angle" [Show,Eq,Ord]
  46 | 
  47 | ||| Convenience constructor for angles.
  48 | |||
  49 | ||| It is safe to invoke `A` directly, but one has to deal
  50 | ||| with the proofing stuff.
  51 | export %inline
  52 | angle : Double -> Angle
  53 | angle v = A (normalize v) v
  54 | 
  55 | export
  56 | halfPi : Angle
  57 | halfPi = angle (pi / 2)
  58 | 
  59 | export
  60 | twoThirdPi : Angle
  61 | twoThirdPi = angle (TwoPi / 3)
  62 | 
  63 | export
  64 | threeHalfPi : Angle
  65 | threeHalfPi = angle (3 * pi / 2)
  66 | 
  67 | export
  68 | pi : Angle
  69 | pi = angle pi
  70 | 
  71 | export
  72 | zero : Angle
  73 | zero = angle 0
  74 | 
  75 | export
  76 | fullSteps : Nat -> Angle
  77 | fullSteps 0 = zero
  78 | fullSteps n = angle (TwoPi / cast n)
  79 | 
  80 | ||| Returns the absolute distance between two angles
  81 | |||
  82 | ||| Unlike `minDelta`, this just subtracts the smaller angle
  83 | ||| from the larger one. Therefore, `delta zero halfPi = halfPi`
  84 | ||| and `delta zero threeHalfPi = threeHalfPi`.
  85 | export
  86 | delta : Angle -> Angle -> Angle
  87 | delta (A x _) (A y _) = angle (abs $ x - y)
  88 | 
  89 | ||| Addition of two angles
  90 | export %inline
  91 | (+) : Angle -> Angle -> Angle
  92 | A x _ + A y _ = angle $ x + y
  93 | 
  94 | ||| Difference between two angles
  95 | export %inline
  96 | (-) : Angle -> Angle -> Angle
  97 | A x _ - A y _ = angle $ x - y
  98 | 
  99 | ||| The inverse of an angle so that `x + negate x == 0`
 100 | ||| (modulo rounding errors).
 101 | export %inline
 102 | negate : Angle -> Angle
 103 | negate (A x _) = angle $ TwoPi - x
 104 | 
 105 | ||| Multiplication of an angle with a constant factor
 106 | export %inline
 107 | (*) : Double -> Angle -> Angle
 108 | v * A x _ = angle $ v * x
 109 | 
 110 | ||| Divides an angle into `n` equal part.
 111 | |||
 112 | ||| Returns the angle unmodified in case `n` equals zero.
 113 | export
 114 | divide : Nat -> Angle -> Angle
 115 | divide 0 a       = a
 116 | divide n (A x _) = angle (x / cast n)
 117 | 
 118 | ||| Returns the shortest distance between two angles.
 119 | ||| (either clockwise or counterclockwise.)
 120 | |||
 121 | ||| Unlike `delta`, this computes the *smaller* angle
 122 | ||| between the two input values.
 123 | ||| Therefore, `minDelta zero halfPi = minDelta zero threeHalfPi = halfPi`.
 124 | export
 125 | minDelta : Angle -> Angle -> Angle
 126 | minDelta x y = min (delta x y) (negate (delta x y))
 127 | 
 128 | ||| Convert and angle to centigrees
 129 | export %inline
 130 | toDegree : Angle -> Double
 131 | toDegree a = a.value * 180 / pi
 132 | 
 133 | ||| Convert an angle in centigrees to one in radians
 134 | export %inline
 135 | fromDegree : Double -> Angle
 136 | fromDegree = angle .  (*) (pi/180)
 137 | 
 138 | ||| Angle bisector counterclockwise
 139 | export
 140 | bisector : Angle -> Angle -> Angle
 141 | bisector x y = x + 0.5 * (y - x)
 142 | 
 143 | ||| From a list of angles, returns the one closest to the given angle
 144 | export %inline
 145 | closestAngle : Angle -> List Angle -> Maybe Angle
 146 | closestAngle = minBy . delta
 147 | 
 148 | shiftZip : List a -> List (a,a)
 149 | shiftZip []        = []
 150 | shiftZip (x :: xs) = zip (x::xs) (xs++[x])
 151 | 
 152 | ||| Given an angle `phi` plus a list of angles, returns from the
 153 | ||| list two angles `a` and `b`, so that `phi` lies between `a` and `b`
 154 | ||| with no other angle closer to `phi`.
 155 | |||
 156 | ||| Note: In case of this being successful, `phi` will always lie counter
 157 | |||       clockwise of the first returned angle and the second returned
 158 | |||       angle will lie counter-clockwise of `phi`.
 159 | export
 160 | enclosingAngles : Angle -> List Angle -> Maybe (Angle,Angle)
 161 | enclosingAngles x xs =
 162 |   case sort xs of
 163 |     []    => Nothing
 164 |     a::as => case find (\(y,z) => y <= x && x <= z) (zip (a::as) as) of
 165 |       Nothing => Just (last $ a::as,a)
 166 |       Just p  => Just p
 167 | 
 168 | export
 169 | largestBisector : List Angle -> Angle
 170 | largestBisector xs =
 171 |   fromMaybe zero . map (uncurry bisector) . maxBy diff . pairs $ sort xs
 172 |   where
 173 |     pairs : List a -> List (a,a)
 174 |     pairs [] = []
 175 |     pairs (h::t) = zip (h::t) (t ++ [h])
 176 | 
 177 |     diff : (Angle, Angle) -> Angle
 178 |     diff (x,y) = y - x
 179 | 
 180 |     maxBy : Ord b => (a -> b) -> List a -> Maybe a
 181 |     maxBy f []     = Nothing
 182 |     maxBy f (h::t) = Just $ foldl (\x,y => if f x >= f y then x else y) h t
 183 | 
 184 | --------------------------------------------------------------------------------
 185 | --          Regular n-gons and regular arcs
 186 | --------------------------------------------------------------------------------
 187 | 
 188 | ||| An arc of `segs` segments connecting two existing nodes.
 189 | |||
 190 | ||| Used for creating rings and bridges when placing molecules.
 191 | public export
 192 | record Arc where
 193 |   constructor MkArc
 194 |   ||| Total angle of the arc
 195 |   angle    : Angle
 196 | 
 197 |   ||| Angle of a single segment of the arc
 198 |   step     : Angle
 199 | 
 200 |   ||| Radius of the arc
 201 |   radius   : Double
 202 | 
 203 |   ||| Distance between the center of the arc and the middle
 204 |   ||| of a segment.
 205 |   distance : Double
 206 | 
 207 | %runElab derive "Arc" [Show,Eq]
 208 | 
 209 | ||| Angle at the corner of a regular n-gon
 210 | export %inline
 211 | ngonAngle : (n : Nat) -> (0 prf : LTE 3 n) => Angle
 212 | ngonAngle n = pi - fullSteps n
 213 | 
 214 | ||| Radius of a regular n-gon with side length `side`.
 215 | |||
 216 | ||| This is the distance from the center of the n-gon to one of its
 217 | ||| corners.
 218 | export %inline
 219 | ngonRadius : (side : Double) -> (n : Nat) -> (0 prf : LTE 3 n) => Double
 220 | ngonRadius side n = 0.5 * side / cos ((ngonAngle n).value / 2)
 221 | 
 222 | segDist : Double -> Double -> Double
 223 | segDist side rad = sqrt (pow rad 2 - pow (side / 2) 2)
 224 | 
 225 | ||| Distance from the center of a regular n-gon with side length `side`
 226 | ||| to the middle of one of its sides.
 227 | export %inline
 228 | ngonDistance : (side : Double) -> (n : Nat) -> (0 prf : LTE 3 n) => Double
 229 | ngonDistance side n = segDist side (ngonRadius side n)
 230 | 
 231 | -- radius of an arc of angle `phi` connecting two points
 232 | -- with a distance of `d`
 233 | arcRadius : (d,phi : Double) -> Double
 234 | arcRadius d phi = d / (2 * sin (phi / 2))
 235 | 
 236 | -- length of a single line segment, out of `n`
 237 | -- segments approximating an arc of angle `phi` and connecting two
 238 | -- points `x` and `y`, with a distance of `d` between `x` and `y`.
 239 | segmentLength : Nat -> (d,phi : Double) -> Double
 240 | segmentLength n d phi = 2 * (arcRadius d phi) * sin (phi / (2 * cast (S n)))
 241 | 
 242 | mkArc : Nat -> (d,phi : Double) -> Arc
 243 | mkArc n d phi =
 244 |  let r := arcRadius d phi
 245 |      a := angle phi
 246 |   in MkArc a (divide (S n) a) r (segDist d r)
 247 | 
 248 | export
 249 | ngon : (side : Double) -> (n : Nat) -> (0 prf : LTE 3 n) => Arc
 250 | ngon side n =
 251 |  let a := fullSteps n
 252 |   in MkArc (negate a) a (ngonRadius side n) (ngonDistance side n)
 253 | 
 254 | ||| Computes the dimensions of an arc that connects two
 255 | ||| existing nodes (with a distance of `d` between the nodes)
 256 | ||| by inserting `n` additional nodes between them.
 257 | |||
 258 | ||| The arc is optimized in such a way that the distance between two
 259 | ||| adjacent nodes is as close to `len` as possible.
 260 | |||
 261 | ||| Note: If `len` is too short, that is `(n+1) * len < d`, this returns
 262 | |||       close to - but not exactly - a straight line. Client code is
 263 | |||       responsible to choose `len` in such a manner, that an arc
 264 | |||       with a reasonable minimal curvature is constructed,
 265 | |||       for instance by setting the minimal `len`
 266 | |||       at `(1+delta)*d / (n+1)` with `delta > 0`.
 267 | export
 268 | arc : (n : Nat) -> (0 prf : IsSucc n) => (len,d : Double) -> Arc
 269 | arc n@(S k) len d =
 270 |   case abs (len-d) < Epsilon of
 271 |     -- this is just (or close to) a regular n-gon
 272 |     True  => ngon len (S $ S $ S k)
 273 |     -- run a binary search to find the ideal arc
 274 |     False => find 64 Epsilon (2*pi - Epsilon)
 275 |   where
 276 |     best : (l,u : Double) -> Arc
 277 |     best l u =
 278 |      let ll := segmentLength n d l
 279 |          lu := segmentLength n d u
 280 |       in if abs (lu-len) <= abs (ll-len) then mkArc n d u else mkArc n d l
 281 | 
 282 |     -- a binary search of at most `iter` iterations to find the
 283 |     -- ideal arc angle `phi` (with current lower and upper bounds `l` and `u`)
 284 |     find : (iter : Nat) -> (l,u : Double) -> Arc
 285 |     find 0     l u = best l u
 286 |     find (S k) l u =
 287 |       case abs (u-l) <= Epsilon of
 288 |         True  => best l u
 289 |         False => case segmentLength n d ((l+u)/2.0) >= len of
 290 |           True  => find k l ((l+u)/2.0)
 291 |           False => find k ((l+u)/2.0) u
 292 | 
 293 | --------------------------------------------------------------------------------
 294 | --          Tests and proofs
 295 | --------------------------------------------------------------------------------
 296 | 
 297 | 0 PI : Double
 298 | PI = pi
 299 | 
 300 | 0 DoubleEq : Double -> Double -> Type
 301 | DoubleEq x y = (abs (x - y) < 0.0000000001) === True
 302 | 
 303 | 0 Norm : Double -> Double
 304 | Norm x = normalize x
 305 | 
 306 | 0 normalize0 : DoubleEq (Norm 0) 0.0
 307 | normalize0 = Refl
 308 | 
 309 | 0 normalize2pi : DoubleEq (Norm $ 2 * PI) 0
 310 | normalize2pi = Refl
 311 | 
 312 | 0 normalizeAny : DoubleEq (Norm $ 21.54 * PI) (1.54 * PI)
 313 | normalizeAny = Refl
 314 | 
 315 | 0 normalizeNeg : DoubleEq (Norm $ negate PI) PI
 316 | normalizeNeg = Refl
 317 | 
 318 | 0 normalizeAnyNeg : DoubleEq (Norm $ (-21.54) * PI) ((2 - 1.54) * PI)
 319 | normalizeAnyNeg = Refl
 320 |