up [pdf]
from rsf.proj import *
from rsf.recipes.beg import server
import math

sliceval = 0.222#0.218

def stackplot(name,stack,title,frame1=100,frame2=300,frame3=4,extra=''):
    Result(name,stack,
           '''
               transp plane=12 | transp plane=23 |
           byte gainpanel=all | 
           grey3 frame1=%d frame2=%d frame3=%d title="%s"
           point1=0.6 point2=0.9 screenratio=0.4 wanttitle=n
               d2num=0.5 n2tic=4 o2num=1.0
                   d1num=0.5 n1tic=7 o1num=4.0
                   d3num=0.03 n3tic=2 o3num=0.21
                   d4num=0.02 n4tic=1 o4num=0.23
               labelsz=14.0 crowd2=0.65 %s
           ''' % (frame1,frame2,frame3,title,extra))

def stackplotc(name,stack,title,frame1=100,frame2=300,frame3=4,extra=''):
    Result(name,stack,
           '''
           transp plane=12 | transp plane=23 |
           byte gainpanel=all clip=0.0356745 | 
           grey3 frame1=%d frame2=%d frame3=%d title="%s"
           point1=0.6 point2=0.9 screenratio=0.4 wanttitle=n
                   d2num=0.5 n2tic=4 o2num=1.0
                   d1num=0.5 n1tic=7 o1num=4.0
                   d3num=0.03 n3tic=2 o3num=0.21
                   d4num=0.02 n4tic=1 o4num=0.23
               labelsz=14.0 crowd2=0.65 %s
           ''' % (frame1,frame2,frame3,title,extra))

def sliceplot(name,slic,title,min1=sliceval,extra=''):
    Result(name,slic,
           '''
           window n1=1 min1=%g |
           grey title="%s" wanttitle=n
           screenratio=0.4
                   d1num=0.5 n1tic=8 o1num=4.0
                   labelsz=12.0
                   %s
           ''' % (min1,title,extra))

# | transp

def halfsliceplot(name,slic,title,min1=sliceval,min2=1.4,max2=1.9,min3=5.4,max3=6.0,extra=''):
    Result(name,slic,
           '''
           put label2=iline label3=xline | window n1=1 min1=%g min2=%g max2=%g min3=%g max3=%g |
           grey title="%s" wanttitle=n
           screenratio=0.6
                   d2num=0.2 n2tic=3 o2num=1.45 labelsz=14.0
                   %s
           ''' % (min1,min2,max2,min3,max3,title,extra))
# 2.0
#  | window min2=1.0 max2=2.0 min3=4.0 max3=5.0

# you can either run scripts in ../p-cable-stack/ and/or comment/uncomment the line below
#SConscript('../p-cable-stack/SConstruct')

Fetch('stackfilled-filt3-packed.rsf','pcable',server)

Flow('stackfilled-filt3','stackfilled-filt3-packed',
                        '''
                        math output="input" |
                        window min3=4.0 max3=7.5 min1=0.21 max1=0.24 min2=1.0 |
                        put unit2=km unit3=km unit1=s label2=iline label3=xline label1=time
                        ''')

stackplot('stack','stackfilled-filt3','Stack',extra='screenratio=0.7 point1=0.7 point2=0.85 labelsz=12.0')

#sliceplot('stack-slice','stackfilled-filt3','Zero Offset')

Flow('vel','stackfilled-filt3','math output="1.5"')

Flow('mig3','stackfilled-filt3 vel','linmig3 adj=y doomp=y vel=${SOURCES[1]} antialias=t apt=40')

#stackplot('mig3','mig3','3D Kirchhoff')

sliceplot('mig3-slice','mig3','3D Kirchhoff')

### PWD

Flow('dip00','stackfilled-filt3','math output=0.0')

Flow('dip0','dip00 dip00','cat axis=4 ${SOURCES[1]}')

