Akima 補間 途中まで書いて放置

akima.hs

type ZipList = [(Double, Double)]

-- akima interpolation
akimaint :: ZipList -> Double -> Maybe Double
akimaint zls x = akimacore zls' x
    where
        zls' = zlistsort zls
        zls'' = selzls zls' n n2
        xls = map fst zls'
        (n,n2) = indsearch xls x

akimacore :: ZipList -> Double -> Maybe Double
akimacore [(x1,y1), (x2,y2), (x3,y3), (x4,y4), (x5,y5), (x6,y6)] x
    = if and [x3 <= x, x4 >= x]
        then Just $ p0 + ax * (p1 + ax * (p2 + ax * p3))
        else Nothing
    where
        p0 = y3
        p1 = t1
        p2 = (3.0 * (y4 - y3) / (x4 - x3) - 2.0 * t1 - t2) / (x4 - x3)
        p3 = (t1 + t2 - 2.0 * (m (x3,y3) (x4,y4))) / ((x4 - x3)^2)
        t1 = t (m (x1,y1) (x2,y2)) (m (x2,y2) (x3,y3)) (m (x3,y3) (x4,y4)) (m (x4,y4) (x5,y5))
        t2 = t (m (x2,y2) (x3,y3)) (m (x3,y3) (x4,y4)) (m (x4,y4) (x5,y5)) (m (x5,y5) (x6,y6))
        t m0 m1 m2 m3 = (abs (m3 - m2) * m1 + abs (m1 - m0) * m2) / (abs (m3 - m2) + abs (m1 - m0))
        m (x,y) (x',y') = (y' - y) / (x' - x)
        ax = (x - x3)
akimacore _ _ = Nothing

zlistsort :: ZipList -> ZipList
zlistsort [] = []
zlistsort (x:xs) = (zlistsort st) ++ [x] ++ (zlistsort lt)
    where
        st = filter ((>)  (fst x) . fst) xs
        lt = filter ((<=) (fst x) . fst) xs

indsearch :: [Double] -> Double -> (Int, Int)
indsearch xs a = indsearch' xs a 0
    where
        indsearch' [] a n = (n-1, n)
        indsearch' (x:xs) a n 
            | a > x     = indsearch' xs a (n+1)
            | otherwise = (n-1,n)

selzls :: ZipList -> Int -> Int -> ZipList
selzls (z:zs) n n2 | and [n >= 2, ]

-- for test
zlist :: ZipList
zlist = zip [0.1, 0.2, 0.3, 0.5, 0.7, 0.8] [8.0, 5.0, 3.0, 5.0, 2.0, 3.0]

main :: IO ()
main = do
    putStrLn $ show $ akimaint zlist 0.4