from rsf.proj import*


## parameters definition
clip=0.85               #display percentage
fraction=0.5    #B=fraction*I
niter=40        #number of iterations
n1=512          #temporal sampling number
n2=256          #spatial sampling number
padno=256       #padding for seislet tranform
r1=10           #smoothing radius
r2=10           #smoothing radius
or1=25
or2=25
oiter=5
ddip=3          #changing dip map interval
fhi=50          #bandpass frequency value
mute1='cp'      #muting function in shaping
mute2='cp'      #muting function in shaping
mode='soft'     #thresholding type (hard thresholding turns to be very bad
thr=30          #thresholding level(percentage)

## module defining
def Grey(hyper,other):
        Result(hyper,'put o2=-1 | grey label2=Trace unit2="" label1=Time labelsz=10 unit1="s" title="" wherexlabel=b wheretitle=t screenratio=1.3 %s'%other)

def Greyz(hyper,other):
        Result(hyper,'math output="input*2" | grey label2=Trace unit2="" label1=Time labelsz=10 unit1="s" title="" wherexlabel=b wheretitle=t screenratio=1.3 %s'%other)

def Greyplot(hyper,other):
        Plot(hyper,'grey label2=Trace unit2="" label1=Time labelsz=10 unit1="s" title="" wherexlabel=b wheretitle=t screenratio=1.3  %s'%other)

def Graph(hyper,other):
        Result(hyper,'graph label1="" label2="" unit1="" labelsz=10 unit2=""  title="" wherexlabel=b wheretitle=t %s' %other)

def Graphplot(hyper,other):
        Plot(hyper,'graph label1="" label2="" unit1="" labelsz=10 unit2=""  title="" wherexlabel=b wheretitle=t %s' %other)


##########################################
#    Make synthetic hyper:hyper1* & hyper2
##########################################
Flow('vrms',None,
     'math d1=0.004 n1=1001 o1=0 output="4500" ')

Flow('synt',None,
     '''
     spike d1=0.004 n1=1001 |
     noise rep=y seed=2006 |
     cut n1=100 | 
     bandpass flo=4 fhi=20 |
     spray axis=2 n=256 d=12.5 o=-1600 label=Offset unit=m 
     ''')

Flow('h','synt vrms',
     'inmo velocity=${SOURCES[1]} half=y | noise seed=2007 var=1e-10 | scale axis=2 | put o2=0 d2=1 | window n1=501')
Flow('h1','h','cp')
Flow('h2','h','cp')

Grey('h1','')
Grey('h2','')

#############################################
#               Experiment
#############################################
## Apply dithering
# var=1 makes the dithering range larger, unit=ms
Flow('dither','h1',
     '''
     window n1=1 |
     noise rep=y seed=122011 var=0.1 | math output="1000*input"
     ''')
Flow('hshottime1','h1','window n1=1 | math output=3*1000*x1')
Flow('hshottime2','hshottime1 dither','add scale=1,1 ${SOURCES[1]}')

## Blend 
Flow('hs','h2 h1 hshottime1 hshottime2','blend shot_time_in=${SOURCES[3]} shot_time_out=${SOURCES[2]} |add scale=1,1 ${SOURCES[1]}' )
Flow('uhs','h1 h2 hshottime1 hshottime2','blend shot_time_in=${SOURCES[2]} shot_time_out=${SOURCES[3]} |add scale=1,1 ${SOURCES[1]}' )

Grey('hs',' clip=%g'%clip)
Grey('uhs',' clip=%g'%clip)
Flow('hsft','hs','math output="input/2"')
Flow('uhsft','uhs','math output="input/2"')
Flow('hsst','hs','math output="input/2"')
Flow('uhsst','uhs','math output="input/2"')
Flow('hsor','hs','math output="input/2"')
Flow('uhsor','uhs','math output="input/2"')
Flow('hsfx','hs','math output="input/2"')
Flow('uhsfx','uhs','math output="input/2"')

## fk transform and filtering
Flow('hsfka','hs','fft1 | fft3 axis=2 pad=1 | cabs')
Flow('hsfkr','hs','fft1 | fft3 axis=2 pad=1 | real')
Flow('hsfki','hs','fft1 | fft3 axis=2 pad=1 | imag')
Flow('hsfk','hs','fft1 | fft3 axis=2 pad=1')
Flow('hsfkr_filt','hsfkr','mutter half=n t0=0 slope0=40 x0=0 ')
Flow('hsfki_filt','hsfki','mutter half=n t0=0 slope0=40 x0=0 ')

