up [pdf]
from rsf.proj import *
#from rsf.recipes.tpx import FPX
import math

# define FPX so that we can use transpose with larger memsize
def FPX(fpx,data,
        np,               # number of slopes
        nw,               # number of frequencies
        p0=-1,            # first slope
        dp=None,          # slope increment
        v0=0,             # velocity continuation
        m=1000            # transpose memsize
        ):

    if not dp:
        dp=-2.0*p0/(np-1)

    # TX -> FX
    fx = 'fx-'+data
    if (v0 > 0):
        Flow(fx,data,
             '''
             fft1 | window n1=%d | fft3 axis=2 | 
             vczo2 v0=0 nv=1 dv=%g | 
             window | fft3 axis=2 inv=y
             ''' % (nw,v0))
    else:
        Flow(fx,data,'fft1 | window n1=%d' % nw)

    # FX -> XPF
    xpf = 'xpf-'+data
    basis = 'basis-'+data
    Flow([xpf,basis],fx,
         '''
         transp memsize=%i|
         cltft basis=${TARGETS[1]} dip=y 
         p0=%g dp=%g np=%d 
         rect=3 niter=1000 verb=n
         ''' % (m,p0,dp,np),split=[1,nw],
         reduce='cat axis=3')

    Flow(fpx,[xpf,basis],
         'mul ${SOURCES[1]} | transp plane=13 memsize=%i'%m,
         split=[2,np])

def TPX(tpx,data,
        nt,               # number of time samples
        np,               # number of slopes
        nw=0,             # number of frequencies
        p0=-1,            # first slope
        dp=None,          # slope increment
        m=1000            # transpose memsize
        ):

    fpx = 'fpx-'+data

    nt2=nt
    if nt2%2:
        nt2 += 1
    nw0=nt2/2+1
    if not nw:
        nw = nw0

    FPX(fpx,data,np,nw,p0,dp,m)

    Flow(tpx,fpx,
         '''
         pad n1=%d | fft1 inv=y
         ''' % nw0,split=[3,'omp'])



##################################
# Usefule plotting functions
##################################

def velplot(title,label1='Depth',unit1='km',max1=2):
    return '''
    window max1=%g |
    grey color=j allpos=y title="%s" scalebar=y
    barlabel=Velocity barunit=km/s
    label1="%s" unit1="%s" label2=Lateral unit2=km
    barreverse=y pclip=100 
    ''' % (max1,title,label1,unit1)

def graph(transp,o2,d2,n2,col,fat,extra=''):
    return '''
    window max1=2 |
    graph transp=%d yreverse=y pad=n min2=%g max2=%g
    plotcol=%d plotfat=%d wantaxis=n wanttitle=n %s
    ''' % (transp,o2,o2+(n2-1)*d2,col,fat,extra)

def graphw(min1,max1,transp,o2,d2,n2,col,fat,extra=''):
    return '''
    window min1=%g max1=%g |
    graph transp=%d yreverse=y pad=n min2=%g max2=%g
    plotcol=%d plotfat=%d wantaxis=n wanttitle=n %s
    ''' % (min1,max1,transp,o2,o2+(n2-1)*d2,col,fat,extra)

def wiggle(title):
    return  '''
    window max1=2 j2=2 |
    wiggle transp=y yreverse=y poly=y 
    title="%s" label2=Offset unit2=km zplot=0.5
    wherexlabel=t wheretitle=b labelsz=12 titlesz=15
    ''' % title

Flow('modl',None,
     '''
     spike n1=1501 o1=-10 d1=0.02 n2=4
     nsp=4 k2=1,2,3,4 mag=0.4,0.7,1,1.5
     ''')
Flow('refl',None,
     '''
     spike n1=1501 n2=4 nsp=4 k2=1,2,3,4
     mag=0.0909091,0.1428570,0.1111110,0.2000000
     ''')
Flow('mod1','modl','window min1=0')

################################################################
#1.25 too shallow
#2.0 works perfectly
#1.75 - little little problems on the last offset
#Flow('rmodl','modl',
#     '''
#     pad n2=100 | noise rep=y seed=112012
#     type=n mean=1.25 range=1
#     ''')
Flow('rmodl','modl',
     '''
     pad n2=100 | noise rep=y seed=112012
     type=n mean=2.0 range=1
     ''')
