up [pdf]
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()

sfmath
sfspike
sfnoise
sfcut
sfbandpass
sfspray
sfinmo
sfscale
sfput
sfwindow
sfcp
sfgrey
sfadd
sfblend
sffft1
sffft3
sfcabs
sfreal
sfimag
sfmutter
sfcmplx
sfdiff
sfthreshold1
sfcat
sftransp
sfgraph
sfpad
sfdip
sfseislet
sfortho
sffxdecon
sfdd