Grey('hsfka',' label1=Frequency unit1=Hz label2=Wavenumber color=j')
Grey('hsfkr',' label1=Frequency unit1=Hz label2=Wavenumber color=j')
Grey('hsfki',' label1=Frequency unit1=Hz label2=Wavenumber color=j')
Grey('hsfkr_filt',' label1=Frequency unit1=Hz label2=Wavenumber color=j')
Grey('hsfki_filt',' label1=Frequency unit1=Hz label2=Wavenumber color=j')

## fk inverse transform -> recon and difference
Flow('hs_re','hsfkr_filt hsfki_filt','cmplx ${SOURCES[1]} | fft3 axis=2 inv=y | fft1 inv=y')
Flow('hs_redif','hs hs_re','add scale=1,-1 ${SOURCES[1]}')
Grey('hs_re','')
Grey('hs_redif','')


def shapeslet(sig1,sig2,dip1,dip2,mute1,mute2,padno,n2,mode,thr,i): #mute1(2) is a string, when there is no mute, mute1(2) is "cp"
    Flow(sig1+"p%g"%(i+1),[sig1,dip1],
         '''
         %s |pad n2=%g  |
         seislet dip=${SOURCES[1]} eps=0.1 
                  adj=y inv=y unit=y type=b |
         threshold1 type=%s ifperc=1 thr=%f |
         seislet dip=${SOURCES[1]} eps=0.1  
                        inv=y unit=y type=b | 
         window n2=%g |%s
         '''%(mute1,padno,mode,thr,n2,mute1))
    Flow(sig2+"p%g"%(i+1),[sig2,dip2],
         ''' 
         %s|pad n2=%g  |
         seislet dip=${SOURCES[1]} eps=0.1 
                  adj=y inv=y unit=y type=b |
         threshold1 type=%s ifperc=1 thr=%f |
         seislet dip=${SOURCES[1]} eps=0.1  
                        inv=y unit=y type=b | 
         window n2=%g | %s
         '''%(mute2,padno,mode,thr,n2,mute2))
    return (sig1+"p%g"%(i+1),sig2+"p%g"%(i+1))

def shapesletor(sig1,sig2,blended1,blended2,diff1,diff2,dip1,dip2,mute1,mute2,padno,n2,mode,thr,i,or1,or2,oiter): #mute1(2) is a string, when there is no mute, mute1(2) is "cp"
    Flow(sig1+"p%g-t"%(i+1),[sig1,dip1],
         '''
         %s |pad n2=%g  |
         seislet dip=${SOURCES[1]} eps=0.1 
                  adj=y inv=y unit=y type=b |
         threshold1 type=%s ifperc=1 thr=%f |
         seislet dip=${SOURCES[1]} eps=0.1  
                        inv=y unit=y type=b | 
         window n2=%g |%s
         '''%(mute1,padno,mode,thr,n2,mute1))
    Flow(sig2+"p%g-t"%(i+1),[sig2,dip2],
         ''' 
         %s|pad n2=%g  |
         seislet dip=${SOURCES[1]} eps=0.1 
                  adj=y inv=y unit=y type=b |
         threshold1 type=%s ifperc=1 thr=%f |
         seislet dip=${SOURCES[1]} eps=0.1  
                        inv=y unit=y type=b | 
         window n2=%g | %s
         '''%(mute2,padno,mode,thr,n2,mute2))

    if(i+1<=oiter):
        Flow([diff1+"p%g"%(i+1),sig1+"p%g"%(i+1)],[blended1,sig1+"p%g-t"%(i+1),sig1+"p%g-t"%(i+1)],
         '''
         add scale=1,-1 ${SOURCES[1]} | 
         ortho niter=100 rect1=%g rect2=%g esp=0.01 
         sig=${SOURCES[1]} sig2=${TARGETS[1]}           
         '''%(or1,or2))
        Flow([diff2+"p%g"%(i+1),sig2+"p%g"%(i+1)],[blended2,sig2+"p%g-t"%(i+1),sig2+"p%g-t"%(i+1)],
         '''
         add scale=1,-1 ${SOURCES[1]} | 
         ortho niter=100 rect1=%g rect2=%g esp=0.01 
         sig=${SOURCES[1]} sig2=${TARGETS[1]}           
         '''%(or1,or2))
    else:
        Flow([diff1+"p%g"%(i+1)],[blended1,sig1+"p%g-t"%(i+1)],
         '''
         add scale=1,-1 ${SOURCES[1]}   
         ''')
        Flow([diff2+"p%g"%(i+1)],[blended2,sig2+"p%g-t"%(i+1)],
         '''
         add scale=1,-1 ${SOURCES[1]}   
         ''')
        Flow(sig1+"p%g"%(i+1),sig1+"p%g-t"%(i+1),'cp')
        Flow(sig2+"p%g"%(i+1),sig2+"p%g-t"%(i+1),'cp')

    return (sig1+"p%g"%(i+1),sig2+"p%g"%(i+1))

