«K12/K13» sim.hs «PDF», «POD»

module Main where

import qualified Data.Map as Map
import qualified Data.Set as Set
import Data.Maybe
import Data.List (sort)
import Control.Monad
import System.IO

data Prep = Prep
  { name :: String
  , nz   :: (Int,Int)
  , half :: Double
  , decs :: Map.Map Prep Double
  } deriving (Show, Eq, Ord)

cores :: Set.Set Prep
inits :: Map.Map Prep Double
(cores,inits) =
  ( Set.fromList [ur1, kind1, kind2, ur2, stab1, stab2]
  , Map.fromList [(ur1, 1000), (ur2, 400), (kind2, 800)]
  )
  where
  ur1   = Prep "A-ur1"   (4,3)  10 $ Map.fromList [(kind1, 0.9), (kind2, 0.1)]
  ur2   = Prep "B-ur2"   (3,4)   2 $ Map.fromList [(kind2, 0.9), (stab1, 0.1)]
  kind1 = Prep "C-kind1" (2,3)   5 $ Map.fromList [(kind2, 1.0)]
  kind2 = Prep "D-kind2" (2,2)  20 $ Map.fromList [(stab1, 0.7), (stab2, 0.3)]
  stab1 = Prep "E-stab1" (1,2) inf $ Map.fromList []
  stab2 = Prep "E-stab2" (1,1) inf $ Map.fromList []

mothers :: Prep -> Set.Set (Prep,Double)
mothers child = Set.map (\m -> (m,fromJust $ child `Map.lookup` (decs m))) $ Set.filter ((child `Map.member`) . decs) cores

numDecays :: Prep -> Double -> Double
numDecays prep t = base - base * exp2 (-t/(half prep))
  where
  base = numAdds prep t

numAdds :: Prep -> Double -> Double
numAdds prep t =
  (maybe 0 id $ prep `Map.lookup` inits) +
  (sumS $ Set.map (\(m,p) -> p * numDecays m t) $ mothers prep)

num :: Prep -> Double -> Double
num prep t = numAdds prep t - numDecays prep t

graph :: Double -> Double -> Double -> [(Double,Map.Map Prep Double)]
graph t0 t1 dt = map (\t -> (t,step t)) $ fromToStep t0 t1 dt
  where
  step t = Set.fold (\core -> Map.insert core (num core t)) Map.empty cores

pretty :: Map.Map Prep Double -> String
pretty = concatMap (\(prep,count) -> show (fst $ nz prep) ++ " " ++ show (snd $ nz prep) ++ " " ++ show count ++ "\n") . valueList . Map.mapWithKey (\prep count -> (prep,count))

main = do
  forM (graph 0 200 0.1) $ \(t,g) -> do
    let fn = "nuclear-" ++ pad10 (show $ round $ t * 100)
    putStrLn $ fn ++ "..."
    writeFile fn $ pretty g
  where
  pad10 xs | length xs == 10 = xs
  pad10 xs | otherwise       = pad10 ('0':xs)

exp2 x = exp $ x * log 2

sumS :: (Num a) => Set.Set a -> a
sumS  = Set.fold (+) 0

valueList :: (Ord v, Ord k) => Map.Map k v -> [v]
valueList = map snd . sort . Map.toList

fromToStep x y d = map (\s -> x + d*s) [0..((y-x) / d)]

forM = flip mapM

inf = 1/0