monte carlo 1d radiative transfer code written in haskell

mc1drad.hs

{-# LANGUAGE ExtendedDefaultRules, NoMonomorphismRestriction #-}

import Conduit

np  = 10000
tau = 1.0 -- cloud optical thickness
sal = 0.5 -- surface albedo
omg = 1.0 -- single scattering albedo
zu  = tau -- atmosphere upper boundary
zb  = 0.0 -- atmosphere lower boundary

main :: IO ()
main = do
    (sumr, sumt) <- world
    print $ [sumr / np, sumt / np]

world = loop 0 0.0 0.0
    where
        loop t sumr sumt
            | t == np = return (sumr, sumt)
            | otherwise = do
                (wr, wt) <- sourceRandom $$ solver 1.0 zu (-1.0) 0.0 0.0
                loop (t+1) (sumr + wr) (sumt + wt)

solver :: (Monad m, MonadIO m) =>
          Double -> Double -> Double
          -> Double -> Double
          -> Sink Double m (Double, Double)
solver w z dir sumr sumt = do
    rand <- await
    case rand of
        Just rand' -> do
            path <- return $ -log(rand')
            nextphase $ z + dir * path
    where
        nextphase zn
            | zn <= zb = solver (w * sal) zb 1.0 sumr (sumt + w)
            | zn > zu = return (sumr + w, sumt)
            | otherwise  = do
                rand2 <- await
                case rand2 of
                    Just rand2' -> do
                        dnew <- return $ if rand2' <= 0.5
                                           then 1.0 * dir
                                           else -1.0 * dir
                        solver (w * omg) zn dnew sumr sumt