################################################################

Flow('amodl','modl rmodl','cat axis=2 ${SOURCES[1]}')

Flow('rrefl','modl',
     '''
     pad n2=100 | noise rep=y seed=122012 |
     math output="input^3"
     ''')
Flow('prefl','rrefl','clip2 lower=10')
Flow('mrefl','rrefl','clip2 upper=-10')
################################################################
#Flow('drefl','prefl mrefl','add scale=0.01,0.01 ${SOURCES[1]}')
Flow('drefl','prefl mrefl','add scale=0.01,0.01 ${SOURCES[1]}')
################################################################

Flow('arefl','refl drefl','cat axis=2 ${SOURCES[1]}')

Flow('unif','mod1','unif2 n1=101 d1=0.02 v00=5,6,8,10,15')
#Result('modl','unif',velplot('Velocity Model',max1=2))

Flow('mod2','unif','math output=1.5+2*x1')

Plot('modl','mod2',velplot('Velocity Model',max1=2))
Plot('modla','mod1',graph(0,0,0.02,101,0,20,'scalebar=y'))
Plot('modlb','mod1',graph(0,0,0.02,101,7,5,'scalebar=y'))
Result('modl2','modl modla modlb','Overlay')

#Flow('data','amodl arefl',
#     '''
#     kirmod nt=600 dt=0.004 freq=25 refl=${SOURCES[1]}
#     nh=51  dh=0.05 h0=0  
#     ns=501 ds=0.02 s0=0 cmp=y
#     vel=1.5 gradz=1 type=v |
#     put label2=Offset unit2=km label3=Midpoint unit3=km
#     ''',split=[1,1501],reduce='add')
nx = 501
dx = 0.02
ox = 0
nt = 600
dt = 0.004
ot = 0
nh = 24
dh = 0.05
oh = 0


nz = 2001
dz = 0.002
zo = 0
# create representation of model
nrefl = 100
modllst = []
# make velocity model for conversion

for i in range(nrefl):
        Flow('diffr-model-ampl-%i'%i,'drefl',
             '''
             window n2=1 f2=%i |
             spray axis=1 n=%i d=%g o=%g
             '''%(i,nz,dz,zo))
        Flow('diffr-model-%i'%i,'rmodl diffr-model-ampl-%i'%i,
             '''
             window n2=1 f2=%i |
             unif2 n1=%i d1=%g o1=%g v00=2,1|
             ai2refl | 
             mul ${SOURCES[1]}
             '''%(i,nz,dz,zo))
        modllst.append('diffr-model-%i'%i)
Flow('diffr-model',modllst,'cat axis=3 d=1 o=0 ${SOURCES[1:%i]} | stack axis=3'%(len(modllst)))
Flow('velo-model-depth','diffr-model','math output="2+x1"')
Flow('diffr-model-t','diffr-model velo-model-depth',
     'depth2time dt=%g t0=%g nt=%i velocity=${SOURCES[1]}'%(dt,ot,nt))
Flow('diffr-response','diffr-model-t','ricker2 | smooth rect2=2| scale dscale=-1')

Flow('diffr-p','rmodl drefl',
     '''
     kirmod nt=600 dt=0.004 freq=25 refl=${SOURCES[1]}
     nh=24  dh=0.05 h0=0  
     ns=501 ds=0.02 s0=0 cmp=y
     vel=2.0 gradz=1 type=v |
     put label2=Offset unit2=km label3=Midpoint unit3=km |
     window | pow pow1=1 
     ''',split=[1,1501],reduce='add')
# add bandpassed noise
# working well with range=.002 var=2e-6
# when mute added, working with range=.002 var=5e-6
#trying to increase
Flow('diffr','diffr-p',
   'noise range=.002 var=1e-5 rep=y | bandpass flo=10 fhi=60 | add  ${SOURCES[0]}')
