monte carlo 1d radiative transfer code written in haskell using random-mwc package.

mc1drad_mwc.hs

import System.Random.MWC

np  = 1000000
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
    seed0 <- createSystemRandom >>= save
    (sumr, sumt) <- world seed0
    print $ [sumr / np, sumt / np]

world :: Seed -> IO (Double, Double)
world seed0 = loop 0 0.0 0.0 seed0
    where
        loop :: Double -> Double -> Double -> Seed -> IO (Double, Double)
        loop t sumr sumt seed
            | t == np = return (sumr, sumt)
            | otherwise = do
                ((wr, wt), newseed') <- solver 1.0 zu (-1.0) 0.0 0.0 seed
                loop (t+1) (sumr + wr) (sumt + wt) newseed'

solver :: Double -> Double -> Double
          -> Double -> Double
          -> Seed -> IO ((Double, Double), Seed)
solver w z dir sumr sumt seed = do
    (n, newseed) <- pureRnd seed
    path <- return $ -log(n)
    nextphase (z + dir * path) newseed
    where
        nextphase zn newseed
            | zn <= zb = solver (w * sal) zb 1.0 sumr (sumt + w) newseed
            | zn > zu = return ((sumr + w, sumt), newseed)
            | otherwise  = do
                (n2, newseed2) <- pureRnd newseed
                dnew <- return $ if n2 <= 0.5
                                   then 1.0 * dir
                                   else -1.0 * dir
                solver (w * omg) zn dnew sumr sumt newseed2

pureRnd :: Seed -> IO (Double, Seed)
pureRnd sd = do
    gen'' <- restore sd
    ret <- uniformR (0, 1) gen''
    seed' <- save gen''
    return (ret, seed')