def queuePerformance(arrivalRate, serviceRate, servers=1):
    """
    Determine the average performance characteristics of a service queue using the M/M/c forumla

    arrivalRate = The rate of new customers in the queue
    serviceRate = The number of customers served in a time interval by all servers
    servers = The number of servers available (defaults to 1)

    returns a dictionary containing the following:

    emptyProbability = The probability the queue will be empty
    fullProbability = The probability the queue will be full
    averageLineLength = The average length of the queue
    averagePopulation = Average number of customers in the system as a whole (waiting and being served)
    averageWaitTime = Average time to wait until service
    averageTime = Average time in the system as a whole

    Doctests follow. To use these, run python interactively then type:
    import mmc
    import doctest
    doctest.testmod(mmc)

    Should return (0, 3) (no failures, 3 tests run)


    >>> queuePerformance(1, 1.1, 1)
    {'averageWaitTime': 9.0909090909090882, 'averageLineLength': 9.0909090909090882, 'averageTime': 9.9999999999999964, 'emptyProbability': 0.090909090909090939, 'averagePopulation': 9.9999999999999964, 'fullProbability': 0.082644628099173584}

    >>> queuePerformance(5, 6.1, 5)
    {'averageWaitTime': 6.3714079566131468e-05, 'averageLineLength': 0.00031857039783065737, 'averageTime': 0.16399814030907434, 'emptyProbability': 0.44055114515780114, 'averagePopulation': 0.81999070154537168, 'fullProbability': 0.0013583632864877702}

    >>> queuePerformance(3.4,4.2,3)
    {'averageWaitTime': 0.0058284681698965819, 'averageLineLength': 0.019816791777648377, 'averageTime': 0.24392370626513465, 'emptyProbability': 0.44281457712583089, 'averagePopulation': 0.82934060130145781, 'fullProbability': 0.039152503642860853}

    """
    # Input validation
    arrivalRate = float(arrivalRate)
    serviceRate = float(serviceRate)
    throughput = arrivalRate/serviceRate

    if arrivalRate >= serviceRate:
        raise ValueError("Arrival rate must be lower than Service rate, otherwise the result is an infinite wait")

    if servers < 1:
        raise ValueError("Servers must be greater than 1, otherwise the result is an infinite wait (or other strangeness)")

    """
    Rather than a for loop, I have used a list comprehension with a sum. I did this because it's fun, but it is functionally
    equivalent to:

    serviceRateSum = 0.0
    for server in range(servers):
        serviceRateSum += ((throughput ** server) / factorial(server))
    """

    serviceRateSum = sum([((throughput ** server) / factorial(server)) for server in range(servers)])

    # I really don't understand this number
    pwait = (throughput ** servers) / (factorial(servers) * (1 - throughput/servers))

    # Work out the queue characteristics
    emptyProbability = 1.0/(serviceRateSum + pwait)
    fullProbability = emptyProbability * throughput**servers / factorial(servers)
    averageLineLength = fullProbability * (throughput / servers) / ((1-throughput / servers) ** 2)
    averagePopulation = averageLineLength + throughput
    averageWaitTime = averageLineLength / arrivalRate
    averageTime = averagePopulation / arrivalRate

    return dict(emptyProbability=emptyProbability, fullProbability=fullProbability, averageLineLength=averageLineLength, averagePopulation=averagePopulation, averageWaitTime=averageWaitTime, averageTime=averageTime)

def mmc(lam, mu, c=1):
    """
    Interface wrapper for queuePerformance
    """
    try:
        result = queuePerformance(lam, mu, c)
    except ValueError:
        return dict()

    return dict(
        lam=lam,
        mu=mu,
        rho=lam/mu,
        c=c,
        L=result['averagePopulation'],
        W=result['averageTime'],
        Lq=result['averageLineLength'],
        Wq=result['averageWaitTime'],
        P0=result['emptyProbability'],
        Pc=result['fullProbability']
    )