module TruncatedSine where

{- Fourier series of a sine of fractional frequency -}

import Graphics.Gnuplot.Simple (plotFuncs, linearScale, )
import qualified FourierSeries as FS

{- these are sine series -}
fourierSeriesRecip, fourierSeriesSincc :: Floating a => a -> [a]
fourierSeriesRecip = fourierSeries recip
fourierSeriesSincc = fourierSeries sincc1

fourierSeriesOddRat :: Floating a => a -> a -> [a]
fourierSeriesOddRat a = fourierSeries (oddRat a)

fourierSeries :: Floating a => (a -> a) -> a -> [a]
fourierSeries rec w =
   zipWith (\d0 d1 -> rec d0 - rec d1)
           (iterate (2+) w) (iterate (subtract 2) w)

{- This continously differentiable function
   is 1/x for all integral x, but 0 at position 0 -}
sincc1 :: Floating a => a -> a
sincc1 0 = 0
sincc1 t =
   let pit = pi*t
   in  (pit - sin pit) / (pit*t)

oddRat :: Fractional a => a -> a -> a
oddRat a t = a*t / (1+a*t*t)

plotRecip, plotSincc, plotSinccAngle, plotOddRat :: IO ()
plotRecip =
   plotFuncs [] (linearScale 300 (0, 2*pi))
                (map (\s -> FS.evalSineL (take 1000
                               (fourierSeriesRecip (s::Double))))
                     (linearScale 6 (1,3)))

plotSincc =
   plotFuncs [] (linearScale 300 (0, 2*pi))
                (map (\s -> FS.evalSineL (take 1000
                               (fourierSeriesSincc (s::Double))))
                     (linearScale 6 (1,3)))

plotSinccAngle =
   plotFuncs [] (linearScale 300 (0, 2*pi))
                (map (\s t -> acos (2/pi *
                             FS.evalSineL (take 1000
                               (fourierSeriesSincc (s::Double))) t) {- - t/2 -})
                     (linearScale 6 (1,3)))


plotOddRat =
   plotFuncs [] (linearScale 300 (0, 2*pi))
                (map (\s -> FS.evalSineL (take 1000
                               (fourierSeriesOddRat 2 (s::Double))))
                     (linearScale 6 (1,3)))