Flow('pwd','stackfilled-filt3 dip0',
        '''
        pwd dip=${SOURCES[1]}
        ''')

Flow('xpwd-notmig','pwd',
        '''
        window n4=1
        ''')

Flow('ypwd-notmig','pwd',
        '''
        window n4=1 f4=1
        ''')

Flow('xpwd','xpwd-notmig vel','linmig3 adj=y doomp=y vel=${SOURCES[1]} antialias=t apt=40')

Flow('ypwd','ypwd-notmig vel','linmig3 adj=y doomp=y vel=${SOURCES[1]} antialias=t apt=40')

#stackplot('xpwd','xpwd','XPWD')

#stackplot('ypwd','ypwd','YPWD')

Flow('xspr-notmig','stackfilled-filt3 dip0 stackfilled-filt3','pwspray2 dip=${SOURCES[1]} ns2=15 ns3=1 | stack | math K=${SOURCES[2]} output="K-input"')

Flow('yspr-notmig','stackfilled-filt3 dip0 stackfilled-filt3','pwspray2 dip=${SOURCES[1]} ns2=1 ns3=15 | stack | math K=${SOURCES[2]} output="K-input"')

stackplotc('xspr-notmig','xspr-notmig','XPWD',frame3=6,extra='screenratio=0.7 point1=0.7 point2=0.85 labelsz=12.0')

#stackplot('yspr-notmig','yspr-notmig','YPWD')

Flow('xspr','xspr-notmig vel','linmig3 adj=y doomp=y vel=${SOURCES[1]} antialias=t apt=40')

Flow('yspr','yspr-notmig vel','linmig3 adj=y doomp=y vel=${SOURCES[1]} antialias=t apt=40')

#stackplot('xspr','xspr','XPWD')

#stackplot('yspr','yspr','YPWD')

### Removing footprints

clipm=0.1

# crossline wavenumber

xlim = 79.9999

num2 = 1125

step2 = 0.142222

# inline wavenumber

start1 = -80.0001

num1 = 600

step1 = 0.266667

# time dimension

num3 = 16

step3 = 0.002

start3 = 0.002

###DOES NOT WORK ON REFLECTIONS 
###THERE MIGHT BE A CONSTANT COMPONENT

Flow('maskart1',None,
        '''
        spike n1=%g d1=%g o1=%g |
        math output="exp(-0.5*(x1*x1)) + exp(-0.05*((%g - x1)*(%g - x1))) + exp(-0.05*((%g - x1)*(%g - x1)))" |
        clip clip=%g |
        smooth rect1=10 |
        math output="1 - input/%g"
        '''%(num2,step2,-1.0*xlim,xlim,xlim,-1.0*xlim,-1.0*xlim,clipm,clipm))

Flow('maskart','maskart1','rtoc | spray axis=1 n=%g d=%g o=%g | spray axis=2 n=%g d=%g o=%g'%(num3,step3,start3,num1,step1,start1))

Flow('xpwd-prec','xpwd maskart',
        '''
        rtoc |
        fft3 axis=2 |
        fft3 axis=3 |
        math K=${SOURCES[1]} output="input*K" |
        fft3 axis=3 inv=y |
        fft3 axis=2 inv=y |
        real
        ''')

Flow('xspr-prec','xspr maskart',
        '''
        rtoc |
        fft3 axis=2 |
        fft3 axis=3 |
        math K=${SOURCES[1]} output="input*K" |
        fft3 axis=3 inv=y |
        fft3 axis=2 inv=y |
        real
        ''')

Flow('yspr-prec','yspr maskart',
        '''
        rtoc |
        fft3 axis=2 |
        fft3 axis=3 |
        math K=${SOURCES[1]} output="input*K" |
        fft3 axis=3 inv=y |
        fft3 axis=2 inv=y |
        real
        ''')

Flow('ypwd-prec','ypwd maskart',
        '''
        rtoc |
        fft3 axis=2 |
        fft3 axis=3 |
        math K=${SOURCES[1]} output="input*K" |
        fft3 axis=3 inv=y |
        fft3 axis=2 inv=y |
        real
        ''')

