# Define delay differential equation for change solute concentrations with
# respect to stream length, assuming steady state with respect to time.
macro stream_model(reactionKinetics, shape)
    if shape == "powerLaw"
        if reactionKinetics == "firstOrder"
        return :( 
            function stream_model!(du,u,p,t)
                τ₋, τ₊, α, k, qPerUnitHZ, sHZ, sCH, v = p
                integral = u[1] * quadgk(τ -> E(τ, α, τ₋, τ₊)*firstOrder(τ, k),τ₋, τ₊)[1]
                du[1] = (qPerUnitHZ*(sHZ/sCH)*(integral - u[1]))/v
            end
                )
        elseif reactionKinetics == "firstOrderProductionAndConsumption"
        return :( 
            function stream_steady_state_withNitrification!(du,u,p,t)
                τ₋, τ₊, α, qPerUnitHZ, sHZ, sCH, v, kprod, kcons, substrate0 = p
                integral =  quadgk(τ -> E(τ, α, τ₋, τ₊)*firstOrderProductionAndConsumption(τ, kprod, kcons, substrate0, u[1]),τ₋, τ₊)[1]
                du[1] = (qPerUnitHZ*(sHZ/sCH)*(integral - u[1]))/v
            end
                )
        end
    end
  end