0 | module Geom.Vector
  1 |
  2 | import Derive.Prelude
  3 | import Data.Refined
  4 | import Geom.Angle
  5 | import Geom.Scale
  6 |
  7 | %default total
  8 | %language ElabReflection
  9 |
 10 | ||| Tries to find solutions to the quadratic equation
 11 | ||| `ax^2 + bx + c = 0`.
 12 | export
 13 | solveQuadratic : (a,b,c : Double) -> Maybe (Double,Double)
 14 | solveQuadratic a b c =
 15 |   let disc := b * b - 4 * a * c
 16 |       True := disc >= 0 && a > 0 | False => Nothing
 17 |       root := sqrt disc
 18 |    in Just (((-b) - root) / (2 * a), ((-b) + root) / (2 * a))
 19 |
 20 | --------------------------------------------------------------------------------
 21 | --          Linear Transformations
 22 | --------------------------------------------------------------------------------
 23 |
 24 | ||| A linear transformation in a two-dimensional vector space
 25 | ||| is a sequence of scalings and rotations.
 26 | |||
 27 | ||| By not using a matrix here we make sure that a transformation
 28 | ||| is invertible, without the risk of a division by zero.
 29 | ||| This might be a performance issue, but only if we create large
 30 | ||| trees of wild mixtures of linear transformations. However,
 31 | ||| in a UI we often work with sequences of the same type of
 32 | ||| transformations such as a sequence of scalings or a sequence of
 33 | ||| rotations. These are merged automatically in the `Semigroup`
 34 | ||| implementation.
 35 | public export
 36 | record LinearTransformation where
 37 |   constructor LT
 38 |   scale  : Scale
 39 |   rotate : Angle
 40 |
 41 | ||| The identity transformation
 42 | public export
 43 | Id : LinearTransformation
 44 | Id = LT 1.0 zero
 45 |
 46 | ||| A clockwise rotation by the given angle.
 47 | public export
 48 | rotation : Angle -> LinearTransformation
 49 | rotation = LT 1.0
 50 |
 51 | ||| A scaling by the given factor.
 52 | public export
 53 | scaling : Scale -> LinearTransformation
 54 | scaling = (`LT` zero)
 55 |
 56 | export
 57 | Semigroup LinearTransformation where
 58 |   LT s1 r1 <+> LT s2 r2 = LT (s1 * s2) (r1 + r2)
 59 |
 60 | export
 61 | Monoid LinearTransformation where neutral = Id
 62 |
 63 | ||| Computes the inverse of a linear transformation, so
 64 | ||| `x <+> inverse x` corresponse to the identity (modulo rounding
 65 | ||| errors).
 66 | export %inline
 67 | inverse : LinearTransformation -> LinearTransformation
 68 | inverse (LT s r) = LT (inverse s) (negate r)
 69 |
 70 | --------------------------------------------------------------------------------
 71 | --          Vectors and VectorSpaces
 72 | --------------------------------------------------------------------------------
 73 |
 74 | ||| A two-dimensional vector.
 75 | public export
 76 | record Vector (t : LinearTransformation) where
 77 |   constructor V
 78 |   x : Double
 79 |   y : Double
 80 |
 81 | %runElab deriveIndexed "Vector" [Show,Eq,Ord]
 82 |
 83 | ||| The zero vector.
 84 | export
 85 | vzero : Vector t
 86 | vzero = V 0 0
 87 |
 88 | ||| Unity vector along the x axis
 89 | export
 90 | vone : Vector t
 91 | vone = V 1 0
 92 |
 93 | ||| Utility to help with type inference when using vectors in `Id`
 94 | export %inline
 95 | vid : (x,y : Double) -> Vector Id
 96 | vid = V
 97 |
 98 | ||| Computes the length of a vector.
 99 | export
100 | length : Vector t -> Double
101 | length (V x y) = sqrt (x * x + y * y)
102 |
103 | ||| Vector addition.
104 | export
105 | (+) : Vector t -> Vector t -> Vector t
106 | V x1 y1 + V x2 y2 = V (x1+x2) (y1+y2)
107 |
108 | ||| Vector subtraction.
109 | export
110 | (-) : Vector t -> Vector t -> Vector t
111 | V x1 y1 - V x2 y2 = V (x1-x2) (y1-y2)
112 |
113 | ||| Inverts the given vector by negating its coordinates.
114 | export
115 | negate : Vector t -> Vector t
116 | negate (V x y) = V (-x) (-y)
117 |
118 | ||| Multiply a vector with a scalar.
119 | export
120 | scale : Double -> Vector t -> Vector t
121 | scale v (V x y) = V (v*x) (v*y)
122 |
123 | ||| Scales the vector to the given length.
124 | export
125 | scaleTo : Double -> Vector t -> Vector t
126 | scaleTo l v =
127 |   let lv := length v in if lv == 0.0 then v else scale (l/lv) v
128 |
129 | ||| Normalize a vector to length 1.
130 | export %inline
131 | normalize : Vector t -> Vector t
132 | normalize = scaleTo 1.0
133 |
134 | ||| Apply a linear transformation to a vector
135 | export
136 | transform : LinearTransformation -> Vector t -> Vector t
137 | transform (LT s r) (V x y) =
138 |   let co := s.value * cos r.value
139 |       si := s.value * sin r.value
140 |    in V (co * x - si * y) (si * x + co * y)
141 |
142 | ||| Define a vector by giving its length and its angle
143 | export
144 | polar : Scale -> Angle -> Vector t
145 | polar s a = transform (LT s a) vone
146 |
147 | ||| Rotate a vector by the specified number of degrees
148 | export
149 | rotate : Angle -> Vector t -> Vector t
150 | rotate = transform . rotation
151 |
152 | ||| Tries to calculate the angle of the given vector.
153 | ||| Fails in case this is the zero vector.
154 | export
155 | angle : Vector t -> Maybe Angle
156 | angle (V x y) =
157 |   if x == 0 && y == 0 then Nothing
158 |   else if x == 0 then Just $ if y > 0 then halfPi else negate halfPi
159 |   else
160 |     let phi := atan (y / x)
161 |      in Just $ if x > 0 then angle phi else angle (phi + pi)
162 |
163 | ||| Like `angle` but returns `zero` in case the vector is the empty vector.
164 | export %inline
165 | angleOrZero : Vector t -> Angle
166 | angleOrZero = fromMaybe zero . angle
167 |
168 | ||| From a given vector `v` computes a vector perpendicular
169 | ||| to `v` so that `dot v (perpendicular v) === 0.0` (up to rounding
170 | ||| errors).
171 | export
172 | perpendicular : Vector t -> Vector t
173 | perpendicular (V x y) = V (-y) x
174 |
175 | ||| Computes the dot product of two vectors.
176 | export
177 | dot : Vector t -> Vector t -> Double
178 | dot (V x1 y1) (V x2 y2) = x1*x2 + y1*y2
179 |