content: {-|
   This is directly implemented from [\"Gear Drawing with Bézier Curves\" by Dr A R Collins](https://www.arc.id.au/GearDrawing.html).
  
   Found Indirectly via Mathew Dockrey 
   ( [Attoparsec](https://www.attoparsec.com/) 
   ( [the YouTube channel](https://www.youtube.com/attoparsec), not the parser library ) )
   and his [Inkscape Plugin](https://github.com/attoparsec/inkscape-extensions)
  
   And David Douard and Jonas Bähr, and their [FreeCAD module](https://github.com/FreeCAD/FreeCAD/blob/0ac0882eeb4e3390aef464e1807a3631c5f2e858/src/Mod/PartDesign/fcgear/involute.py)
  
  <<models/gear.glb>>
  -}
  module GearExample 
  ( gearExample
  ) where 
  
  import qualified Waterfall.Solids as Solids
  import qualified Waterfall.TwoD.Shape as Shape
  import qualified Waterfall.TwoD.Path2D as Path2D
  import Waterfall.TwoD.Transforms (rotate2D)
  import Control.Lens ((^.))
  import Linear (V2 (..), _y)
  import Data.Maybe (catMaybes)
  
  
  chebyExpnCoeffs :: Int -> (Double -> Double) -> Double
  chebyExpnCoeffs j f = 
      let n = 50 :: Int
          jf = fromIntegral j
          nf = fromIntegral n 
          c = sum [let kf = fromIntegral k in f (cos (pi * (kf - 0.5)/nf)) * cos (pi * jf * (kf - 0.5)/nf)| k <- [1..n]]
       in 2 * c / nf
  
  cheby :: [[Double]] 
  cheby = [ [ 1,  0,  0,  0,  0,  0]
          , [ 0,  1,  0,  0,  0,  0]
          , [-1,  0,  2,  0,  0,  0]
          , [ 0, -3,  0,  4,  0,  0]
          , [ 1,  0, -8,  0,  8,  0]
          , [ 0,  5,  0,-20,  0, 16]
          ]
  
  -- limited to p' = 5, but in practice p' = 4
  chebyApprox :: (Double -> Double) -> Int -> [Double]
  chebyApprox f p' = 
      let fnCoeffs = [chebyExpnCoeffs k f | k <- [0..p'] ]
          adjust 0 = head fnCoeffs /2
          adjust _ = 0
       in [ sum [fnCoeffs!!k  * (cheby !! k !! pwr) | k <- [0..p'] ] - adjust pwr | pwr <- [0..p'] ]
  
  binom :: Int -> Int -> Double 
  binom n k = fromIntegral (product [n - k + 1 .. n]) / fromIntegral (product [1..k])
  
  involuteBezCoeffs :: Double -> Double -> Double -> Double -> (V2 Double, V2 Double, V2 Double, V2 Double)
  involuteBezCoeffs rA rB fStart fStop = 
      let 
          p = 3
          ta = sqrt (rA * rA - rB * rB)/ rB -- involute angle at addendum
          ts = ta * sqrt fStart
          te = ta * sqrt fStop
  
          involuteXbez t =
              let x = t*2 -1
                  theta = x * (te - ts) / 2 + (ts + te)/2
              in rB * ( cos theta + theta * sin theta )
          involuteYbez t = 
              let x = t*2 - 1
                  theta = x * (te - ts) / 2 + (ts + te)/2
              in rB * ( sin theta - theta * cos theta )
  
          bezCoeff i f = 
              let polyCoeffs = chebyApprox f p
              in sum [binom i j * (polyCoeffs !! j) / binom p j | j<- [0..i]]
  
          v i = V2 (bezCoeff i involuteXbez) (bezCoeff i involuteYbez)
  
      in (v 0, v 1, v 2, v 3)
  
  genInvolutePolar :: Double -> Double -> Double
  genInvolutePolar rb r = let ra = sqrt (r * r - rb * rb) in ra/rb - acos (rb / r)
  
  polarToCart :: Double -> Double -> V2 Double
  polarToCart rad angle = V2 (rad * cos angle) (rad * sin angle)
  
  genGearToothData :: Double -> Int -> Double -> Path2D.Path2D
  genGearToothData m z phi = 
      let addendum = m
          dedendum = 1.25 * m
          clearance = dedendum - addendum
          rPitch = fromIntegral z  * m / 2
          rb = rPitch * cos phi
          ra = rPitch + addendum
          rRoot = rPitch - dedendum
          fRad = 1.5 * clearance 
          pitchAngle = 2 * pi / fromIntegral z
          baseToPitchAngle = genInvolutePolar rb rPitch
          rf' = sqrt ((rRoot + fRad) * (rRoot + fRad) - (fRad * fRad))
          rf = if rb < rf'
                  then rRoot + clearance  
                  else rf'
          pitchToFilletAngle  = 
              if rf > rb 
                  then baseToPitchAngle - genInvolutePolar rb rf 
                  else baseToPitchAngle 
          filletAngle = atan (fRad / (fRad + rRoot))
          fe = 1
          fs = if rf > rb 
              then (rf * rf - rb *rb) / (ra*ra - rb*rb)
              else 0.01 -- fraction of length offset from base to avoid singularity
          fm = fs + (fe - fs)/ 4
          (dbz1, dbz2, dbz3, dbz4) = involuteBezCoeffs ra rb fs fm 
          (_, abz2, abz3, abz4) = involuteBezCoeffs ra rb fm fe 
          rotateBez = rotate2D (-baseToPitchAngle-pitchAngle/4)
          rotateBez' = (* V2 1 (-1))  . rotateBez
          fillet = polarToCart rf (-pitchAngle / 4 - pitchToFilletAngle)
          arcMiddle = polarToCart ra 0
          filletR = (* V2 1 (-1)) fillet
          rootR = polarToCart rRoot (pitchAngle/4 +pitchToFilletAngle + filletAngle)
          rootNext = polarToCart rRoot (3*pitchAngle/4 - pitchToFilletAngle - filletAngle)
          filletNext = rotate2D pitchAngle fillet
  
      in Path2D.pathFrom fillet $ 
          catMaybes 
              [ if rf < rb 
                  then Just $ Path2D.lineTo (rotateBez dbz1)
                  else Nothing
              , Just $ Path2D.bezierTo (rotateBez dbz2) (rotateBez dbz3) (rotateBez dbz4)
              , Just $ Path2D.bezierTo (rotateBez abz2) (rotateBez abz3) (rotateBez abz4)
              , Just $ Path2D.arcViaTo arcMiddle (rotateBez' abz4)
              , Just $ Path2D.bezierTo (rotateBez' abz3) (rotateBez' abz2) (rotateBez' dbz4)
              , Just $ Path2D.bezierTo (rotateBez' dbz3) (rotateBez' dbz2) (rotateBez' dbz1)
              , if rf < rb 
                  then Just $ Path2D.lineTo filletR
                  else Nothing
              , if rootNext ^. _y > rootR ^. _y
                  then Just $ Path2D.pathFromTo 
                      [ Path2D.arcTo Path2D.Counterclockwise fRad rootR
                      , Path2D.arcTo Path2D.Counterclockwise rRoot rootNext -- these lines should be arcs
                      ]
                  else Nothing
              , Just $ Path2D.arcTo Path2D.Counterclockwise fRad filletNext
          ]
  -- Thickness, Module, Number Teeth, pressure Angle 
  gearExample :: Double -> Double -> Int -> Double -> Solids.Solid
  gearExample thickness moduleLength nGears pressureAngle =
      let segment = genGearToothData moduleLength nGears pressureAngle
          path = Path2D.repeatLooping segment
          --path = mconcat [rotate2D (-fromIntegral n * pi * 2 / fromIntegral nGears) segment | n <- [0..nGears]]
      in Solids.prism thickness . Shape.makeShape $ path