#stackplot('xpwd-prec','xpwd-prec','XPWD')

#stackplot('ypwd-prec','ypwd-prec','YPWD')

#stackplot('xspr-prec','xspr-prec','XPWD')

sliceplot('xspr-prec-slice','xspr-prec','XPWD')

#stackplot('yspr-prec','yspr-prec','YPWD')

# Generating a structure tensor in a slice by slice fashion 

Flow('delx','xpwd-prec','window n1=1 min1=%g'%sliceval)

Flow('dely','ypwd-prec','window n1=1 min1=%g'%sliceval)

Flow('delx-delx','delx','math output="input*input" | smooth rect1=10 rect2=10')
Flow('dely-dely','dely','math output="input*input" | smooth rect1=10 rect2=10')
Flow('delx-dely','delx dely','math K=${SOURCES[1]} output="K*input" | smooth rect1=10 rect2=10')
Flow('dely-delx','dely delx','math K=${SOURCES[1]} output="K*input" | smooth rect1=10 rect2=10')

Flow('eig10 ux0 uy0 eig20 vx0 vy0','delx-delx delx-dely dely-dely',
        '''
        pwdtensorh in2=${SOURCES[1]} in3=${SOURCES[2]} eps=0.0
        ux=${TARGETS[1]} uy=${TARGETS[2]} out2=${TARGETS[3]}
        vx=${TARGETS[4]} vy=${TARGETS[5]} normalize=n
        ''')

Flow('az0','vx0',
        '''
        math output="(180/%g)*asin(input)" |
        math output="input"'''%(math.pi))

# Spraying to the original dimensions

Flow('vx','vx0','spray axis=1 n=16 | put o1=0.21')

Flow('vy','vy0','spray axis=1 n=16 | put o1=0.21')

Flow('az00','az0','spray axis=1 n=16 | math output="input+90" | put o1=0.21')

#Flow('az','az00 vx vy','anisodiffuse vx=${SOURCES[1]} vy=${SOURCES[2]} eps=10 niter=5 repeat=1')

Flow('az','az00','smooth rect1=5 rect2=10 rect3=10')

#Result('az','window n1=1 min1=0.222 | grey color=j scalebar=y allpos=y')

Flow('xpwd-sm','xpwd-prec vx vy','anisodiffuse vx=${SOURCES[1]} vy=${SOURCES[2]} eps=10 niter=5 repeat=1 verb=n')

Flow('ypwd-sm','ypwd-prec vx vy','anisodiffuse vx=${SOURCES[1]} vy=${SOURCES[2]} eps=10 niter=5 repeat=1 verb=n')

# Linear combinations of sprays

Flow('azspr','xspr-prec yspr-prec az','azspr spry=${SOURCES[1]} az=${SOURCES[2]}')

#stackplot('azspr','azspr','AzSPR')

Flow('azpwd0','stackfilled-filt3 az dip0','azpwd az=${SOURCES[1]} dip=${SOURCES[2]}')

Flow('azpwd','stackfilled-filt3 az dip0 maskart vel',
        '''
        azpwd az=${SOURCES[1]} dip=${SOURCES[2]} |
        linmig3 adj=y doomp=y vel=${SOURCES[4]} antialias=t apt=40|
        rtoc |
        fft3 axis=2 |
        fft3 axis=3 |
        math K=${SOURCES[3]} output="input*K" |
        fft3 axis=3 inv=y |
        fft3 axis=2 inv=y |
        real
        ''')

#stackplot('azpwd','azpwd','AzPWD')

#sliceplot('azpwd-slice','azpwd','AzPWD')

Flow('azpwd-sm','azpwd vx vy','anisodiffuse vx=${SOURCES[1]} vy=${SOURCES[2]} eps=10 niter=5 repeat=1 verb=n')