#vel=1.5 type=c |
#gradz=1 and type=v
Result('diffr','transp plane=23 | grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 wanttitle=n')
Result('diffr0-n','diffr','window n2=1 min2=0 | grey gainpanel=all min1=0.8 max1=1.6 min2=1.0 max2=9.0 title="Noisy Diffraction Data"')
pmin1 = 0.8
pmax1 = 1.6
pmin2 = 1.0
pmax2 = 9.0
point1 = .8
point2 = .7
Result('diffr3-n','diffr',
   '''
   transp plane=23 |
   window min1=%g max1=%g
   min2=%g max2=%g |
   byte gainpanel=a |
   grey3 flat=n
   title="Noisy Diffraction Data"
   frame1=%i frame2=%i frame3=%i
   point1=%g point2=%g
   '''%(pmin1,pmax1,pmin2,pmax2,
       .5*(pmax1-pmin1)/dt,(3.62-pmin2)/dx,0,
       point1,point2))

# Making offsets to be half offsets
# creating CMP_Y additional dimension for
# sfpreconstkirchhoff

# Padding is extremely important
# Looks like if you have edge effects slope decomp
# might be wrong

#| pad end1=100
Flow('tiffr','diffr','costaper nw1=400 nw3=40')

#Result('diffr','transp plane=23 | grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 wanttitle=n')
#Result('tiffr','transp plane=23 | grey gainpanel=all wanttitle=n')

# dh=0.05 => 0.025
Flow('data','tiffr',
            '''
            spray axis=4 n=1 |
            put label4=CMP_Y |
            transp plane=23 |
            transp plane=34 |
            put d4=0.025 
            ''')

# VC range

v0=2.0#1.5
nv=151
dv=0.01

# first step
# 4th dimension is offset
# 3rd dimension is CMP_Y n3=1

Flow('mig_no_halfint','data','preconstkirch aal=y zero=n vel=%g' %v0, split=[4,24])

# correct phase

Flow('mig','mig_no_halfint','transp plane=24 | transp plane=34 | halfint inv=y adj=y | transp plane=23 | transp plane=34')
#Result('mig','window | grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 wanttitle=n')
#Result('mig_nw','mig','window | grey gainpanel=all wanttitle=n')

# slope decomposition

pad=1000
padft=501#0
nw=231#501

# get rid of CMP_Y by simple window command
Flow('warp','mig','window | t2warp pad=%i'%(pad))
Flow('diffr-response-warp','diffr-response','window | t2warp pad=%i'%(pad))
#Flow('warp','mig','window n4=1 | t2warp pad=%i'%(pad))

# Padding is extremely important
# Looks like if you have edge effects slope decomp
# might be wrong
# so it is a really good idea to check the spectrum
# a really good way to plot spectrum
#<warp.rsf sffft1 | sfmath output="abs(input)" | sfreal | sfgraph | sfpen

# I slope decompose all the offsets at once
np=351
#<mig.rsf sfwindow min1=0.8 max1=1.2 min2=3.5 max2=5.0 f4=-1 | sfgrey | sfpen
p0=-1.7
dp = -p0/((np-1)/2)
FPX('fpx','warp',np=np,p0=p0,nw=nw,v0=0.0,m=10000)
#FPX('fpx','warp',np=np,p0=p0,nw=501,v0=0.0)
FPX('diffr-response-fpx','diffr-response-warp',np=np,p0=p0,nw=nw,v0=0.0,m=10000)
Flow('diffr-response-txp','diffr-response-fpx','pad n1=%i|fft1 inv=y | t2warp inv=y'%padft)
# check slope decomposition
offset_num = 0#20
Flow('txp','fpx','pad n1=%i | fft1 inv=y | t2warp inv=y | transp plane=23 memsize=1000'%padft)
#Result('txp','window n4=1 f4=%d | grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 title="Slope Decomposed Data offset#=%d" '%(offset_num,offset_num))
#Result('migcomp','txp','stack axis=3 | grey title="Stacked Slope Decomposed Data"')

# PSOVC

# if you use psovcp instead of psovc you can avert transp highlighted by ^^^

                                                                     #^^^
#Flow('f_hall_pk','fpx','fft3 axis=3 | transp plane=24 memsize=1000 | transp plane=34 memsize=1000')
Flow('f_hall_kp','fpx','fft3 axis=3 | transp plane=24 memsize=30000')
#Flow('f_hs_pk','fpx','sfwindow n4=1 f4=%d | fft3 axis=3 | transp plane=24 memsize=1000 | transp plane=34 memsize=1000'%(offset_num))

#sfwindow n4=1 min4=0.5
#Flow('f_hs_kp','fpx','sfwindow n4=1 f4=%d | fft3 axis=3 | transp plane=24 memsize=1000'%(offset_num))

