2 | import Derive.Prelude
10 | %language ElabReflection
32 | record AffineTransformation where
34 | transform : LinearTransformation
35 | translate : Vector Id
39 | Id : AffineTransformation
44 | scaling : Scale -> AffineTransformation
45 | scaling sc = AT (scaling sc) vzero
48 | Semigroup AffineTransformation where
49 | AT tl vl <+> AT tr vr =
50 | let vr2 := transform tl vr in AT (tl <+> tr) (vr2 + vl)
53 | Monoid AffineTransformation where neutral = Id
58 | inverse : AffineTransformation -> AffineTransformation
59 | inverse (AT tr v) = let i := inverse tr in AT i $
negate (transform i v)
66 | translate : Vector Id -> AffineTransformation
67 | translate v = AT Id v
85 | record Point (t : AffineTransformation) where
90 | %runElab deriveIndexed "Point" [Show,Eq,Ord]
100 | (-) : Point t -> Point t -> Vector t.transform
101 | (-) (P x1 y1) (P x2 y2) = V (x1 - x2) (y1 - y2)
106 | convert : {s,t : _} -> Point s -> Point t
108 | let AT l (V dx dy) := t <+> inverse s
109 | V x2 y2 := transform l (p - origin)
110 | in P (x2 + dx) (y2 + dy)
113 | {s,t : _} -> Cast (Point s) (Point t) where cast = convert
122 | interface GetPoint (0 a : Type) where
123 | gtrans : AffineTransformation
124 | point : a -> Point gtrans
127 | 0 GPoint : (a : Type) -> (g : GetPoint a) => Type
128 | GPoint a = Point (gtrans @{g})
131 | 0 GVect : (a : Type) -> (g : GetPoint a) => Type
132 | GVect a = Vector (transform $
gtrans @{g})
135 | {t : _} -> GetPoint (Point t) where
142 | interface ModPoint (0 a : Type) where
143 | mtrans : AffineTransformation
144 | modPoint : (Point mtrans -> Point mtrans) -> a -> a
147 | {t : _} -> ModPoint (Point t) where
152 | 0 MPoint : (a : Type) -> (m : ModPoint a) => Type
153 | MPoint a = Point (mtrans @{m})
156 | 0 MVect : (a : Type) -> (m : ModPoint a) => Type
157 | MVect a = Vector (transform $
mtrans @{m})
161 | setPoint : ModPoint a => MPoint a -> a -> a
162 | setPoint = modPoint . const
166 | pointId : GetPoint a => a -> Point Id
167 | pointId = convert . point
171 | translate : ModPoint a => MVect a -> a -> a
172 | translate (V dx dy) = modPoint (\(P x y) => P (x + dx) (y + dy))
176 | scale : ModPoint a => Scale -> a -> a
177 | scale v = modPoint (\(P x y) => P (x * v.value) (y * v.value))
181 | rotateAt : ModPoint a => MPoint a -> Angle -> a -> a
182 | rotateAt o a = modPoint $
\p => let v := p-o in translate (rotate a v - v) p
186 | rotate : ModPoint a => Angle -> a -> a
187 | rotate = rotateAt origin
191 | distance : GetPoint a => a -> a -> Double
192 | distance x y = length (point x - point y)
196 | distanceFromZero : GetPoint a => a -> Double
197 | distanceFromZero = distance origin . point
201 | near : GetPoint a => (x,y : a) -> (delta : Double) -> Bool
202 | near x y = (distance x y <=)
204 | record Center2d where
210 | toCenter : Center2d -> Point t
211 | toCenter (C2D _ _ 0) = P 0 0
212 | toCenter (C2D x y n) = let dn := cast n in P (x/dn) (y/dn)
214 | addPoint : GetPoint a => Center2d -> a -> Center2d
215 | addPoint (C2D sx sy c) v = let P x y := point v in C2D (sx+x) (sy+y) (S c)
219 | center2d : GetPoint a => Foldable t => t a -> GPoint a
220 | center2d = toCenter . foldl addPoint (C2D 0 0 0)
228 | apply : Matrix -> Point t -> Point t
229 | apply (M a b c d) (P x y) = P (a * x + b * y) (c * x + d * y)
235 | perpendicularPoint :
237 | -> (p1,p2 : Point k)
238 | -> (distance : Double)
239 | -> (positive : Bool)
241 | perpendicularPoint p1 p2 d flag =
242 | let v := scaleTo d $
p2 - p1
243 | phi := if flag then halfPi else negate halfPi
244 | in translate (rotate phi v) p1
254 | intersect : (p11,p12,p21,p22 : Point t) -> Maybe (Point t)
255 | intersect (P x11 y11) (P x12 y12) (P x21 y21) (P x22 y22) =
256 | let p := P {t} (x21 - x11) (y21 - y11)
257 | m := M (x12 - x11) (x21 - x22) (y12 - y11) (y21 - y22)
258 | Just m' := inverse m | Nothing => Nothing
259 | P r s := apply m' p
260 | in Just $
P (x11 + r * (x12 - x11)) (y11 + r * (y12 - y11))
265 | distanceToLine : {t : _} -> (p, pl1, pl2 : Point t) -> Maybe Double
266 | distanceToLine p pl1 pl2 =
267 | let pp := perpendicularPoint p (translate (pl1 - pl2) p) 1 True
268 | in distance p <$> intersect pl1 pl2 p pp
280 | alignBond : {t : _} -> (pr1, pr2, p1, p2 : Point t) -> Point t -> Point t
281 | alignBond pr1 pr2 p1 p2 =
282 | let Just ar := angle (pr2 - pr1) | Nothing => id
283 | Just a := angle (p2 - p1) | Nothing => id
284 | dr := distance pr2 pr1
285 | d := distance p2 p1
286 | True := d > 0.0 | False => id
287 | sc := Scale.scale (dr/d)
289 | in \p => translate (origin - p1) p |> scale sc |> rotate da |> translate (pr1-origin)
300 | circularFreeSweep :
303 | -> {auto prf : IsSucc n}
307 | circularFreeSweep n p [] = (zero, fullSteps n)
308 | circularFreeSweep n p [x] = (angleOrZero $
x - p, fullSteps $
S n)
309 | circularFreeSweep n p ps@(x::_) =
310 | let ang := angleOrZero $
p - center2d ps
311 | angs := map (\y => angleOrZero $
y - p) ps
312 | in case enclosingAngles ang angs of
313 | Nothing => (angleOrZero $
x - p, fullSteps $
S n)
314 | Just (a1,a2) => (a1, divide (1+n) (a2-a1))
323 | perpendicularTo : {t : _} -> (x,y,c : Point t) -> Vector (transform t)
324 | perpendicularTo x y c =
325 | let cl := center2d (the (List _) [x,y])
327 | v := perpendicular (y-x)
328 | in if dot v vc >= 0 then v else negate v
332 | perpendicularFrom : {t : _} -> (x,y,c : Point t) -> Vector (transform t)
333 | perpendicularFrom x y = negate . perpendicularTo x y