#stackplot('azpwd-sm','azpwd-sm','AzPWD')
#stackplot('xpwd-sm','xpwd-sm','XPWD')
#stackplot('ypwd-sm','ypwd-sm','YPWD')

# Path-summation integral evalaution on crossline spray

Flow('linpi-data0','xspr-notmig','linpi v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n')

Flow('linpi-data','linpi-data0 maskart',
        '''
        rtoc |
        fft3 axis=2 |
        fft3 axis=3 |
        math K=${SOURCES[1]} output="input*K" |
        fft3 axis=3 inv=y |
        fft3 axis=2 inv=y |
        real
        ''')

#stackplot('linpi-data-stack','linpi-data','AzPWD')

sliceplot('linpi-data-slice','linpi-data','linpi-data',sliceval)

#halfsliceplot('linpi-data-hslice','linpi-data','linpi-data')

# by default antialiasing is t
"""
Flow('modl-5 snaps-5','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=n thr=100 mode=hard doanisodiff=n anisoeps=10 niter=5 repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=5 oniter=1 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        ''')

Flow('modl-2 snaps-2','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=n thr=100 mode=hard doanisodiff=n anisoeps=10 niter=5 repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=2 oniter=1 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        ''')

Flow('modl-3 snaps-3','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=n thr=100 mode=hard doanisodiff=n anisoeps=10 niter=5 repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=3 oniter=1 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        ''')

Flow('modl-10 snaps-10','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=n thr=100 mode=hard doanisodiff=n anisoeps=10 niter=5 repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=10 oniter=1 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        ''')
"""
#Flow('snaps-sm','snaps vx vy',
#       '''
#       window n4=1 f4=0 |
#       thr thr=0.015 mode=hard |
#       anisodiffuse vx=${SOURCES[1]} vy=${SOURCES[2]}
#       eps=10 niter=10 repeat=3 verb=n
#       ''')

# two sets of parameters to look at
# strong
thrs = 0.01
epss = 20
niters = 30

# weak
thrw = 0.008
epsw = 20
niterw = 30

"""
Flow('modl-str snaps-str','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=y thr=%g mode=hard doanisodiff=y anisoeps=%g niter=%d repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=2 oniter=10 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        '''%(thrs,epss,niters))

Flow('modl-wk snaps-wk','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=y thr=%g mode=hard doanisodiff=y anisoeps=%g niter=%d repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=2 oniter=10 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        '''%(thrw,epsw,niterw))

Flow('modl-str-100 snaps-str-100','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=y thr=%g mode=hard doanisodiff=y anisoeps=%g niter=%d repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=2 oniter=100 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        '''%(thrs,epss,niters))

Flow('modl-wk-100 snaps-wk-100','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=y thr=%g mode=hard doanisodiff=y anisoeps=%g niter=%d repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=2 oniter=100 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        '''%(thrw,epsw,niterw))

"""

Flow('modl-wk-100 snaps-wk-100','linpi-data dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=y thr=%g mode=hard doanisodiff=y anisoeps=%g niter=%d repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=y sm=y initer=2 oniter=20 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        '''%(thrw,epsw,niterw))

# removing anisodiffuse instabilities

filtering = '''
transp plane=12 |
bandpass fhi=50 |
transp plane=13 |
bandpass fhi=50 |
transp plane=31 |
transp plane=21
'''
# bandpass fhi=50 |

Flow('modl-wk-100ft','modl-wk-100',filtering)

#Flow('modl-str-100ft','modl-str-100',filtering)

sliceplot('modl-wk-100','modl-wk-100ft','Inversion Weak',sliceval,'pclip=97')

#sliceplot('modl-str-100','modl-str-100ft','Inversion Strong')

# looking at different iterations

