# coding: utf-8
import sys

def SPQ(elv,ELini,hmax):
    # To obtain the discharge of overflow crest from the water level
    ELcr=ELini
    b=20.5
    n=4
    x=0.1387
    Cd=1.971+0.498*x+6.63*x*x
    a=(1.6-Cd)/(Cd-3.2)
    h=elv-ELcr
    q=0.0
    if 0.0<h:
        if hmax<0.0: Hd=h
        if 0.0<hmax: Hd=hmax
        c=1.6*(1.0+2.0*a*h/Hd)/(1.0+a*h/Hd)
        kp=0.112703962703963*(h/Hd)**2-0.237885780885781*(h/Hd)+0.126496503496503
        ka=-0.0813160731673404*(h/Hd)**2+0.255143813465017*(h/Hd)
#        kp=0.02
#        ka=10.0*kp
        bx=b*float(n)-2.0*(float(n-1)*kp+ka)*h
        q=c*bx*h**1.5
    return q

def RET_V(nn,rh,rv,elv):
    # To obtain the cumulative volume from the water level
    for i in range(0,nn-1):
        if rh[i]<=elv and elv<=rh[i+1]: break
    if rh[nn-1]<elv: i=nn-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    v=(elv-b)/a
    return v

def RET_H(nn,rh,rv,v):
    # To obtain the water level from cumulative volume
    for i in range(0,nn-1):
        if rv[i]<=v and v<=rv[i+1]: break
    if rv[nn-1]<v: i=nn-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    elv=a*v+b
    return elv


#Main routine
param=sys.argv
fnameR1=param[1]
fnameR2=param[2]
fnameW =param[3]

# Input H-V data
rh=[]
rv=[]
fin=open(fnameR1,'r')
text=fin.readline()
text=text.strip()
text=text.split()
vcoef=float(text[0])
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    rh=rh+[float(text[0])]
    rv=rv+[float(text[1])*vcoef]
fin.close()
nn=len(rh)

# Input time sequence of inflow
ti=[]
q_in=[]
fin=open(fnameR2,'r')
text=fin.readline()
text=text.strip()
text=text.split()
ELini=float(text[0])
outlet=float(text[1])
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    ti=ti+[float(text[0])]
    q_in=q_in+[float(text[1])]
fin.close()
nd=len(ti)

hmax=-99.9
for iii in range(0,10):
    pEL=[]
    pQo=[]
    fout=open(fnameW,'w')
    # Initial outflow
    elv=ELini
    EL=elv
    VOL=RET_V(nn,rh,rv,elv)
    q_out=SPQ(elv,ELini,hmax)+outlet
    i=0
    iUD=0
    print('{0:>5s} {1:>5s} {2:>10s} {3:>10s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'.format('i','iUD','time','EL','EL-elv','VOL','q_in','q_out'),file=fout)
    print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i,iUD,ti[i],EL,EL-elv,VOL,q_in[i],q_out),file=fout)
    # Iterative calculation
    iUD=1
    dh=0.0001
    eps=0.0001
    itmax=int(1.0/dh)*10
    for i in range(0,nd-1):
        Qin=0.5*(q_in[i]+q_in[i+1])*(ti[i+1]-ti[i])*3600.0
        hh=0.0
        j=0
        while 1:
            j=j+1;
            hh=hh+float(iUD)*dh
            elv=EL;    q1=SPQ(elv,ELini,hmax)+outlet
            elv=EL+hh; q2=SPQ(elv,ELini,hmax)+outlet
            Qout=0.5*(q1+q2)*(ti[i+1]-ti[i])*3600.0
            R=VOL+Qin-Qout
            elv=RET_H(nn,rh,rv,R)
            if iUD==1 and j==10:
                if EL+hh>elv:
                    iUD=-1
                    hh=0.0
                    j=0
            if iUD==-1 and j==10:
                if EL+hh<elv:
                    iUD=1
                    hh=0.0
                    j=0
            if abs(EL+hh-elv)<eps: break
#            if itmax<j: break
        VOL=R    # Cumulative volume
        EL=EL+hh # Elevation
        q_out=q2 # Outflow
        print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i+1,iUD,ti[i+1],EL,EL-elv,VOL,q_in[i+1],q_out),file=fout)
        pEL=pEL+[EL]
        pQo=pQo+[q_out]
        sys.stdout.write('\r%d / %d' % (i+1,nd-1))
        sys.stdout.flush()
    fout.close()
    if abs(ELini+hmax-max(pEL))<0.001: break
    hmax=max(pEL)-ELini
    sys.stdout.write('\n')
    sys.stdout.write('Time: %10.3f\n'% max(ti))
    sys.stdout.write('h   : %10.3f\n'% hmax)
    sys.stdout.write('Qin : %10.3f %10.3f\n'% (min(q_in),max(q_in)))
    sys.stdout.write('Qout: %10.3f %10.3f\n'% (min(pQo),max(pQo)))
    sys.stdout.write('EL  : %10.3f %10.3f\n'% (min(pEL),max(pEL)))
    del pEL,pQo
sys.stdout.write('\n')
sys.stdout.write('Time: %10.3f\n'% max(ti))
sys.stdout.write('h   : %10.3f\n'% hmax)
sys.stdout.write('Qin : %10.3f %10.3f\n'% (min(q_in),max(q_in)))
sys.stdout.write('Qout: %10.3f %10.3f\n'% (min(pQo),max(pQo)))
sys.stdout.write('EL  : %10.3f %10.3f\n'% (min(pEL),max(pEL)))