{-|
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
, playgroundGearExample
) 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 p
| p == 0, (c : _) <- fnCoeffs = c / 2
| otherwise = 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
-- | This is intended to be viewed via the WASM playground
playgroundGearExample :: Solids.Solid
playgroundGearExample = gearExample 10 5 20 (20*pi/180)