Flow('0snaps-wk-100','snaps-wk-100 vx vy',
        '''
        window n4=1 |
        thr thr=%g mode=hard |
        anisodiffuse eps=%g niter=%d repeat=1
        vx=${SOURCES[1]} vy=${SOURCES[2]}
        '''%(thrw,epsw,niterw) + ' | ' + filtering)

Flow('00snaps-wk-100','snaps-wk-100 vx vy',
        '''
        window n4=1
        ''')

Flow('1snaps-wk-100','snaps-wk-100 vx vy',
        '''
        window n4=1 f4=1 |
        thr thr=%g mode=hard |
        anisodiffuse eps=%g niter=%d repeat=1
        vx=${SOURCES[1]} vy=${SOURCES[2]}
        '''%(thrw,epsw,niterw) + ' | ' + filtering)

Flow('10snaps-wk-100','snaps-wk-100 vx vy',
        '''
        window n4=1 f4=9 |
        thr thr=%g mode=hard |
        anisodiffuse eps=%g niter=%d repeat=1
        vx=${SOURCES[1]} vy=${SOURCES[2]}
        '''%(thrw,epsw,niterw) + ' | ' + filtering)
"""
Flow('100snaps-wk-100','snaps-wk-100 vx vy',
        '''
        window n4=1 f4=99 |
        thr thr=%g mode=hard |
        anisodiffuse eps=%g niter=%d repeat=1
        vx=${SOURCES[1]} vy=${SOURCES[2]}
        '''%(thrw,epsw,niterw) + ' | ' + filtering)
"""
Flow('100snaps-wk-100','snaps-wk-100 vx vy',
        '''
        window n4=1 f4=19 |
        thr thr=%g mode=hard |
        anisodiffuse eps=%g niter=%d repeat=1
        vx=${SOURCES[1]} vy=${SOURCES[2]}
        '''%(thrw,epsw,niterw) + ' | ' + filtering)

#stackplot('0snaps-wk-100','0snaps-wk-100','ZO',100,300,6)

#stackplot('1snaps-wk-100','1snaps-wk-100','ZO',100,300,6)

#stackplot('10snaps-wk-100','10snaps-wk-100','ZO',100,300,6)

#stackplot('100snaps-wk-100','100snaps-wk-100','ZO',100,300,6)

#stackplot('00snaps-wk-100','00snaps-wk-100','ZO',100,300,6)

# advantage of anisotropic smoothing

halfsliceplot('hs-modl-wk-100','modl-wk-100ft','Inversion Weak')

Flow('modl-wk-100ft-thr','snaps-wk-100','window n4=1 f4=19 | ' + filtering + ' | thr thr=%g mode=hard'%thrw)

halfsliceplot('hs-modl-wk-100-thr','modl-wk-100ft-thr','Inversion Weak')

# going back to data domain

Flow('modl-wk-100ft-zo','modl-wk-100ft dip0 az vx vy vel',
        '''
        piazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=n sm=y adj=n doomp=y
        ''')

#stackplot('modl-wk-100ft-zo','modl-wk-100ft-zo','ZO',100,300,6)

# restoring w#

# using our result as a starting model

Flow('d0','modl-wk-100ft-zo xspr-notmig','add scale=-1,1 ${SOURCES[1]}')

thrr = thrw+0.003
epsr = epsw
niterr = niterw

# path-summation is disabled
Flow('y y-snaps','d0 dip0 az vx vy vel',
        '''
        lspiazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dothr=n thr=0.0 mode=hard doanisodiff=n anisoeps=0.0 niter=1 repeat=1
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=n sm=y initer=2 oniter=1 doomp=y snaps=y dsnaps=1 snapsf=${TARGETS[1]}
        ''')

# we utilize snaps since regularization should be applied on the sum only

Flow('diffractivity','y-snaps modl-wk-100ft vx vy',
        '''
        add scale=1,1 ${SOURCES[1]} |
        thr thr=%g mode=hard |
        anisodiffuse vx=${SOURCES[2]} vy=${SOURCES[3]}
        eps=%g niter=%d repeat=1 verb=n
        '''%(thrr,epsr,niterr) + ' | ' + filtering)