def shapefft(sig1,sig2,mute1,mute2,mode,thr,i):#mute1(2) is a string, when there is no mute, mute1(2) is "cp"
    Flow(sig1+"p%g"%(i+1),sig1,
         '''
         %s|fft1 | fft3 axis=2 pad=1|
         threshold1 type=%s ifperc=1 thr=%f |
         fft3 inv=y axis=2 | fft1 inv=y |%s
         '''%(mute1,mode,thr,mute1))
    Flow(sig2+"p%g"%(i+1),sig2,
         '''
         %s|fft1 | fft3 axis=2 pad=1|
         threshold1 type=%s ifperc=1 thr=%f |
         fft3 inv=y axis=2 | fft1 inv=y | %s
         '''%(mute2,mode,thr,mute2))
    return (sig1+"p%g"%(i+1),sig2+"p%g"%(i+1))

def shapefxdecon(sig1,sig2,mute1,mute2,n2,i): #mute1(2) is a string, when there is no mute, mute1(2) is "cp"
    Flow(sig1+"p%g"%(i+1),[sig1],
         '''
         %s | fxdecon n2w=%g 
         |%s
         '''%(mute1,n2,mute1))
    Flow(sig2+"p%g"%(i+1),[sig2],
         ''' 
         %s | fxdecon n2w=%g 
         |%s
         '''%(mute2,n2,mute2))
    return (sig1+"p%g"%(i+1),sig2+"p%g"%(i+1))

def step(sig1,sig2,blended1,blended2,init1,init2,shottime1,shottime2,i,n2,fraction):
#    Flow(init1+'sigtemp'+'%g'%(i+1),[sig1,sig2],'cat axis=2 ${SOURCES[1]}')
#    Flow(init1+'shottimetemp'+'%g'%(i+1),[shottime1,shottime2],'cat axis=2 ${SOURCES[1]}')
#    Flow(init1+'blendedtemp'+'%g'%(i+1),[init1+'sigtemp'+'%g'%(i+1),init1+'shottimetemp'+'%g'%(i+1)],
#       '''
#       blend shot_time_in=${SOURCES[1]} shot_time_out=${SOURCES[1]} 
#       ''')
#       
#    Flow(init1+"%g"%(i+1),[init1+'blendedtemp'+'%g'%(i+1),blended1,sig1], 
#         '''
#         window n2=%g | add scale=-1,1 ${SOURCES[1]} | add scale=%g,%g ${SOURCES[2]}
#         '''%(n2,fraction,1))
#    Flow(init2+"%g"%(i+1),[init1+'blendedtemp'+'%g'%(i+1),blended2,sig2], 
#         '''
#         window n2=%g f2=%g | add scale=-1,1 ${SOURCES[1]} | add scale=%g,%g ${SOURCES[2]}
#         '''%(n2,n2-1,fraction,1))

    Flow(init1+"%g"%(i+1),[sig2,shottime1,shottime2,blended1,sig1],
         '''
          blend shot_time_in=${SOURCES[2]} shot_time_out=${SOURCES[1]}
            | add scale=%f,%f,%f ${SOURCES[3]} ${SOURCES[4]}
         '''%(-fraction,fraction,1-fraction))

    Flow(init2+"%g"%(i+1),[sig1,shottime2,shottime1,blended2,sig2],
         '''
          blend shot_time_in=${SOURCES[2]} shot_time_out=${SOURCES[1]}
            | add scale=%f,%f,%f ${SOURCES[3]} ${SOURCES[4]}
         '''% (-fraction,fraction,1-fraction))
    return (init1+"%g"%(i+1),init2+"%g"%(i+1))