#####################################################################################################
# checking the amplitudes
#<fpx.rsf sfstack | sffft1 inv=y | sfattr
#<fpx.rsf sfwindow n4=1 | sfstack | sffft1 inv=y | sfattr
# amplitudes have one order
# l2 norm is higher for full offset - not sure if it matters
# l2 for zo = 0.004
# l2 for all offsets = 0.018
# l2 is dependent on the # of coefficients
#####################################################################################################

# parallelezation is done for k - lateral wavenumber
#Flow('vc_f_hall_pk_psovc','f_hall_pk','psovc nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[4,1024])
Flow('vc_f_hall_kp_psovc','f_hall_kp','psovcp nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[4,np])
#Flow('vc_f_hs_pk_psovc','f_hs_pk','psovc nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[4,1024])
#Flow('vc_f_hs_kp_psovc','f_hs_kp','psovcp nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[4,np])

#Result('vc_psovc','grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 wanttitle=n')
#Result('vc_psovc15','vc_psovc','window n3=1 min3=1.5 | grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 wanttitle=n')
#Result('vc_psovc_hs','grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 wanttitle=n')
#Result('vc_psovc_hs15','vc_psovc_hs','window n3=1 min3=1.5 | grey gainpanel=all min1=0.5 max1=2.0 min2=1.0 max2=9.0 wanttitle=n')

### semblance estimation #####################################################
#^^^
#Flow('fvkp_hs','vc_f_hs_pk_psovc',
#     '''
#     transp memsize=50000 plane=34
#     ''')

# export OMP_NUM_THREADS=1 and memsize=100 for fft1 or you will get memory allocation problem and overhead 
# not sure how to make it optimal for parallel
#Flow('tvxp_hs','fvkp_hs','fft3 axis=3 inv=y | pad n1=%i | fft1 memsize=100 inv=y | t2warp inv=y'%padft,split=[4,np])
#Flow('tvxp_hs','fvkp_hs','fft3 axis=3 inv=y | pad n1=%i | fft1 memsize=100 inv=y | t2warp inv=y')
Flow('tvxp','vc_f_hall_kp_psovc',
   '''
   fft3 axis=3 inv=y | 
   pad n1=%i | 
   fft1 memsize=100 inv=y | 
   t2warp inv=y
   '''%padft,split=[4,np])
# the operation below takes a very long time (~6.5 hrs)
# doesnt seem to make too much of a difference either
#Flow('tpxv-taper','tvxp','transp plane=24 memsize=10000 | costaper nw2=60')
#Flow('tvx2_hs','tvxp_hs','mul $SOURCE | stack axis=4 norm=n')
#Flow('tvx2','tpxv-taper','mul $SOURCE | stack axis=2 | transp plane=23 memsize=10000')
# normalize so max is 1
Flow('tvxp-n','tvxp','norm apply=${SOURCE} ')
Flow('tvx2','tvxp-n','mul $SOURCE | stack axis=4 ')

# Flow('tvx_hs','tvxp_hs','stack axis=4 norm=n')
#Flow('tvx','tpxv-taper','stack axis=2 | transp plane=23 memsize=10000')

Flow('tvx','tvxp-n','stack axis=4 ')

# rect3 since x is the third dimension
# rect1 since t is the first dimension
#Flow('semb_hs','tvx_hs tvx2_hs',
#     '''
#     mul ${SOURCES[0]} |
#     divn den=${SOURCES[1]} niter=50 rect3=1000 rect1=5 |
#     clip2 lower=0
#     ''')


# new stuff by luke for probibilistic diffraction imaging.

drect1 = 5 # t 10 was a little much
drect2 = 3 # v 5
drect3 = 3 # x
#drect1 = 2 # t
#drect2 = 2 # v
#drect3 = 2 # x

Flow('mask','tvx',
   '''
   window n3=1 | 
   math output=1 | 
   mutter x0=2.0 inner=y v0=0.5 half=n t0=0.5  |
   mutter x0=2.2 inner=n v0=0.7 half=n t0=-0.5 |
   smooth rect1=5 rect2=5 |
   spray axis=3 n=501 d=.02 o=0
   ''')
Flow('semb-pre','tvx tvx2 mask',
     '''
     mul ${SOURCES[0]} |
     divnp den=${SOURCES[1]} rect1=%i rect2=%i rect3=%i |
     add mode=p ${SOURCES[2]} | 
     clip2 lower=0
     '''%(drect1,drect2,drect3))
#Flow('semb-pre','tvx tvx2',
#     '''
#     mul ${SOURCES[0]} |
#     divnp den=${SOURCES[1]} rect1=%i rect2=%i rect3=%i |
#     clip2 lower=0
#     '''%(drect1,drect2,drect3))
# normalize semblance
Flow('semb-a','semb-pre','stack axis=2 ')
Flow('semb-agc','semb-a','agc rect1=50 rect2=200')
Flow('semb-fact','semb-agc semb-a','divn den=${SOURCES[1]} rect1=%i rect2=%i'%(drect1,drect3))

Flow('semb-c','semb-pre',
   "softclip upper=`<${SOURCE} $RSFROOT/bin/sfscale dscale=3| $RSFROOT/bin/sfattr want=std | awk '{print $4}'`")

Flow('semb','semb-c',
   """
   clip2 lower=0  
   """)

#Flow('semb','semb-fact semb-pre','spray axis=2 n=%i o=%g d=%g |add mode=p ${SOURCES[1]} | costaper nw2=25'%(nv,v0+dv,dv))
# unsued muting (for now)
#     mutter x0=2.0 inner=y v0=0.5 half=n t0=0.25  |
#     mutter x0=2.2 inner=n v0=0.7 half=n t0=-0.25 |
Flow('vtrue','diffr',
   'window n2=1 | math output="2.0*sqrt((exp(x1)-1)/x1)"')
Flow('slice_true','tvx vtrue',
     'slice pick=${SOURCES[1]}')
Flow('vtrue-spray','vtrue','spray axis=3 n=%i d=%g o=%g'%(np,dp,p0))
Flow('slice-tpx','tvxp-n vtrue-spray','slice pick=${SOURCES[1]} | transp plane=23 memsize=10000')
vrect1 = 25   # in t
vrect2 = 200  # in x
# stack semb
Flow('semb-stk','semb-pre','stack axis=2')
# determine expectation value of velocity
# might have to do this assuming v of z media like in ovc paper
Flow('v-exp','semb-pre semb-stk',
   '''
   math output="input*x2" | 
   stack axis=2 | 
   divn den=${SOURCES[1]} rect1=%i rect2=%i | 
   clip2 lower=%g
   '''%(vrect1,vrect2,v0+dv))
# determine the variance
Flow('v-exp-var','v-exp semb-pre semb-stk',
   '''
   spray axis=2 n=%i d=%g o=%g |
   math A=${SOURCES[1]} output="A*(input-x2)^2" |
   stack axis=2 |
   divn den=${SOURCES[2]} rect1=%i rect2=%i
   '''%(nv,dv,v0+dv,vrect1,vrect2))
Flow('v-var-spray','v-exp-var','spray axis=2 n=%i d=%g o=%g'%(nv,dv,v0+dv))
Flow('d-vel','v-exp v-var-spray',
   '''
   spray axis=2 n=%i d=%g o=%g |
   math output="(input-x2)^2" |
   divn den=${SOURCES[1]} rect1=%i rect2=%i rect3=%i |
   scale dscale=.5 | 
   math output="exp(-1*input)" 
   '''%(nv,dv,v0+dv,drect1,drect2,drect3))
Flow('d-vel-stk','d-vel','math output="abs(input)" | stack axis=2')

Flow('bay-semb','semb v-exp','slice pick=${SOURCES[1]} ')

Flow('semb-dz','semb-pre',' deriv | add mode=abs')

Flow('semb-dx','semb-pre',
   '''
   transp plane=13 memsize=10000 | 
   deriv | 
   add mode=abs | 
   transp plane=13 memsize=10000
   ''')

Flow('semb-pre-eps','semb-pre',
    '''math output="input+`<${SOURCE} $RSFROOT/bin/sfattr want=std | awk '{print $4}'`" ''')
Flow('dsemb','semb-dz semb-dx semb-pre-eps',
   '''
   math A=${SOURCES[0]} B=${SOURCES[1]} output="sqrt(A*A+B*B)" |
   divnp den=${SOURCES[2]} rect1=%i rect2=%i rect3=%i |
   clip2 lower=0
   '''%(drect1,drect2,drect3))
#Flow('dsemb','semb-dz semb-dx ',
#   'math A=${SOURCES[0]} B=${SOURCES[1]} output="sqrt(A*A+B*B)"')
Flow('dsemb-stk','dsemb','add abs=$SOURCE | stack axis=2')



# normalize over 2nd axis
def normalize(file,nv,vo,dv,r1,r2,r3):
    Flow(file+'-d',file,
       'math output="abs(input)" | stack axis=2 | spray axis=2 n=%i d=%g o=%g'%(nv,dv,vo))
    Flow(file+'-n',[file,file+'-d'],'divnp den=${SOURCES[1]} rect1=%i rect2=%i rect3=%i'%(r1,r2,r3))

normalize('d-vel',nv,v0,dv,drect1,drect2,drect3)
# combined weights
Flow('wts','semb-pre d-vel dsemb d-vel-d',
   '''
   add mode=p ${SOURCES[1]} |
   add mode=p ${SOURCES[2]} |
   divnp den=${SOURCES[3]} rect1=%i rect2=%i rect3=%i |
   clip2 lower=0 
   '''%(drect1,drect2,drect3))

# make wieghted gathers
Flow('wtd-gath','wts tvxp-n',
   '''
   spray axis=4 n=%i d=%g o=%g |
   add mode=p ${SOURCES[1]} |
   stack axis=2|
   transp plane=23 memsize=10000
   '''%(np,dp,p0))
Flow('const-gath','tvxp-n','stack axis=2 | transp plane=23 memsize=10000')
# and the ideal gather
Flow('ideal-gath','diffr-response','scale dscale=-1| spray axis=2 n=%i d=%g o=%g'%(np,dp,p0))

Flow('wtd-img','wts tvx','add mode=p ${SOURCES[1]}')
Flow('prob-dimage','wtd-img','stack axis=2 | bandpass flo=10 ')

Flow('prob-dimage-squared','prob-dimage','add mode=p ${SOURCE}')

# multiply weights with images
#Flow('prob-dimage-numer','tvx dsemb semb d-vel',
#   '''
#   add mode=p ${SOURCES[1]} |
#   add mode=p ${SOURCES[2]} |
#   add mode=p ${SOURCES[3]} |
#   stack axis=2             
#   ''')


# create probibilistic diffraction image variance
Flow('prob-dimage-var','prob-dimage wtd-img prob-dimage-squared',
   '''
   spray axis=2 n=%i d=%g o=%g |
   math B=${SOURCES[1]} output="(input-B)^2" |
   stack axis=2  |
   divn den=${SOURCES[2]} rect1=%i rect2=%i |
   clip2 lower=0
   '''%(nv,dv,v0+dv,drect1,drect3))
#Flow('prob-dimage-var','dimage-numerator-var dimage-var-denom',
#   '''
#   divn den=${SOURCES[1]} rect1=%i rect2=%i |
#   sfmath output="abs(input)"
#   '''%(drect1,drect3))
#Flow('wtd-img','semb d-vel dsemb tvx d-vel-d',
#   '''
#   add mode=p ${SOURCES[1]}| 
#   add mode=p ${SOURCES[2]}|
#   add mode=p ${SOURCES[3]} |
#   divnp den=${SOURCES[4]} rect1=%i rect2=%i rect3=%i 
#   '''%(drect1,drect2,drect3))
# figures of the I(t,v,x) by midpoint and associated weights
wtlst = ['tvx','wts','wtd-img','semb-pre','d-vel-n','dsemb']
titles = ['I(t,v,x)','Combined Weights','Weighted Image','W1(t,v,x)','W2(t,v,x)','W3(t,v,x)']
gathers = ['wtd-gath','slice-tpx','const-gath','diffr-response-txp']
gathtitles = ['Probabilistic Weight Gather','Deterministic Gather','Equal Weight Gather','Ideal Gather']
#titles = ['i','ii','iii','iv','v']
colorlst = [' ','j',' ','j','j','j']
allpos = ['n','y','n','y','y','y']
clipss = [1e-4,.075,5e-5,.4,4,.1]

gmax1 = .8
gmin1 = 1.6
gmax2 = nv*dv+v0
gmin2 = v0+dv
pmin = -1
pmax = 1

for n in range(5):
    #x = (3,5,7)[n]
    x = (3.62,4.82,6.8,5.46,4.8)[n]#(3.75,4.9,6.8,5.5,4.8)
    plst = []
    plst3 = []
    gathlst = []
    for i in range(len(gathers)):
        gath = gathers[i]
        gtitle = gathtitles[i]
        Result('synth-n-'+gath+'-%i'%n,gath,
             '''
             window n3=1 min3=%g|
             grey min1=%g max1=%g title="%s"
             label1=Time unit1=s label2=Slope unit2="s\^2\_/km" 
             pclip=99
             '''%(x,gmax1,gmin1,gtitle))
    # get the velocity, graph
    Flow('v-exp-%i'%n,'v-exp','window n2=1 min2=%g '%x)
    Plot('v-exp-%i'%n,
       '''
       graph min1=%g max1=%g min2=%g max2=%g 
       label1= label2= unit1= unit2=  title= 
       transp=y n1tic=0 plotfat=20 plotcol=0
       '''%(gmin1,gmax1,gmin2,gmax2))
    # get variance
    Flow('v-exp-var-%i'%n,'v-exp-var','window n2=1 min2=%g'%x)
    # add +/- standard deviation to expectation velocity
    Flow('v-exp-low-%i'%n,['v-exp-var-%i'%n,'v-exp-%i'%n],
       'clip2 lower=0 | math A=${SOURCES[1]} output="A-sqrt(input)"')
    Plot('v-exp-low-%i'%n,
       '''
       graph min1=%g max1=%g min2=%g max2=%g 
       label1= label2= unit1= unit2=  title= 
       transp=y n1tic=0 dash=1 plotfat=20 plotcol=7
       '''%(gmin1,gmax1,gmin2,gmax2))
    Flow('v-exp-high-%i'%n,['v-exp-var-%i'%n,'v-exp-%i'%n],
       'clip2 lower=0 | math A=${SOURCES[1]} output="A+sqrt(input)"')
    Plot('v-exp-high-%i'%n,
       '''
       graph min1=%g max1=%g min2=%g max2=%g 
       label1= label2= unit1= unit2=  title= 
       transp=y n1tic=0 dash=1 plotfat=20 plotcol=7
       '''%(gmin1,gmax1,gmin2,gmax2))
#    Plot('lines-pre-%i'%n,['v-exp-low-%i'%n,'v-exp-high-%i'%n,'Overlay'])
#    Plot('lines-%i',['lines-pre-%i'%n,'v-exp-%i'%n],'Overlay')

    for k in range(len(wtlst)):
        item = wtlst[k]
        Flow(item+'-%i'%n,item,'window n3=1 min3=%g'%x)
        Plot(item+'-%i'%n,
           '''
           grey title="%s" at x=%g label2=Velocity
           unit2="km/s" min1=.8 max1=1.6
           color=%s allpos=%s  %g
           '''%(titles[k],x,colorlst[k],allpos[k],clipss[k]))
        point1d = 0.8
        point2d = 0.5
        Plot(item+'3-%i'%n,item,
           '''
           window min1=%g max1=%g  min3=%g max3=%g|
           byte gainpanel=2 allpos=%s |
           grey3 unit2="km/s" label2=Velocity 
           color=%s title="%s" flat=n
           frame1=%i frame2=%i frame3=%i 
           point1=%g point2=%g 
           allpos=%s
           screenht=15 screenratio=1.2
           titlesz=16 labelsz=8 
           n1tic=3 o1num=2.25 d1num=0.5
           '''%(pmin1,pmax1,pmin2,pmax2,
                allpos[k],colorlst[k],titles[k],
                .5*(pmax1-pmin1)/dt,nv/2-1,(x-pmin2)/dx,
                point1d,point2d,
                allpos[k]))
        Result(item+'3-%i'%n,item,
           '''
           window min1=%g max1=%g  min3=%g max3=%g|
           byte gainpanel=2 allpos=%s |
           grey3 unit2="km/s" label2=Velocity 
           color=%s title="%s" flat=n
           frame1=%i frame2=%i frame3=%i 
           point1=%g point2=%g 
           allpos=%s
           screenht=28 screenratio=2
           larnersz=85 titlesz=16
           '''%(pmin1,pmax1,pmin2,pmax2,
                allpos[k],colorlst[k],titles[k],
                .5*(pmax1-pmin1)/dt,nv/2,(x-pmin2)/dx,
                point1,point2,
                allpos[k]))
        plst.append(item+'-%i'%n)
        plst3.append(item+'3-%i'%n)
    Plot(plst[4]+'-o',[plst[4],'v-exp-%i'%n,'v-exp-low-%i'%n,'v-exp-high-%i'%n],'Overlay')
    Plot(plst[1]+'-o',[plst[1],'v-exp-%i'%n,'v-exp-low-%i'%n,'v-exp-high-%i'%n],'Overlay')
#    Plot(plst[4]+'-o',[plst[4],'lines-%i'%n],'Overlay')
    Result('synth-weights-n-%i-a'%n,[plst[0],plst[1],plst[2]],'SideBySideAniso')
    Result('synth-weights-n-%i-b'%n,[plst[3],plst[4]+'-o',plst[5]],'SideBySideAniso')
    Result('synth-weights-n3-%i-a'%n,[plst3[0],plst3[1],plst3[2]],'SideBySideIso')
    Result('synth-weights-n3-%i-b'%n,[plst3[3],plst3[4],plst3[5]],'SideBySideIso')
    Result('synth-weights-n3a-%i'%n,plst3,'TwoRows')
# and some image figures
pclip1=99.9
Plot('synth-prob-dimage','prob-dimage','grey min1=.8 max1=1.6 min2=1 max2=9 pclip=%g title="Probabilistic Weight Image"'%pclip1)
Plot('synth-prob-dimage-var','prob-dimage-var','grey pclip=99.995 min1=.8 max1=1.6 min2=1 max2=9 title="Image Variance" color=j scalebar=y allpos=y')
Plot('synth-top','synth-prob-dimage synth-prob-dimage-var','SideBySideAniso')
Plot('synth-pathint-img','tvx','stack axis=2 | grey pclip=%g min1=.8 max1=1.6 min2=1 max2=9 title="Equal Weight Image"'%pclip1)
Plot('synth-det-img','slice_true',' grey pclip=%g min1=.8 max1=1.6 min2=1 max2=9 title="True Velocity Image" '%pclip1)
Plot('synth-bot','synth-det-img synth-pathint-img','SideBySideAniso')
Result('synthss-n','synth-prob-dimage synth-prob-dimage-var synth-det-img synth-pathint-img','TwoRows')

Result('synth-n-prob-dimage','prob-dimage',
   'grey pclip=%g min1=.8 max1=1.6  min2=1 max2=9  title="Probabilistic Weight Image"'%pclip1)
Result('synth-n-prob-dimage-var','prob-dimage-var',
   'grey pclip=%g  title="Image Variance" color=j scalebar=y allpos=y'%pclip1)
Result('synth-n-top','synth-prob-dimage synth-prob-dimage-var','SideBySideIso')
Result('synth-n-pathint-img','tvx','stack axis=2 | grey pclip=%g min1=.8 max1=1.6 min2=1 max2=9 title="Equal Weight Image"'%pclip1)
Result('synth-n-det-img','slice_true',' grey pclip=%g min1=.8 max1=1.6 min2=1 max2=9 title="Deterministic Image" '%pclip1)
Result('synth-n-bot','synth-det-img synth-pathint-img','SideBySideIso')

Result('synth-n-diffr-response','diffr-response',
     'grey pclip=%g min1=.8 max1=1.6 label2=Midpoint unit2=km label1=Time unit1=s  min2=1 max2=9 title="Ideal Image"'%(pclip1))

End()

sfspike
sfwindow
sfpad
sfnoise
sfcat
sfmath
sfclip2
sfadd
sfunif2
sfgrey
sfgraph
sfspray
sfai2refl
sfmul
sfstack
sfdepth2time
sfricker2
sfsmooth
sfscale
sfkirmod
sfput
sfpow
sfbandpass
sftransp
sfbyte
sfgrey3
sfcostaper
sfpreconstkirch
sfhalfint
sft2warp
sffft1
sfcltft
sffft3
sfpsovcp
sfnorm
sfmutter
sfdivnp
sfagc
sfdivn
sfsoftclip
sfslice
sfderiv