#stackplot('diffractivity','diffractivity','ZO',100,300,6)

#stackplot('x0','modl-wk-100ft','ZO',100,300,6)

Flow('diffractions','diffractivity dip0 az vx vy vel',
        '''
        piazpwdmig3 v_1=1.3 v_2=1.45 v_3=1.55 v_4=1.75 fftDo=n apt=40
        dip=${SOURCES[1]} az=${SOURCES[2]} vx=${SOURCES[3]} vy=${SOURCES[4]} vel=${SOURCES[5]}
        domod=y dopi=n sm=y adj=n doomp=y
        ''')

#stackplotc('diffractions','diffractions','ZO',100,300,6)

#stackplot('azpwd0','azpwd0','ZO',100,300,6)

Flow('noise','xspr-notmig diffractions','add scale=1,-1 ${SOURCES[1]}')

#stackplotc('noise','noise','ZO',100,300,6)

Flow('noise-ortho diffractions-ortho','noise diffractions',
        '''
        ortho rect1=5 rect2=10 rect3=10 niter=30 sig=${SOURCES[1]} sig2=${TARGETS[1]}
        ''')

stackplotc('noise-ortho','noise-ortho','ZO',100,300,6,extra='screenratio=0.7 point1=0.7 point2=0.85 labelsz=12.0')

stackplot('diffractions-ortho','diffractions-ortho','ZO',100,300,6,extra='screenratio=0.7 point1=0.7 point2=0.85 labelsz=12.0')

# plots for the 4th reviewer: asked for ``some'' frequency plots
def FK(pad=1000):
        return '''
    window n2=1 f2=100 |
        costaper nw1=5 |
    pad end1=%d |
        bandpass flo=40.0 fhi=200 |
    fft1 | fft3 axis=2 |
    window min2=-40 max2=40 |
    math output="abs(input)" | real
    '''%pad

def F(pad=1000):
        return '''
    window n2=1 f2=100 |
        costaper nw1=5 |
    pad end1=%d |
        bandpass flo=40.0 fhi=200 |
    fft1 | stack | scale axis=2 |
    math output="abs(input)" | real
    '''%pad

for d in [('stackfilled-filt3','Stack (refl + diffr + noise))'),('xspr-notmig','Stack (diffr + noise)'),('diffractions-ortho','Diffr (from inversion)')]:
        nm,ti = d
        Flow(nm+'-fk',nm,FK())
        Flow(nm+'-fft',nm,F())
    # screenratio = 0.4 vpconvert format=png pen=gd fat=3 serifs=n bgcolor=w xshift=-2.4 n1=1570 ./Fig/sf-f3-fft.vpl
        # vpconvert format=png pen=gd fat=3 serifs=n bgcolor=w ./Fig/*-fk.vpl 
    # vpconvert format=png pen=gd fat=3 serifs=n bgcolor=w ./Fig/*-fft.vpl
        Result(nm+'-fk','grey title="%s"'%ti)
        Result(nm+'-fft','graph label2="Abs A" unit2=  title="%s"'%ti)

Flow('sf-f3-nw-fft','stackfilled-filt3-packed',FK(pad=0))
Result('sf-f3-nw-fft','grey title="Stack (refl + diffr + noise) no win"')

End()

sfmath
sfwindow
sfput
sftransp
sfbyte
sfgrey3
sflinmig3
sfgrey
sfcat
sfpwd
sfpwspray2
sfstack
sfspike
sfclip
sfsmooth
sfrtoc
sfspray
sffft3
sfreal
sfpwdtensorh
sfanisodiffuse
sfazspr
sfazpwd
sflinpi
sflspiazpwdmig3
sfbandpass
sfthr
sfpiazpwdmig3
sfadd
sfortho
sfcostaper
sfpad
sffft1
sfscale
sfgraph