unblended1='h1'
unblended2='h2'
blended1='hs'
blended2='uhs'
init1='hsft'
init2='uhsft'
deblended1='hdbft1'
deblended2='hdbft2'
shottime1='hshottime1'
shottime2='hshottime2'
thr=20
sig1, sig2 = (init1,init2)
sigs1=[]
sigs2=[]
snrsa=[]
snrsb=[]
Flow('zero',unblended1,'math output=0')
Flow('energy1',[unblended1, 'zero'],'diff match=${SOURCES[1]}')
Flow('energy2',[unblended2, 'zero'],'diff match=${SOURCES[1]}')

for i in range(niter):
        old1,old2 = (sig1,sig2)
        snra=init1+'-snra%g'%i
        snrb=init2+'-snrb%g'%i

        Flow(snra,[unblended1,sig1,'energy1'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )
        Flow(snrb,[unblended2,sig2,'energy2'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )

        (nsig1,nsig2)=step(sig1,sig2,blended1,blended2,init1,init2,shottime1,shottime2,i,n2,fraction)
        (sig1,sig2)=shapefft(nsig1,nsig2,mute1,mute2,mode,thr,i)


        Greyplot(sig1,'title="Esti R 1 (iter=%g)" clip=%g'% (i+1,clip))
        Greyplot(sig2,'title="Esti R 2 (iter=%g)" clip=%g' % (i+1,clip))

        sigs1.append(sig1)
        sigs2.append(sig2)
        snrsa.append(snra)
        snrsb.append(snrb)

Flow(init1+'-snrsa',snrsa,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp | put n1=%g o1=0 d1=1 '%(len(snrsa),len(snrsa)))
Flow(init2+'-snrsb',snrsb,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp |  put n1=%g o1=0 d1=1'%(len(snrsb),len(snrsb)))

Graph(init1+'-snrsa','title="Data error 1" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)
Graph(init2+'-snrsb','title="Data error 2" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)

#Making movie
Plot(init1+'-sigs1',sigs1,'Movie')
Plot(init2+'-sigs2',sigs2,'Movie')
Flow(deblended1,sig1,'cp')
Flow(deblended2,sig2,'cp')

thr=32
unblended1='h1'
unblended2='h2'
blended1='hs'
blended2='uhs'
init1='hsst'
init2='uhsst'
deblended1='hdbst1'
deblended2='hdbst2'
shottime1='hshottime1'
shottime2='hshottime2'

Flow('mask',blended1,'math output=1 | pad n2=%g'%(padno))
sig1, sig2 = (init1,init2)
sigs1=[]
sigs2=[]
snrsa=[]
snrsb=[]
Flow('zero',unblended1,'math output=0')
Flow('energy1',[unblended1, 'zero'],'diff match=${SOURCES[1]}')
Flow('energy2',[unblended2, 'zero'],'diff match=${SOURCES[1]}')
for i in range(niter):
#######################################################################################
# update dip map every "ddip" iterations
#######################################################################################
        if i % ddip ==0 :
                dip1=init1+'dip%g'%int(i/ddip)
                dip2=init2+'udip%g'%int(i/ddip)
                Flow(dip1,[sig1,'mask'],
                        '''
                        bandpass fhi=%g | pad n2=%g | 
                        dip mask=${SOURCES[1]} rect1=%g rect2=%g
                        '''%(fhi,padno,r1,r2))
                Flow(dip2,[sig2,'mask'],
                        '''
                        bandpass fhi=%g | pad n2=%g | 
                        dip mask=${SOURCES[1]} rect1=%g rect2=%g
                        '''%(fhi,padno,r1,r2))
#######################################################################################
        old1,old2 = (sig1,sig2)
        snra=init1+'-snra%g'%i
        snrb=init2+'-snrb%g'%i

        Flow(snra,[unblended1,sig1,'energy1'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )
        Flow(snrb,[unblended2,sig2,'energy2'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )

        (nsig1,nsig2)=step(sig1,sig2,blended1,blended2,init1,init2,shottime1,shottime2,i,n2,fraction)
        (sig1,sig2)=shapeslet(nsig1,nsig2,dip1,dip2,mute1,mute2,padno,n2,mode,thr,i)



        Greyplot(sig1,'title="Esti R 1 (iter=%g)" clip=%g'% (i+1,clip))
        Greyplot(sig2,'title="Esti R 2 (iter=%g)" clip=%g' % (i+1,clip))

        sigs1.append(sig1)
        sigs2.append(sig2)
        snrsa.append(snra)
        snrsb.append(snrb)

Flow(init1+'-snrsa',snrsa,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp | put n1=%g o1=0 d1=1 '%(len(snrsa),len(snrsb)))
Flow(init2+'-snrsb',snrsb,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp |  put n1=%g o1=0 d1=1'%(len(snrsb),len(snrsb)))

Graph(init1+'-snrsa','title="Data error 1" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)
Graph(init2+'-snrsb','title="Data error 2" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)

#Making movie
Plot(init1+'-sigs1',sigs1,'Movie')
Plot(init2+'-sigs2',sigs2,'Movie')

Flow(deblended1,sig1,'cp')
Flow(deblended2,sig2,'cp')

thr=32
unblended1='h1'
unblended2='h2'
blended1='hs'
blended2='uhs'
diff1='hdifor1'
diff2='hdifor2'
init1='hsor'
init2='uhsor'
deblended1='hdbor1'
deblended2='hdbor2'
shottime1='hshottime1'
shottime2='hshottime2'

Flow('mask',blended1,'math output=1 | pad n2=%g'%(padno))
sig1, sig2 = (init1,init2)
sigs1=[]
sigs2=[]
snrsa=[]
snrsb=[]
Flow('zero',unblended1,'math output=0')
Flow('energy1',[unblended1, 'zero'],'diff match=${SOURCES[1]}')
Flow('energy2',[unblended2, 'zero'],'diff match=${SOURCES[1]}')
for i in range(niter):
#######################################################################################
# update dip map every "ddip" iterations
#######################################################################################
        if i % ddip ==0 :
                dip1=init1+'dip%g'%int(i/ddip)
                dip2=init2+'udip%g'%int(i/ddip)
                Flow(dip1,[sig1,'mask'],
                        '''
                        bandpass fhi=%g | pad n2=%g | 
                        dip mask=${SOURCES[1]} rect1=%g rect2=%g
                        '''%(fhi,padno,r1,r2))
                Flow(dip2,[sig2,'mask'],
                        '''
                        bandpass fhi=%g | pad n2=%g | 
                        dip mask=${SOURCES[1]} rect1=%g rect2=%g
                        '''%(fhi,padno,r1,r2))
#######################################################################################
        old1,old2 = (sig1,sig2)
        snra=init1+'-snra%g'%i
        snrb=init2+'-snrb%g'%i

        Flow(snra,[unblended1,sig1,'energy1'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )
        Flow(snrb,[unblended2,sig2,'energy2'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )

        (nsig1,nsig2)=step(sig1,sig2,blended1,blended2,init1,init2,shottime1,shottime2,i,n2,fraction)
        (sig1,sig2)=shapesletor(nsig1,nsig2,blended1,blended2,diff1,diff2,dip1,dip2,mute1,mute2,padno,n2,mode,thr,i,or1,or2,oiter)



        Greyplot(sig1,'title="Esti R 1 (iter=%g)" clip=%g'% (i+1,clip))
        Greyplot(sig2,'title="Esti R 2 (iter=%g)" clip=%g' % (i+1,clip))

        sigs1.append(sig1)
        sigs2.append(sig2)
        snrsa.append(snra)
        snrsb.append(snrb)

Flow(init1+'-snrsa',snrsa,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp | put n1=%g o1=0 d1=1 '%(len(snrsa),len(snrsb)))
Flow(init2+'-snrsb',snrsb,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp |  put n1=%g o1=0 d1=1'%(len(snrsb),len(snrsb)))

Graph(init1+'-snrsa','title="Data error 1" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)
Graph(init2+'-snrsb','title="Data error 2" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)

#Making movie
Plot(init1+'-sigs1',sigs1,'Movie')
Plot(init2+'-sigs2',sigs2,'Movie')

Flow(deblended1,sig1,'cp')
Flow(deblended2,sig2,'cp')

unblended1='h1'
unblended2='h2'
blended1='hs'
blended2='uhs'
init1='hsfx'
init2='uhsfx'
deblended1='hdbfx1'
deblended2='hdbfx2'
shottime1='hshottime1'
shottime2='hshottime2'

sig1, sig2 = (init1,init2)
sigs1=[]
sigs2=[]
snrsa=[]
snrsb=[]
Flow('zero',unblended1,'math output=0')
Flow('energy1',[unblended1, 'zero'],'diff match=${SOURCES[1]}')
Flow('energy2',[unblended2, 'zero'],'diff match=${SOURCES[1]}')
for i in range(niter):
#######################################################################################
        old1,old2 = (sig1,sig2)
        snra=init1+'-snra%g'%i
        snrb=init2+'-snrb%g'%i

        Flow(snra,[unblended1,sig1,'energy1'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )
        Flow(snrb,[unblended2,sig2,'energy2'],'diff match=${SOURCES[1]} | math a=${SOURCES[2]} output="a/input"' )

        (nsig1,nsig2)=step(sig1,sig2,blended1,blended2,init1,init2,shottime1,shottime2,i,n2,fraction)
        (sig1,sig2)=shapefxdecon(nsig1,nsig2,mute1,mute2,n2,i)


        Greyplot(sig1,'title="Esti R 1 (iter=%g)" clip=%g'% (i+1,clip))
        Greyplot(sig2,'title="Esti R 2 (iter=%g)" clip=%g' % (i+1,clip))

        sigs1.append(sig1)
        sigs2.append(sig2)
        snrsa.append(snra)
        snrsb.append(snrb)

Flow(init1+'-snrsa',snrsa,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp | put n1=%g o1=0 d1=1 '%(len(snrsa),len(snrsb)))
Flow(init2+'-snrsb',snrsb,'cat axis=2 ${SOURCES[1:%g]}|math output="10*log(input)/log(10)" | transp |  put n1=%g o1=0 d1=1'%(len(snrsb),len(snrsb)))

Graph(init1+'-snrsa','title="Data error 1" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)
Graph(init2+'-snrsb','title="Data error 2" label1=Iteration label2="c(n)" symbosz=5 symbol="o" min1=0 max1=%g d1=1'%niter)

#Making movie
Plot(init1+'-sigs1',sigs1,'Movie')
Plot(init2+'-sigs2',sigs2,'Movie')

Flow(deblended1,sig1,'cp')
Flow(deblended2,sig2,'cp')

## Ploting difference, error, deblended sections
Flow('hnft1','hs hdbft1','add scale=1,-1 ${SOURCES[1]}')
Flow('hnft2','uhs hdbft2','add scale=1,-1 ${SOURCES[1]}')
Flow('hnst1','hs hdbst1','add scale=1,-1 ${SOURCES[1]}')
Flow('hnst2','uhs hdbst2','add scale=1,-1 ${SOURCES[1]}')
Flow('hnor1','hs hdbor1','add scale=1,-1 ${SOURCES[1]}')
Flow('hnor2','uhs hdbor2','add scale=1,-1 ${SOURCES[1]}')
Flow('hnfx1','hs hdbfx1','add scale=1,-1 ${SOURCES[1]}')
Flow('hnfx2','uhs hdbfx2','add scale=1,-1 ${SOURCES[1]}')

Flow('heft1','h1 hdbft1','add scale=1,-1 ${SOURCES[1]}')
Flow('heft2','h2 hdbft2','add scale=1,-1 ${SOURCES[1]}')
Flow('hest1','h1 hdbst1','add scale=1,-1 ${SOURCES[1]}')
Flow('hest2','h2 hdbst2','add scale=1,-1 ${SOURCES[1]}')
Flow('heor1','h1 hdbor1','add scale=1,-1 ${SOURCES[1]}')
Flow('heor2','h2 hdbor2','add scale=1,-1 ${SOURCES[1]}')
Flow('hefx1','h1 hdbfx1','add scale=1,-1 ${SOURCES[1]}')
Flow('hefx2','h2 hdbfx2','add scale=1,-1 ${SOURCES[1]}')

Grey('hnft1','title="" clip=%g'%clip)
Grey('hnft2','title="" clip=%g'%clip)
Grey('hnst1','title="" clip=%g'%clip)
Grey('hnst2','title="" clip=%g'%clip)
Grey('hnor1','title="" clip=%g'%clip)
Grey('hnor2','title="" clip=%g'%clip)
Grey('hnfx1','title="" clip=%g'%clip)
Grey('hnfx2','title="" clip=%g'%clip)
Grey('heft1','title="" clip=%g'%clip)
Grey('heft2','title="" clip=%g'%clip)
Grey('hest1','title="" clip=%g'%clip)
Grey('hest2','title="" clip=%g'%clip)
Grey('heor1','title="" clip=%g'%clip)
Grey('heor2','title="" clip=%g'%clip)
Grey('hefx1','title="" clip=%g'%clip)
Grey('hefx2','title="" clip=%g'%clip)
#Grey('hdbft1','title="Deblended 1 (ft)"clip=%g'%clip)
#Grey('hdbft2','title="Deblended 1 (ft)"clip=%g'%clip)
#Grey('hdbst1','title="Deblended 1 (Seislet)" clip=%g'%clip)
#Grey('hdbst2','title="Deblended 1 (Seislet)" clip=%g'%clip)
Grey('hdbft1',' clip=%g'%clip)
Grey('hdbft2',' clip=%g'%clip)
Grey('hdbst1',' clip=%g'%clip)
Grey('hdbst2',' clip=%g'%clip)
Grey('hdbor1',' clip=%g'%clip)
Grey('hdbor2',' clip=%g'%clip)
Grey('hdbfx1',' clip=%g'%clip)
Grey('hdbfx2',' clip=%g'%clip)

## Ploting
Flow('hsnrsa','hsft-snrsa hsst-snrsa hsor-snrsa hsfx-snrsa','cat axis=2 ${SOURCES[1:4]}')
Flow('hsnrsb','uhsft-snrsb uhsst-snrsb uhsor-snrsb uhsfx-snrsb','cat axis=2 ${SOURCES[1:4]}')

Graph('hsnrsa','dash=0,1,0 title=""  symbol="o+@*" symbolsz=8 label1="Iteration no. #" label2="SNR" unit2="dB"  min1=0 max1=%g min2=0 max2=30 d1=1 screenratio=1 title="SNR Comparison" '%niter)
Graph('hsnrsb','dash=0,1,0 title=""  symbol="o+@*" symbolsz=8 label1="Iteration no. #" label2="SNR" unit2="dB"  min1=0 max1=%g min2=0 max2=30 d1=1 screenratio=1 title="SNR Comparison" '%niter)

Greyplot('hs',' title="Iter # = %g"'%(0))
Greyplot('uhs',' title="Iter # = %g"'%(0))
deblendfts1=['hs']
deblendfts2=['uhs']
deblendsts1=['hs']
deblendsts2=['uhs']
deblendfxs1=['hs']
deblendfxs2=['uhs']

for i in range(niter):
        deblendft1='hsft%gp%g'%(i+1,i+1)
        deblendft2='uhsft%gp%g'%(i+1,i+1)
        deblendst1='hsst%gp%g'%(i+1,i+1)
        deblendst2='uhsst%gp%g'%(i+1,i+1)
        deblendfx1='hsfx%gp%g'%(i+1,i+1)
        deblendfx2='uhsfx%gp%g'%(i+1,i+1)
        deblendfts1.append(deblendft1)
        deblendfts2.append(deblendft2)
        deblendsts1.append(deblendst1)
        deblendsts2.append(deblendst2)
        deblendfxs1.append(deblendfx1)
        deblendfxs2.append(deblendfx2)

        Greyplot(deblendft1,' title="Iter # = %g"'%(i+1))
        Greyplot(deblendft2,' title="Iter # = %g"'%(i+1))
        Greyplot(deblendst1,' title="Iter # = %g"'%(i+1))
        Greyplot(deblendst2,' title="Iter # = %g"'%(i+1))
        Greyplot(deblendfx1,' title="Iter # = %g"'%(i+1))
        Greyplot(deblendfx2,' title="Iter # = %g"'%(i+1))
Plot('deblendft1',deblendfts1,'Movie')
Plot('deblendft2',deblendfts2,'Movie')
Plot('deblendst1',deblendsts1,'Movie')
Plot('deblendst2',deblendsts2,'Movie')
Plot('deblendfx1',deblendfxs1,'Movie')
Plot('deblendfx2',deblendfxs2,'Movie')


Flow('h1-z','h1','window min1=1.1 max1=1.6 min2=90 max2=180')
Flow('hs-z','hs','window min1=1.1 max1=1.6 min2=90 max2=180')
Flow('hdbft1-z','hdbft1','window min1=1.1 max1=1.6 min2=90 max2=180')
Flow('hdbst1-z','hdbst1','window min1=1.1 max1=1.6 min2=90 max2=180')
Flow('hdbor1-z','hdbor1','window min1=1.1 max1=1.6 min2=90 max2=180')
Flow('hdbfx1-z','hdbfx1','window min1=1.1 max1=1.6 min2=90 max2=180')


Greyz('h1-z','screenratio=0.75 clip=%g'%clip)
Greyz('hs-z','screenratio=0.75 clip=%g'%clip)
Greyz('hdbft1-z','screenratio=0.75 clip=%g'%clip)
Greyz('hdbst1-z','screenratio=0.75 clip=%g'%clip)
Greyz('hdbor1-z','screenratio=0.75 clip=%g'%clip)
Greyz('hdbfx1-z','screenratio=0.75 clip=%g'%clip)

## Creating framebox
x=140
y=1.3
w=30
w1=0.2

Flow('framez.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        string.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('framez','framez.asc',
        '''
        dd type=complex form=native |
        graph min1=90 max1=180 min2=1.1 max2=1.6 pad=n plotfat=12 plotcol=3 screenratio=0.75
        wantaxis=n wanttitle=n yreverse=y 
        ''')

Result('hdbor1-z-0','Fig/hdbor1-z.vpl framez','Overlay')
Result('hdbst1-z-0','Fig/hdbst1-z.vpl framez','Overlay')

## Creating framebox
x=90
y=1.1
w=90
w1=0.5

Flow('frame.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        string.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame','frame.asc',
        '''
        dd type=complex form=native |
        graph min1=1 max1=256 min2=0 max2=2 pad=n plotfat=15 plotcol=2 screenratio=1.3 
        wantaxis=n wanttitle=n yreverse=y 
        ''')


# adding reference trace
Flow('trace.asc',None,
         '''
         echo %d 0 %d 1 n1=4 in=$TARGET data_format=ascii_float
         ''' % (25,25)) #25 is the trace number
Plot('trace','trace.asc',
         '''
         dd form=native type=complex | 
         graph min1=1 max1=256 min2=0 title="" wantaxis=n scalebar=n pad=n plotfat=8 dash=4 screenratio=1.3
         ''') #250 is the number of traces

Result('h1-1','Fig/h1.vpl frame trace','Overlay')
Result('hs-1','Fig/hs.vpl frame','Overlay')
Result('hdbft1-1','Fig/hdbft1.vpl frame trace','Overlay')
Result('hdbst1-1','Fig/hdbst1.vpl frame trace','Overlay')
Result('hdbor1-1','Fig/hdbor1.vpl frame trace','Overlay')
Result('hdbfx1-1','Fig/hdbfx1.vpl frame trace','Overlay')

Result('h1-0','Fig/h1.vpl frame','Overlay')
Result('hdbft1-0','Fig/hdbft1.vpl frame','Overlay')
Result('hdbst1-0','Fig/hdbst1.vpl frame','Overlay')
Result('hdbor1-0','Fig/hdbor1.vpl frame','Overlay')
Result('hdbfx1-0','Fig/hdbfx1.vpl frame','Overlay')

Flow('htc','h1 hdbft1 hdbst1 hdbor1 hdbfx1','cat axis=3 ${SOURCES[1:5]} | window n2=1 f2=25 min1=1.35 max1=1.38')
Graph('htc','label2=Amplitude unit2= label1=Time unit1=s screenratio=1 dash=0,3,5,7,8 plotfat=4 plotcol="7,3,5,2,1" symbolsz=12 max2=-0.1 min2=-0.5 title="Amplitude Comparison"')

End()