content: 
  module Solve.Density
          (densitySteps)
  where
  import Stage.Diffusion
  import Stage.Advection
  import Stage.Sources
  import Config
  import Model
  
  -- | Run the stages for processing the density field in one time step
  densitySteps
          :: Config
          -> Int
          -> DensityField
          -> Maybe (SourceDensity Float)
          -> VelocityField
          -> IO DensityField
  
  densitySteps config step df ds vf
   = {-# SCC "Solve.densitySteps" #-}
     do   df1     <- addSources   (configDelta config) (configDensity   config)
                                  ds df
  
          let diff = if  configDiffAfter config /= 0
                      && step >= configDiffAfter config
                       then 0.0005
                       else configDiff config
  
          let diffSolver
                  = if configUnstable config
                          then DiffUnstable
                          else DiffStable (configIters config)
  
          df2     <- diffusion    diffSolver (configDelta config) diff
                                  df1
  
          df'     <- advection    (configDelta config) vf df2
  
          return  df'
  