Something

import numpydef iterate_scenario_1(Q23i, QX, V2i, VAL):    if Q23i > 0.001:        if V2i > VAL/2.:            return QX        else:            return 0    elif Q23i < 0.001:        if V2i < 2 * VAL:            return 0        else:            return QXdef set_QX(V24i, C24i, k24, QX1, QX2):    if V24i < 1.5 * C24i/k24: return QX1    else: return QX2def model_sewer(**kwargs):    iterations = int(kwargs.get('iterations', 800))    try: scenario   = int(kwargs.get('scenario', 1))    except: scenario = 1    CM12 = float(kwargs.get('CM12', 25))    CM24 = float(kwargs.get('CM24', 30))    k12  = float(kwargs.get('k12', 70))    k24  = float(kwargs.get('k24', 30))    VAL  = float(kwargs.get('VAL', 1))    QX1  = float(kwargs.get('QX1', 60))    QX2  = float(kwargs.get('QX2', 120))    b12 = float(kwargs.get('b12', 10./k12))    b24 =  float(kwargs.get('b24', 5./k24))        V12 = numpy.zeros((iterations,))    V24 = numpy.zeros((iterations,))    V2  = numpy.zeros((iterations,)); V2[0] = VAL    C12 = numpy.zeros((iterations,)); C12[0] = CM12    C24 = numpy.zeros((iterations,)); C24[0] = CM24    Q12 = numpy.zeros((iterations,))    Q23 = numpy.zeros((iterations,))    Q24 = numpy.zeros((iterations,))    R12 = numpy.ones((iterations,)) * 20.    r0 = int(kwargs.get('R12_start', 36))    r1 = r0 + int(kwargs.get('R12_duration', 72))    R12[r0:r1] = float(kwargs.get('R12_amount', 50.))    R24 = numpy.ones((iterations,)) * 10.    r0 = int(kwargs.get('R24_start', 36))    r1 = r0 + int(kwargs.get('R24_duration', 72))    R24[r0:r1] = float(kwargs.get('R24_amount', 30.))    QX = QX1    delta_t = 1./144    for i in range(iterations - 1):        if V12[i] < C12[i]/k12:            Q12[i] = k12 * V12[i]         else:            Q12[i] = C12[i]                if V24[i] < C24[i]/k24:            Q24[i] = k24 * V24[i]        else:            Q24[i] = C24[i]        if scenario == 2:            QX = set_QX(V24[i], C24[i], k24, QX1, QX2)        Q23[i + 1] = iterate_scenario_1(Q23[i], QX, V2[i], VAL)        # ADDED: max() so cant go below zero.        V12[i + 1] = V12[i] + (R12[i] - Q12[i]) * delta_t        V24[i + 1] = V24[i] + (R24[i] - Q24[i]) * delta_t        V2[i + 1]  = max(V2[i] + (Q12[i] + Q24[i] - Q23[i]) * delta_t, 0.)        if V2[i + 1] < VAL:            C12[i + 1] = CM12            C24[i + 1] = CM24        else:            C12[i + 1] = CM12 - b12 * ( V2[i + 1] - VAL )            C24[i + 1] = CM24 - b24 * ( V2[i + 1] - VAL )    if kwargs.get('scale'):        Q12/=CM12; Q23/=100; Q24/=CM24        V2/=5.    # NOTE V2 gets scaled, but note V12, V24    return {'L12': Q12, 'L23': Q23, 'L24': Q24}\            , {'V12': V12, 'V2': V2, 'V24': V24} \            , {'R12': R12, 'R24': R24 }import pylabdef plot(qs, vs, rs, filename=None):    colors = {'12': 'r', '23': 'g', '24': 'b', '2': 'k'}    pylab.figure(figsize=(5,6))    pylab.subplot(311)    labels = sorted(rs.keys())    for l in labels:        pylab.plot(rs[l], c=colors[l[1:]])    pylab.legend(labels)    pylab.ylim(0, 60)    pylab.ylabel('Runoff (Mg/day)')    labels = sorted(qs.keys())    pylab.subplot(312)    for l in labels:        pylab.plot(qs[l], c=colors[l[1:]])    pylab.legend(['Q' + l[1:] for l in labels])    pylab.ylim(0, 100)    pylab.ylabel('Flow (Mg/day)')    pylab.subplot(313)    labels = sorted(vs.keys())    for l in labels:        pylab.plot(vs[l], c=colors[l[1:]])    pylab.legend(labels)    pylab.ylim(0, 20)    pylab.ylabel('Volume (Mg)')    pylab.xlabel('T')    pylab.subplots_adjust(top=.99)    if not filename:        pylab.show()    else:        r = random.random()        fn = '/var/www/ms_tmp/model%f.png' % r        pylab.savefig(fn)        pylab.close()        return fnif __name__ == "__main__":    import pylab    qs, vs, rs = model_sewer(C23=3, iterations=144, scenario=1)    img = plot(qs, vs, rs, False)