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

mc1drad_mersenne.hs

import System.Random.Mersenne

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

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

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