from rsf.proj import *
import math

# define local FPX, mormalize
# 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))

# 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=5000            # 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=5000            # 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'])




### plotting functions


def section(title,label1='Time',unit1='s',min1=5.8,max1=8.0,extra=" "):
    return '''
    window min1=%g max1=%g |
    grey title="%s"
    label1="%s" unit1="%s" label2=Distance unit2=m %s
    ''' % (min1,max1,title,label1,unit1,extra)

def sectionw(title,label1='Time',unit1='s',min1=5.5,max1=7.0):
    return '''
    window min1=%g max1=%g max2=5500 |
    grey title="%s"
    label1="%s" unit1="%s" label2=Distance unit2=m
    ''' % (min1,max1,title,label1,unit1)

def plotvel(title,min1=5.5,max1=8.0):
    return '''
    grey min1=%g max1=%g minval=1480 maxval=2000 bias=1700 
    clip=200 allpos=n color=j scalebar=y barreverse=y barunit=m/s barlabel="Velocity"
    title="%s" label2="Distance" unit2=m 
    ''' %(min1,max1,title)

def plotdip(title,min1=5.5,max1=8.0):
    return '''
    grey min1=%g max1=%g color=j scalebar=y barlabel="Dip"
    title="%s" label2="Distance" unit2=m 
    ''' %(min1,max1,title)

def plotdipw(title,label1='Time',unit1='s',min1=5.5,max1=7.0):
    return '''
    window min1=%g max1=%g max2=5500 |
    grey title="%s" color=j scalebar=y barlabel="Dip"
    label1="%s" unit1="%s" label2=Distance unit2=m
    ''' % (min1,max1,title,label1,unit1)

frame1=(int)( (7.0-5.5)/0.002 )
frame2=(int)( (2000)/16.667 )
frame3=(int)( (1500-1400)/20 )

#print frame1
#print frame2
#print frame3

def plotdip3d(title,min1=5.5,max1=8.0,fr1=frame1,fr2=frame2,fr3=frame3):
    return '''
    window min1=%g max1=%g max2=5500 | byte | grey3 flat=n color=j scalebar=n barlabel="Dip"
    title="%s" label2="Distance" label3="Velocity" unit3="m/s" unit2=m frame1=%d frame2=%d frame3=%d
    point1=.6 point2=.7
    ''' %(min1,max1,title,fr1,fr2,fr3)

def plotpick(min2=1400,max2=2600):
    return '''
    graph yreverse=y transp=y plotcol=7 plotfat=7 
    pad=n min2=%g max2=%g wantaxis=n wanttitle=n
    '''%(min2,max2)

def plot3d(title,min1=5.5,max1=8.0,fr1=frame1,fr2=frame2,fr3=frame3):
    return '''
    window min1=%g max1=%g max2=5500 | byte | grey3 flat=n scalebar=n barlabel="Dip"
    title="%s" label2="Distance" label3="Velocity" unit3="m/s" unit2=m frame1=%d frame2=%d frame3=%d
    point1=.6 point2=.7
    ''' %(min1,max1,title,fr1,fr2,fr3)

####### Download data

Fetch('Nshots.su','nankai')

#sampling info
ot1 = 4
dt = .004
dx = 16.667
nx = 401
ox = 900

####### Convert to RSF

Flow('shots tshots','Nshots.su',
     'suread suxdr=y tfile=${TARGETS[1]}')

####### DC Removal

Flow('mean','shots','stack axis=1 | spray axis=1 n=5500 o=0.0 d1=0.002')

Flow('shotsdc','shots mean','add scale=1,-1 ${SOURCES[1]}')

####### Bandpass Filtering

Flow('shotsf','shotsdc','bandpass flo=10')

####### Windowing Data

Flow('shotsw','shotsf','window min1=5.5')

# windowing 5.5 and bandpass filtering significantly
# improves surface-consistent analysis

####### Mask zero traces

Flow('mask0','shotsw','mul $SOURCE | stack axis=1 | mask min=1e-20')

Flow('shots0','shotsw mask0','headerwindow mask=${SOURCES[1]}')

# update a database
Flow('tshots0','tshots mask0','headerwindow mask=${SOURCES[1]}')

####### Surface consistent

# Average trace amplitude
Flow('arms','shots0',
     'mul $SOURCE | stack axis=1 | math output="log(input)"')

# shot/offset indeces: fldr and tracf

Flow('indexshot','tshots0','window n1=1 f1=2')

Flow('offsets4index','tshots0',' headermath output=offset | dd type=float | window')

Flow('offsetindex','offsets4index','math output="abs(input) - 170" | dd type=int')

# receiver/midpoint

Flow('midpoint','tshots0','window n1=1 f1=5')

Flow('cmps4index','tshots0',' headermath output=cdp | dd type=float | math output="input*16.667" | window')

Flow('recv','cmps4index offsets4index','add scale=1,0.5 ${SOURCES[1]} | math output="input - 13799" | dd type=int')

Flow('index','indexshot offsetindex',
        '''
        cat axis=2 ${SOURCES[1]}
        ''')

Flow('extindex','index midpoint',
        '''
        cat axis=2 ${SOURCES[1]}
        ''')

Flow('extindrecv','extindex recv',
        '''
        cat axis=2 ${SOURCES[1]}
        ''')

def plot(title):
    return '''
    spray axis=1 n=1 | 
    intbin head=${SOURCES[1]} yk=fldr xk=tracf | window | 
    grey title="%s" label2="Shot Number" unit2= 
    label1="Offset Number" unit1= scalebar=y
    ''' % (title)

def plotb(title,bias=-5):
    return '''
    spray axis=1 n=1 | 
    intbin head=${SOURCES[1]} yk=fldr xk=tracf | window | 
    grey title="%s" label2="Shot Number" unit2= 
    label1="Offset Number" unit1= scalebar=y clip=3 bias=%g
    ''' % (title,bias)

# Display in shot/offset coordinates
Flow('varms','arms tshots0','spray axis=1 n=1 | intbin head=${SOURCES[1]} yk=fldr xk=tracf | window')
Plot('varms','arms tshots0',plotb('Log-Amplitude'))

#prog = Program('surface-consistent.c')
#sc = str(prog[0])

# recv index

# get model dimensions
#Flow('recvmodel',['arms','extindrecv',sc],
#     './${SOURCES[2]} index=${SOURCES[1]} verb=y')
Flow('recvmodel',['arms','extindrecv'],
     'surface-consistent index=${SOURCES[1]} verb=y')

# find a term
Flow('recvsc',['arms','extindrecv','recvmodel'],
     '''
     conjgrad surface-consistent index=${SOURCES[1]} 
     mod=${SOURCES[2]} niter=150
     ''')

# project to a data space
Flow('recvscarms',['recvsc','extindrecv'],
     'surface-consistent index=${SOURCES[1]} adj=n')

Plot('recvvscarms','recvscarms tshots0',
       plotb('Source, Offset, CDP, Recv S-C Log(A)'))

# compute difference
Flow('recvadiff','arms recvscarms','add scale=1,-1 ${SOURCES[1]}')

Plot('recvadiff','recvadiff tshots0',plot('s,h,cdp,r difference'))

### apply to traces to all times - no windowing is considered

Flow('ampl','recvscarms',
     'math output="exp(-input/2)" | spray axis=1 n=5500 d=0.002 o=0')

Flow('shotsf0','shotsf mask0','headerwindow mask=${SOURCES[1]}')

Flow('shots-preproc','shotsf0 ampl','mul ${SOURCES[1]}')

Plot('shots-preproc','shots-preproc','window n2=100 | grey min1=6.0 max1=8.0 title="Shots Preproc"')

Plot('shots-raw','shots0','window n2=100 | grey min1=6.0 max1=8.0 title="Shots Raw"')

### make stack figures to compare

# Resample to 4 ms

Flow('subsampled','shots-preproc',
     '''
     bandpass fhi=125 | window j1=2
     ''')

#Result('spectra','subsampled',
#       'spectra all=y | graph title="Subsampled Spectra"')

#Result('spectra-check','shots-preproc',
#       'spectra all=y | graph max1=160 title="Spectra Check"')

Flow('cmps mask','subsampled tshots0',
        '''
        intbin xk=tracf yk=cdp head=${SOURCES[1]} mask=${TARGETS[1]} |
        put o3=900 d3=16.667 label3=Distance unit3=m |
        pow pow1=2
        ''')

Flow('cmps-raw mask-raw','shots0 tshots0',
        '''
        intbin xk=tracf yk=cdp head=${SOURCES[1]} mask=${TARGETS[1]} |
        put o3=900 d3=16.667 label3=Distance unit3=m
        ''')

Flow('offset-file','tshots0',
     '''
     sfheadermath output=offset | dd type=float |
     intbin xk=tracf yk=cdp head=$SOURCE
     ''')

Flow('watervel','cmps','window n2=1 | math output="1500"')

Flow('nmo-wat-vel-preproc','cmps offset-file mask watervel',
     '''
     nmo half=n offset=${SOURCES[1]} mask=${SOURCES[2]}
     velocity=${SOURCES[3]}
     ''')

Flow('nmo-wat-vel-raw','cmps-raw offset-file mask watervel',
     '''
     nmo half=n offset=${SOURCES[1]} mask=${SOURCES[2]}
     velocity=${SOURCES[3]}
     ''')

Flow('stack-wat-vel-preproc','nmo-wat-vel-preproc','stack')

Flow('stack-wat-vel-raw','nmo-wat-vel-raw','stack')

#Result('stack-preproc','stack-wat-vel-preproc',section('Stack 1.5 km/s Preproc'))

#Result('stack-raw','stack-wat-vel-raw',section('Stack 1.5 km/s Raw'))

### DMO stacking 

nv=60

Flow('stacks','cmps offset-file mask',
     '''
     window min1=4 |
     stacks half=n v0=1400 nv=%g dv=20 
     offset=${SOURCES[1]} mask=${SOURCES[2]}
     '''%nv)#,split=[3,'omp'])

#>>>
# Consider stacking without time windowing to check
# if slope decomposition is affected
#Flow('nw-stacks','cmps offset-file mask',
#     '''
#     stacks half=n v0=1400 nv=%g dv=20 
#     offset=${SOURCES[1]} mask=${SOURCES[2]}
#     '''%nv,split=[3,'omp'])

### interface

p = 100

min1=5.51
datamax1=7.99
max1 = 7.6
mute = '''
        mutter v0=130 x0=1300 t0=4.0 half=n inner=n |
        mutter x0=1400 v0=20 t0=5.0 half=n inner=y | 
        mutter v0=2500 x0=1400 t0=5.8 half=n inner=n |
        mutter v0=500 x0=1400 t0=7.0 half=n inner=y 
        '''

pick = '''
        pick rect1=50 rect2=20 vel0=1480
        '''

mutec = '''
        mutter v0=130 x0=1300 t0=4.0 half=n inner=n |
        mutter x0=1400 v0=20 t0=5.0 half=n inner=y | 
        mutter v0=2500 x0=1400 t0=5.8 half=n inner=n |
        mutter v0=500 x0=1400 t0=7.0 half=n inner=y 
        '''

pickc = '''
        pick rect1=80 rect2=20 vel0=1400
        '''

# Compute envelope for picking (stack without DMO)
Flow('envelope-stacks','stacks','envelope | scale axis=2')#,split=[3,'omp'])

Flow('vpick-stacks','envelope-stacks',mutec + ' | ' + pickc)

# Take a slice
Flow('slice-stacks','stacks vpick-stacks','slice pick=${SOURCES[1]}')

#Result('vpick-stacks',plotvel('NMO velocity'))

#Result('stacks-check-w','stacks','window n2=1 min2=1499 | window max2=5500 | grey title="Stack 1500 m/s (extracted)" min1=5.5 max1=7.0')

#Result('stacks-check','stacks','window n2=1 min2=1499 | grey title="Stack 1500 m/s (extracted)" min1=4.0')

Flow('stackst','stacks','costaper nw3=20')

# Apply double Fourier transform (cosine transform)
# pad n3=601 | 
Flow('cosft','stackst','cosft sign1=1 sign3=1')

# Transpose f-v-k to v-f-k
Flow('transp','cosft','transp')#,split=[3,'omp'])

# Fowler DMO: mapping velocities
Flow('map','transp',
     '''
     math output="x1/sqrt(1+0.25*x3*x3*x1*x1/(x2*x2))" | 
     cut n2=1
     ''')

Flow('fowler','transp map','iwarp warp=${SOURCES[1]} | transp')#,
     #split=[3,'omp'])

# Inverse Fourier transform
# | window n3=401 # does not improve the result severely

Flow('dmo','fowler','cosft sign1=-1 sign3=-1')

#>>>
#Flow('nw-stackst','nw-stacks','costaper nw3=20')
#Flow('nw-cosft','nw-stackst','cosft sign1=1 sign3=1')
#Flow('nw-transp','nw-cosft','transp',split=[3,'omp'])
#Flow('nw-map','nw-transp',
#     '''
#     math output="x1/sqrt(1+0.25*x3*x3*x1*x1/(x2*x2))" | 
#     cut n2=1
#     ''')
#Flow('nw-fowler','nw-transp nw-map','iwarp warp=${SOURCES[1]} | transp',
#     split=[3,'omp'])
#Flow('nw-dmo','nw-fowler','cosft sign1=-1 sign3=-1')

# Looks like DMO helps to preserve diffractions 
# dipping flanks from right to left right part of a section

#Result('dmo-wat-vel','dmo','window n2=1 min2=1500 | ' + section('DMO 1500 m/s (extracted)'))

#Result('dmo-wat-vel-w','dmo','window n2=1 min2=1500 | ' + sectionw('DMO 1500 m/s (extracted)'))         

# Compute envelope for picking
Flow('envelope','dmo','envelope | scale axis=2')#,split=[3,'omp'])

Flow('vpick','envelope',mutec + ' | ' + pickc)

#Result('vpick',plotvel('DMO velocity'))

# Take a slice
Flow('slice','dmo vpick','slice pick=${SOURCES[1]}')

# 'grey title="Nankai DMO Stack" '
Result('slice',section('DMO stack'))

# 'window max2=5500 | grey title="DMO stack (extracted)" min1=5.5 max1=7.0'
#Result('slice-w','slice',sectionw('DMO Stacked Data'))

#>>>
#Flow('nw-envelope','nw-dmo','envelope | scale axis=2',split=[3,'omp'])
#Flow('nw-vpick','nw-envelope',mutec + ' | ' + pickc)
#Result('nw-vpick',plotvel('nw DMO velocity'))
#Flow('nw-slice','nw-dmo nw-vpick','slice pick=${SOURCES[1]}')
#Result('nw-slice',section('nw DMO stack'))

### PWD on 1.5 km/s image

Flow('vc15','slice',
        '''
        cosft sign2=1 |
        vczo v0=0.0 nv=1 dv=1500 |
        window |
        cosft sign2=-1
        ''')

#>>>
#Flow('nw-vc15','nw-slice',
#       '''
#       cosft sign2=1 |
#       vczo v0=0.0 nv=1 dv=1500 |
#       window |
#       cosft sign2=-1
#       ''')


#Result('vc15',section("VC 1.5 km/s"))

#Result('vc15-w','vc15',sectionw("VC 1.5 km/s"))

#Result('vc15-w','vc15','window max2=5500 | grey title="VC 1500 m/s" min1=5.5 max1=7.0')

diprect1=30
diprect2=20
diprect3=30

Flow('dip-vc','vc15','dip rect1=%d rect2=%d' % (diprect1, diprect2))

#Result('dip-vc',plotdip("VC 1.5 km/s Dip"))

#Result('dip-vc-w','dip-vc',plotdipw("VC 1.5 km/s Dip"))

# extract diffractions with pwd
#Flow('pwd-vc','vc15 dip-vc','pwd dip=${SOURCES[1]}')
# plane wave spray radius
ns = 2
# amplify reflections
Flow('refl-vc','vc15 dip-vc','pwspray dip=${SOURCES[1]} ns=%i | stack axis=2'%ns)
# and remove them
Flow('difr-vc','vc15 refl-vc','math A=${SOURCES[1]} output="input-A"')
#Result('pwd-vc-w','pwd-vc',sectionw("VC 1.5 km/s PWD"))

#Result('pwd-vc',section("VC 1.5 km/s PWD"))

#Flow('dif','pwd-vc',
Flow('dif','difr-vc',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')
# full wavefield
Flow('full','vc15',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')


#Result('inv-vc15-w','dif',sectionw("inv VC 1.5 km/s PWD"))

#Result('inv-vc15','dif',section("inv VC 1.5 km/s PWD"))

Result('dif',section('Separated Diffractions'))

#>>>
#Flow('nw-dip-vc','nw-vc15','dip rect1=%d rect2=%d' % (rect1, rect2))
#Result('nw-dip-vc',plotdip("nw VC 1.5 km/s Dip"))
#Result('nw-dip-vc-w','nw-dip-vc',plotdipw("nw VC 1.5 km/s Dip"))
#Flow('nw-pwd-vc','nw-vc15 nw-dip-vc','pwd dip=${SOURCES[1]}')
#Result('nw-pwd-vc-w','nw-pwd-vc',sectionw("nw VC 1.5 km/s PWD"))
#Result('nw-pwd-vc',section("nw VC 1.5 km/s PWD"))
#Flow('nw-dif','nw-pwd-vc',
#       '''
#       cosft sign2=1 |
#       vczo v0=1500.0 nv=1 dv=-1500 |
#       window |
#       cosft sign2=-1
#       ''')
#Result('nw-inv-vc15-w','nw-dif',sectionw("nw inv VC 1.5 km/s PWD"))
#Result('nw-inv-vc15','nw-dif',section("nw inv VC 1.5 km/s PWD"))

### Lets perform slope decomposition

# it is deep not sure if padding is needed
pad=1800#3000

Flow('warp','dif','t2warp pad=%i'%(pad))

Flow('warp-full','full','t2warp pad=%i'%(pad))

Plot('dif-spectrum','dif','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="Spectrum before warping"')

Plot('spectrum','warp','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="Spectrum after warping"')

#Result('warp','grey title="Warp"')

#>>>
#Flow('nw-warp','nw-dif','t2warp') 
#Plot('nw-dif-spectrum','nw-dif','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="nw Spectrum before warping"')
#Plot('nw-spectrum','nw-warp','sffft1 | sfmath output="abs(input)" | sfreal | stack | sfgraph title="nw Spectrum after warping"')

np=151
p0=-0.014
nw=901
dp = -2*p0/(np-1)

#?# -0.014 is it reasonable?
#?# do not understand why it is so small
#?# dt/dx = 0.5:0.002/1000:16.667 = tgalpha = 5
#?# is it a problem that we cut a dataset?
#?# Luke actually windows Barrolka
#?# I pad to 3000 and get 3000 instead of 3500
#?# may be should pad more than 3500 - 5000 and window frequencies
#?# try the above tests
#?# start with 0.005 may be - looks this is sufficient

#?# might make sence to perform v0 migration before slope decomposition
#?# it should reduce p0
FPX('fpx','warp',np=np,p0=p0,nw=nw,v0=0.0)

FPX('fpx-full','warp-full',np=np,p0=p0,nw=nw,v0=0.0)

Flow('txp','fpx','fft1 inv=y | t2warp inv=y | transp plane=23 memsize=1000')

Flow('txp-full','fpx-full','fft1 inv=y | t2warp inv=y | transp plane=23 memsize=1000')

#Result('txp',section("TXP"))

#Result('tx','txp','stack axis=3 | ' + section("TX"))

#>>> looks like slope decomposition does not care if we cut data or not
#FPX('nw-fpx','nw-warp',np=np,p0=p0,nw=1441,v0=0.0)
#Flow('nw-txp','nw-fpx','fft1 inv=y | t2warp inv=y | transp plane=23 memsize=1000')
#Result('nw-txp',section("nw TXP"))
#Result('nw-tx','nw-txp','stack axis=3 | ' + section("nw TX"))

### Perform OVC with 1.5 km/s and see if it makes sence on CIGs

Flow('fkp','fpx','transp plane=23 memsize=1000 | fft3 axis=2')

Flow('fkp-full','fpx-full','transp plane=23 memsize=1000 | fft3 axis=2')

#>>>
#Flow('nw-fkp','nw-fpx','transp plane=23 memsize=1000 | fft3 axis=2')

### trying to use psovc to avert additional transp ###
#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))

### trying to use psovc to avert additional transp ###
# creating offset dimension
#Flow('fhpk','fpx','spray axis=2 n=1 o=0.0 d=0.1 | fft3 axis=4')

#?# had some troubles to parallelize it
Flow('ovc-fkp-15','fkp','ovczo nv=%d dv=%g v0=%g' % (1,1500,0.0))

Flow('ovc-15','ovc-fkp-15','window | fft3 axis=2 inv=y | fft1 inv=y | t2warp inv=y | transp plane=23 memsize=5000')

#?# looks reasonable and even p0=-0.014 looks to be too wide
#Result('ovc-15-cig','ovc-15','window n3=1 min3=2516.7 | grey min1=6.0 max1=8.0 title="OVC 1.5 km/s CIG"')

Flow('ovc-15-image','ovc-15','stack')

#Result('ovc-15-image',section("OVC 1.5 km/s"))

#>>> we still have sharp slope change on the left
#Flow('nw-ovc-fkp-15','nw-fkp','ovczo nv=%d dv=%g v0=%g' % (1,1500,0.0))
#Flow('nw-ovc-15','nw-ovc-fkp-15','window | fft3 axis=2 inv=y | fft1 inv=y | t2warp inv=y | transp plane=23 memsize=5000')
#Result('nw-ovc-15-cig','nw-ovc-15','window n3=1 min3=2516.7 | grey min1=6.0 max1=8.0 title="nw OVC 1.5 km/s CIG"')

### Performing actual OVC

nv = 81 # was 60
dv = 5 # work with 10, was 20
v0 = 1400 # was 1400, then 1480
# minimum velocity allowed, used to look ok with v0+dv
vmin = 1485

Flow('ovc-1st-step','fkp','ovczop nv=%d dv=%g v0=%g | window' % (1,v0,0.0))
Flow('ovc-1st-step-full','fkp-full','ovczop nv=%d dv=%g v0=%g | window' % (1,v0,0.0))
#np
#Flow('ovc-fkp','ovc-1st-step','ovczo nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[3,'omp'],reduce='cat axis=4')

Flow('ovc-fkp','ovc-1st-step','ovczop nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[3,np],reduce='cat axis=4')

Flow('ovc-fkp-full','ovc-1st-step-full','ovczop nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[3,np],reduce='cat axis=4')

### trying to use psovc to avert additional transp - not sure if it is worth it ###
#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_hall_kp_psovc','f_hall_kp','psovcp nv=%d dv=%g v0=%g' % (nv,dv,v0),split=[4,np])

#?# f-v-k-p -> t-v-x-p, t-x-v-p - because fft3 seems to have a problem when applied to a different axis

Flow('txvp','ovc-fkp',
        '''
        transp plane=23 memsize=1000 |
        fft3 axis=2 inv=y |
        fft1 memsize=1000 inv=y |
        t2warp inv=y
        ''')

Flow('txvp-full','ovc-fkp-full',
        '''
        transp plane=23 memsize=1000 |
        fft3 axis=2 inv=y |
        fft1 memsize=1000 inv=y |
        t2warp inv=y
        ''')

# this is where fft3 fails
#Flow('tvxp','ovc-fkp','fft3 axis=3 inv=y | fft1 memsize=100 inv=y | t2warp inv=y',split=[4,np])

Flow('ovc','txvp','stack axis=4')

Flow('tvx','ovc','transp plane=23 memsize=5000')

#Flow('ovc','ovc-fkp','stack axis=4 | fft3 axis=3 inv=y | fft1 inv=y | t2warp inv=y')

#Result('ovc',section("OVC"))

#?# sharp boundary on CIG - zero on the left (from around +/-0.005) and non-zero on the right
#Result('ovc-cig','txvp','window n2=1 min2=2516.7 | transp plane=23 | grey min1=6.0 max1=8.0 title="OVC CIG"')

### Semblance calculation

Flow('txv1','txvp','stack norm=n axis=4')
Flow('txv2','txvp','mul $SOURCE | stack norm=n axis=4')

#Result('txvp',
#             '''
#             window n3=1 f3=1 | 
#             pow pow1=1 |
#             byte | 
#             grey3 flat=n frame1=250 frame2=250 frame3=100
#             ''')

#Plot('txv1','grey gainpanel=all title="Velocity Continuation" ',view=1)
#Plot('txv2','grey gainpanel=all title="What is txv2 check" ',view=1)

# optimal rect1=5 rect2=5
# rect1=10
Flow('sem','txv1 txv2',
     '''
     mul ${SOURCES[0]} | divn den=${SOURCES[1]} rect1=5 rect2=5 niter=25 |
     clip2 lower=0 
     ''')

# I dont know how to automate it

#Flow('max.txt','sem','attr want=max | awk "{print $2}" --out=$TARGET',stdout=0)

#Flow('max',None,'spike mag=`attr <sem.rsf want=max | awk "{print $4}"`')

# scaling actually improves picking
Flow('semb','sem','transp plane=23 | math output="input/36.4387"')

# Semblance picking parameters


# create mute
nx = 401
dx = 16.667
ox = 900

# determine expectation value of velocity
# smoothing for velocity
vrect1 = 25   #  in t
vrect2 = 25   #  in x

# smoothing for division
drect1 = 2 # t 
drect2 = 2 # v 
drect3 = 2 # x

# smoothing for dix inversion
dixrect1 = 5  # or 15
dixrect2 = 2  # or 2


# old mute
#Flow('mute-pre','semb',
#   '''
#   window n3=1 | math output=1 |
#   mutter x0=1400 v0=100 t0=5.0 half=n inner=y | 
#   mutter v0=2500 x0=1400 t0=5.8 half=n inner=n |
#   mutter v0=2000 x0=1500 t0=5.8  inner=n |
#   smooth rect1=20 rect2=5 |
#   spray axis=3 n=%i d=%g o=%g |
#   add mode=p ${SOURCE}
#   '''%(nx,dx,ox))

mutesemb = '''
        mutter v0=130 x0=1300 t0=4.0 half=n inner=n |
        mutter x0=1400 v0=20 t0=5.0 half=n inner=y | 
        mutter v0=2500 x0=1400 t0=5.8 half=n inner=n |
        mutter v0=500 x0=1400 t0=7.0 half=n inner=y 
        '''

Flow('mute-pre','semb',
   '''
   window n3=1 | math output=1 |
   '''+mutesemb+''' |
   smooth rect1=10 rect2=5 |
   spray axis=3 n=%i d=%g o=%g |
   add mode=p ${SOURCE}
   '''%(nx,dx,ox))

# use agc to normalize the semblance over time values
Flow('mute-a','mute-pre','stack axis=2 ')
Flow('mute-agc','mute-a','agc rect1=25 rect2=25')
Flow('mute-fact','mute-agc mute-a','divn den=${SOURCES[1]} rect1=%i rect2=%i'%(drect1,drect3))
#Flow('mute-fact','mute-a mute-agc ','scale dscale=1000 | divn den=${SOURCES[1]} rect1=%i rect2=%i'%(drect1,drect3))
Flow('mute','mute-fact mute-pre','spray axis=2 n=%i o=%g d=%g |add mode=p ${SOURCES[1]} '%(nv,v0+dv,dv))



Result('mute','window n3=1 min3=1200 |grey color=j min1=%g max1=%g allpos=y'%(min1,max1))

# optimal rect1=50 rect2=20
# rect1=30 rect2=10
picksemb = '''
        pick rect1=50 rect2=20 vel0=1480
        '''

Flow('semb-mute','semb',mutesemb)

Flow('vpick-semb','mute',picksemb)

Result('vpick-semb',plotvel('OVC Picked Velocity'))

# do some probabilistic stuff
Flow('mute-stk','mute','stack axis=2')


Flow('v-exp','mute mute-stk',
   '''
   math output="input*x2" | 
   stack axis=2 | 
   divnp den=${SOURCES[1]} rect1=%i rect2=%i | 
   clip2 lower=%g 
   '''%(vrect1,vrect2,vmin))
Result('nankai-v-exp','v-exp',
   '''
   grey title="Expectation Velocity" 
   color=j allpos=y scalebar=y bias=%g
   min1=%g max1=%g barunit="m/s" barlabel="RMS Velocity"
   '''%(vmin,min1,max1))
# convert to dix

Flow('v-dix','v-exp','dix rect1=%g rect2=%g | clip2 lower=%g'%(dixrect1,dixrect2,vmin))
Result('nankai-v-dix','v-dix',
   '''
   grey title="Dix Velocity"
   min1=%g max1=%g
   scalebar=y allpos=y bias=%g color=j barunit="m/s" barlabel="Velocity"
   '''%(min1,max1,vmin))
Result('nankai-v-dix-nbar','v-dix',
   '''
   grey title= 
   min1=%g max1=%g
   scalebar=n allpos=y bias=%g color=j barunit="m/s" barlabel="Velocity"
   '''%(min1,max1,vmin))
# determine the variance
Flow('v-exp-var','v-exp mute mute-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 |
   clip2 lower=%g
   '''%(nv,dv,v0+dv,vrect1,vrect2,dv*dv))

Result('nankai-v-exp-var','v-exp-var',
   '''
   grey title="Expectation Velocity Variance" 
   color=j allpos=y scalebar=y barunit="m\^2\_/s\^2" barlabel="Variance"
   min1=%g max1=%g 
   '''%(min1,max1))

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" |
   divnp 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','mute v-exp','slice pick=${SOURCES[1]} ')

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

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

Flow('semb-eps','semb',
    '''math output="input+`<${SOURCE} $RSFROOT/bin/sfattr want=std | awk '{print $4}'`" ''')
Flow('dsemb','semb-dz semb-dx semb-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')

# combined weights
Flow('wts','mute d-vel dsemb d-vel',
   '''
   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))

# weighted gathers
Flow('wtd-gath','wts tvxp',
   '''
   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))

# constant gather
Flow('const-gath','tvxp','stack axis=2 | transp plane=23 memsize=10000')

Flow('wtd-img','wts tvx',' add mode=p ${SOURCES[1]}')
#Flow('wtd-img','wts tvx','math output="input^(1/3)"| add mode=p ${SOURCES[1]}')
Flow('prob-dimage','wtd-img','stack axis=2')

pdimgclip = 99.95
Result('nankai-prob-dimage','prob-dimage',
   '''
   grey title="Probabilistic Diffraction Image" 
   min1=%g max1=%g 
   '''%(min1,max1))
Flow('prob-dimage-squared','prob-dimage','add mode=p ${SOURCE}')

# 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=.05
   '''%(nv,dv,v0+dv,drect1,drect3))

Result('nankai-prob-dimage-var','prob-dimage-var',
   '''
   grey color=j allpos=y scalebar=y
   title="Image Variance" 
   min1=%g max1=%g pclip=99.9 bias=.05
   '''%(min1,max1))

# create expectation velocity slices
Flow('vexp-slice','ovc v-exp','transp plane=23 memsize=1000 | slice pick=${SOURCES[1]}')

# create slices for deterministic images
Flow('vel-spray','v-exp','spray axis=3 n=%i d=%g o=%g'%(np,dp,p0))
Flow('slice-tpx','tvxp vel-spray','slice pick=${SOURCES[1]} | transp plane=23 memsize=10000')
# full image
Flow('tvxp-full','txvp-full','transp plane=23 memsize=10000')
Flow('full-slice-tpx','tvxp-full vel-spray','slice pick=${SOURCES[1]} | transp plane=23 memsize=10000')
# stack to make image
Flow('full-img','full-slice-tpx','stack axis=2')
# compute dip from image

# warp image to squared time
Flow('full-img2','full-img','t2warp pad=%i'%pad)
dt2 =0.0583169 # from warp to squared time
dt2dx = dt2/dx
# compute dip, warp back
Flow('full-img-slope','full-img2',
   '''
   dip rect1=%i rect2=%i|
   t2warp inv=y |
   scale dscale=%g
   '''%(20,10,dt2dx))
#   math output="sqrt(%g*%g)*sqrt(abs(input))*input/(abs(input)+0.0000000000000001)"
wtlst = ['tvx','wts','wtd-img','mute','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','full-slice-tpx']
gathtitles = ['Probabilistic Weight Gather','Deterministic Gather','Equal Weight Gather','Complete Image 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]
clipss = [99,99,99,99,100,99]

gmax1 = 6
gmin1 = 7.6
gmax2 = nv*dv+v0
gmin2 = v0+dv

nz = 1601
dz = 4
zo = vmin*(min1)
zo = 0
dx = 16.667
ox = 900

point1 = .8
point2 = .7
point1d = 0.8
point2d = 0.5
minslope = 0.01
# plotcol 4 is fuschia
for n in range(6):
    #x = (3,5,7)[n]
    x = (1200, 2700, 3000, 4100, 5000, 6500)[n]
    plst  = []
    plst3 = []
    gathlst = []
#math.sqrt(dt)/dx
    Plot('img-slope-%i'%n,'full-img-slope',
       '''
       window min2=%g n2=1 |
       graph min1=%g max1=%g min2=%g max2=%g 
       label1= label2= unit1= unit2=  title= 
       transp=y n1tic= 0 plotfat=8 plotcol=4 dash=1
       '''%(x,gmin1,gmax1,-1*minslope,minslope))
    for i in range(len(gathers)):
        gath = gathers[i]
        gtitle = gathtitles[i]
        Result('nankai-'+gath+'-%i'%n,gath,
             '''
             window n3=1 min3=%g min2=%g max2=%g|
             grey min1=%g max1=%g title="%s"
             label1=Time unit1=s label2=Slope unit2="s\^2\_/m" 
             pclip=99.8
             '''%(x,-1*minslope, minslope, gmax1,gmin1,gtitle))
        Plot('nankai-'+gath+'-%i'%n,gath,
             '''
             window n3=1 min3=%g min2=%g max2=%g|
             grey min1=%g max1=%g title="%s"
             label1=Time unit1=s label2=Slope unit2="s\^2\_/m" 
             pclip=99.8
             '''%(x,-1*minslope, minslope, gmax1,gmin1,gtitle))
        Result('nankai-'+gath+'-slope-%i'%n,['nankai-'+gath+'-%i'%n,'img-slope-%i'%n],'Overlay')

    # get the velocity, graph
    Flow('vpick-semb-%i'%n,'vpick-semb','window n2=1 min2=%g '%x)
    Plot('vpick-semb-%i'%n,
       '''
       graph min1=%g max1=%g min2=%g max2=%g 
       label1= label2= unit1= unit2=  title= 
       transp=y n1tic=0 dash=2 plotfat=20 
       '''%(gmin1,gmax1,gmin2,gmax2))
    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=7
       '''%(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
           min1=%g max1=%g unit2="m/s"
           color=%s allpos=%s  pclip=%g
           '''%(titles[k],x,gmax1,gmin1,colorlst[k],allpos[k],clipss[k]))
        Plot(item+'3-%i'%n,item,
           '''
           window min1=%g max1=%g |
           byte gainpanel=2 pclip=%g allpos=%s |
           grey3 unit2="m/s" label2=Velocity 
           color=%s title="%s" flat=n
           frame1=%i frame2=%i frame3=%i 
           point1=%g point2=%g 
           color=%s allpos=%s
           screenht=15 screenratio=1.2
           titlesz=16 labelsz=8 
           '''%(gmax1,gmin1,clipss[k],allpos[k],allpos[k],
                titles[k],140,nv/3-1,(x-ox)/dx+1,point1d,point2d,
                colorlst[k],allpos[k]))
        Result(item+'3-%i'%n,item,
           '''
           window min1=%g max1=%g |
           byte gainpanel=2 pclip=%g allpos=%s |
           grey3 unit2="m/s" label2=Velocity 
           color=%s title="%s" flat=n
           frame1=%i frame2=%i frame3=%i 
           point1=%g point2=%g 
           color=%s allpos=%s
           screenht=28 screenratio=2
           larnersz=85 titlesz=16
           '''%(gmax1,gmin1,clipss[k],allpos[k],allpos[k],
                titles[k],140,nv/3-1,(x-ox)/dx,point1,point2,
                colorlst[k],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,'vpick-semb-%i'%n],'Overlay')
    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('nankai-weights-%i-a'%n,[plst[0],plst[1],plst[2]],'SideBySideAniso')
    Result('nankai-weights-%i-b'%n,[plst[3],plst[4]+'-o',plst[5]],'SideBySideAniso')
    Result('nankai-weights3-%i-a'%n,[plst3[0],plst3[1],plst3[2]],'SideBySideIso')
    Result('nankai-weights3-%i-b'%n,[plst3[3],plst3[4],plst3[5]],'SideBySideIso')
    Result('nankai-weights3a-%i'%n,plst3,'TwoRows')
# 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)


picking = []

pickingm = []

picking1 = []

picking1m = []

for loc in (1200, 2700, 3000, 5000, 6500):

        Flow('semb-%g'%loc,'semb','window n3=1 min3=%g min1=%g max1=%g' % (loc,min1,max1))

        Flow('semb-m-%g'%loc,'semb-mute','window n3=1 min3=%g min1=%g max1=%g' % (loc,min1,max1))

        Plot('semb-%g'%loc,
                '''
                grey color=j allpos=y title="CMP %g m" 
                label2=Velocity unit2=m/s labelsz=12. titlesz=12.
                d1num=250 n1tic=5 o1num=1450 d2num=0.5 n2tic=6 o2num=5.5
                '''%loc)

        Plot('semb-m-%g'%loc,
                '''
                grey color=j allpos=y title="CMP %g m" 
                label2=Velocity unit2=m/s labelsz=12. titlesz=12. d2num=0.25 n2tic=4 o2num=5.5
                '''%loc)

        Flow('vpick-semb-%g'%loc,'vpick-semb','window n2=1 min2=%g min1=%g max1=%g' % (loc,min1,max1))

        Plot('vpick-semb-%g'%loc,plotpick())

        Plot('vpick-semb-plot-%g'%loc,'semb-%g vpick-semb-%g'%(loc,loc),'Overlay')

        picking.append('vpick-semb-plot-%g'%loc)

        Plot('vpick-semb-plot-m-%g'%loc,'semb-m-%g vpick-semb-%g'%(loc,loc),'Overlay')

        pickingm.append('vpick-semb-plot-m-%g'%loc)

        # picking on a single location without smoothing

        Flow('vpick-semb-1-%g'%loc,'semb-mute','window n3=1 min3=%g  | '%loc + picksemb + '| window min1=%g max1=%g'% (min1,max1))

        Plot('vpick-semb-1-%g'%loc,plotpick())

        Plot('vpick-semb-plot-1-%g'%loc,'semb-%g vpick-semb-1-%g'%(loc,loc),'Overlay')

        picking1.append('vpick-semb-plot-1-%g'%loc)

        Plot('vpick-semb-plot-m-1-%g'%loc,'semb-m-%g vpick-semb-1-%g'%(loc,loc),'Overlay')

        picking1m.append('vpick-semb-plot-m-1-%g'%loc)

#Result('picking',picking,'SideBySideAniso')
#Result('picking-m',pickingm,'SideBySideAniso')

Result('g-picking','vpick-semb-plot-1-2700 vpick-semb-plot-1-5000 vpick-semb-plot-1-6500','SideBySideAniso')

# suspicious behaviour for the first cmp being considered
#Result('picking-1',picking1,'SideBySideAniso')
#Result('picking-1-m',picking1m,'SideBySideAniso')      

# Semblance velocity
#<ovc.rsf sfwindow min2=4100 max2=4400 min1=6.0 max1=7.0 | sfgrey screenratio=6.0 labelsz=2 xll=5.0 | sfpen

Flow('vc-slice-semb','ovc vpick-semb','transp plane=23 memsize=1000 | slice pick=${SOURCES[1]}')

Result('vc-slice-semb',section('Diffraction Image'))

# path summation image
Flow('path-sum','ovc','stack axis=3')
Result('nankai-path-sum','path-sum',
   '''
   grey title="Equal Weight Diffraction Image"
   min1=%g max1=%g
   '''%(min1,max1))


Flow('full-image','full-slice-tpx','stack axis=2 | bandpass flo=5 fhi=45')

Result('nankai-deterministic','vexp-slice',
   '''
   grey title="Deterministic Diffraction Image"
   min1=%g max1=%g
   '''%(min1,max1))
### Perform velocity continuation on dmo stack
Flow('vc-fw','slice',
        '''
        cosft sign2=1 |
        vczo nv=1 dv=%g v0=0.0 | window |
        vczo nv=%d dv=%g v0=%g | 
        transp plane=23 memsize=1000 |
        cosft sign2=-1
        '''%(v0,nv,dv,v0))

Flow('vc-fw-slice-semb','vc-fw vpick-semb','transp plane=23 memsize=1000 | slice pick=${SOURCES[1]}')

Result('vc-fw-slice-semb',section('Conventional Image'))

#Flow('vc-full-slice','vc-fw v-exp','transp plane=23 memsize=1000 | slice pick=${SOURCES[1]}')
Flow('vc-full-slice','full-image','cp')
Result('nankai-dif','dif',
   '''
   grey title="Diffraction Data"
   min1=%g max1=%g
   '''%(min1,datamax1))
Result('nankai-data','slice',
   '''
   grey title="Complete Data"
   min1=%g max1=%g
   '''%(min1,datamax1))
Result('nankai-full','vc-full-slice',
   '''
   grey title="Complete Image"
   min1=%g max1=%g
   '''%(min1,max1))
Result('nankai-full-2','full-img',
   '''
   grey title="Complete Image"
   min1=%g max1=%g
   '''%(min1,max1))
### Plotting gathers as in BEI

#plot3 = ' put d1=0.004 d2=0.0335 o2=0 label2=Midpoint unit2=km unit1=s| window  min2=.500 max2=7.750 max1=1.7 min1=1 max3=2 min3=-2 |transp plane=23 memsize=5000| byte|grey3 frame1=60 frame3=95 frame2=69 label1=Time flat=n unit2="s\^2\_/km"'

# frames for interesting location
# cmp 2700 [m] - 108th sample
# TWTT 6.5 [s] - 621th sample  

ovcmin1=6.0
ovcmax1=7.0

ovcmin2=3000#2000
ovcmax2=5000#4000

pmin=-0.004#-0.005
pmax=0.004# 0.005

#dp = 0.000186667

frameovc1 = (int)((6.5 - ovcmin1)/0.004)
frameovc2 = (int)((pmax)/(dp))#(int)(np/2)
frameovc3 = (int)((4100.06 - ovcmin2)/16.667) #2700


# List of interesting places: diffractions
# fault zone left part
# v = 1500 m/s, t = 6.5 s
x1 = 2666.70 #m
# diffraction central part
# v = 1500 m/s, t = 6.4 s
x2 = 4100.06 #m
# weaker diffraction central part
# v = 1480 m/s, t = 6.2 s
x3 = 4216.73 #m 




# Just slope decomposed volume

Result('tpx','fpx',
        '''
        fft1 inv=y | t2warp inv=y | 
        window min1=%g max1=%g min2=%g max2=%g min3=%g max3=%g |
        byte |
        grey3 frame1=%d frame2=%d frame3=%d label1=Time flat=n unit2="s\^2\_/m"
        title="Slope Decomposed Diffraction Data" d1num=0.004 n1tic=3 o1num=-0.004 o3num=%g n3tic=4 d3num=500
        '''%(ovcmin1,ovcmax1,pmin,pmax,ovcmin2,ovcmax2,frameovc1,frameovc2,frameovc3,ovcmin2+500))

Flow('vpick-semb-slice-mig','vpick-semb','window min2=%g max2=%g | spray axis=3 n=%d d=0.000186667 o=%g'%(ovcmin2,ovcmax2,np,p0))

Flow('tvxp','txvp','transp plane=23 memsize=5000')

Flow('txp-migrated','tvxp vpick-semb-slice-mig','window min3=%g max3=%g | slice pick=${SOURCES[1]}'%(ovcmin2,ovcmax2))

Result('txp-migrated',
        '''
        transp plane=23 |
        window min1=%g max1=%g min2=%g max2=%g min3=%g max3=%g |
        byte |
        grey3 frame1=%d frame2=%d frame3=%d label1=Time flat=n unit2="s\^2\_/m"
        title="Slope Decomposed Diffraction Image" d1num=0.004 n1tic=3 o1num=-0.004 o3num=%g n3tic=4 d3num=500
        '''%(ovcmin1,ovcmax1,pmin,pmax,ovcmin2,ovcmax2,frameovc1,frameovc2,frameovc3,ovcmin2+500))

# time-to-depth conversion (Sripanich and Fomel, 2017)
# from migration velocity to Dix velocity

padbeg1 = ot1/dt
Flow('v-exp-mir','v-exp','reverse which=2 opt=i')
Flow('v-exp-ext','v-exp-mir v-exp','cat axis=2 ${SOURCES[1]} ${SOURCES[0]}')
vdixmin = 1485
Flow('v-exp-ext-pd','v-exp-ext','window n1=1 | spray axis=1 n=%i d=%g o=%g | cat axis=1 ${SOURCES[0]}'%(padbeg1,dt,0))
Flow('vdix-pre','v-exp-ext-pd',
     '''
     dix rect1=%i rect2=%i  niter=80 |
     window n2=%i f2=%i
     '''%(dixrect1,dixrect2,nx,nx))
# make a superior vdix figure, the limits of diffractions in the deep lead to funny looking and distracting velocities
# therefore, just clip the lower portion to make sure things dont get too weird.
dixmin1 = 7.
dixlowminv = 2380
Flow('vdix-bot','vdix-pre','window min1=%g | clip2 lower=%g'%(dixmin1,dixlowminv))
Flow('vdix','vdix-pre vdix-bot','window max1=%g | cat axis=1 ${SOURCES[1]} | put o2=%g'%(dixmin1-dt,ox))
#Flow('v-dix-bot','v-dix','window min1=%g | clip2 lower=%g'%(dixmin1,dixlowminv))
#Flow('v-dix2','v-dix v-dix-bot','window max1=%g | cat axis=1 ${SOURCES[1]}'%(dixmin1-dt))
Flow('v-dix2','vdix','window min1=%g'%ot1)
timedepth = '''
   time2depth velocity=${SOURCES[1]} nz=%i dz=%g z0=%g intime=y twoway=y |
   put label1=Depth unit1=m
   '''%(nz,dz,zo)
# image to depth
timeimg = 'full-img'#'vc-full-slice'
Flow(timeimg+'-d',[timeimg,'vdix'],timedepth)
# convert vdix to depth
Flow('dixdepth','vdix vdix',timedepth)
# derivatives of dix velocity squared     
Flow('dv2dt0','vdix','math output="input^2" | smoothder')
Flow('dv2dx0','vdix','math output="input^2" | transp | smoothder | transp')
Flow('beta','dv2dt0 vdix',timedepth)
Flow('alpha','dv2dx0 vdix',timedepth)
# Reference velocity squared
Flow('refdix','dixdepth','math output="input^2" ')
Flow('refvz','dixdepth',
   '''
   window n2=1 f2=%i| 
   math output="input^2" | 
   spray axis=2 n=%i o=%g d=%g
   '''%((nx-1)/2,nx,ox,dx))
Flow('depth dx0 dt0 dv','refdix refvz alpha beta',
        '''
        time2depthweak zsubsample=50 nsmooth=16
        velocity=$SOURCE refvelocity=${SOURCES[1]} dvdx0=${SOURCES[2]} dvdt0=${SOURCES[3]}
        outdx0=${TARGETS[1]} outdt0=${TARGETS[2]} outdv=${TARGETS[3]}
        ''')
Flow('finalv','refvz dv','math est=${SOURCES[1]} output="sqrt(input+est)" | put d3=1 o3=0 ')
Flow('diff1','finalv dixdepth','math est=${SOURCES[1]} output="input-est" | put d3=1 o3=0 ')
Flow('diff','finalv dixdepth','math est=${SOURCES[1]} output="input^2-est^2" | put d3=1 o3=0 ')
Flow('reft0','refvz','math output="2*1/sqrt(input)*%g" | causint '%(dz)) # two-way
Flow('finalt0','reft0 dt0','math dt=${SOURCES[1]} output="input+dt"')

Flow('refx0','refvz','math output="x2"')
Flow('finalx0','refx0 dx0','math dx=${SOURCES[1]} output="input+dx"')

Flow('dixcoord','reft0 refx0',
     '''
     cat axis=3 ${SOURCES[1]} |
     transp plane=23 | transp plane=12
     ''')

Flow('finalcoord','finalt0 finalx0',
     '''
     cat axis=3 ${SOURCES[1]} |
     transp plane=23 | transp plane=12
     ''')

timeimg = 'vc-full-slice'
difimg = 'prob-dimage'
name = 'nankai'
dmin1 = 4.2
dmax1 = 6.2
dmin2 = ox+5*dx
dmax2 = (nx-5)*dx+ox
screenratio = (dmax1-dmin1)*1000/(dmax2-dmin2)
screenht = 4.4
screenht2 = 4
imglst = [timeimg,difimg,'v-dix2','vdix']
titlst = ['Complete Depth Image','Diffraction Depth Image','Dix Velocity in Depth','Dix Velocity in Depth']
colorlst = ['i','i','j','j']
allposlst = ['n','n','y','y']
biaslst = [0,0,vmin,vmin]
extra = ['','',' pclip=96 ',' clip=96']
pcliplst = [99,99,99,99]
# contour plot
Plot('finalcoord-t','finalcoord',
   'put d2=.004 |window n1=1 | contour title= min1=%g max1=%g min2=%g max2=%g unit1= label1= dc=.250 plotcol=1'%(dmin1,dmax1,dmin2,dmax2))
Plot('finalcoord-x',
   'finalcoord','put d2=.004 |window n1=1 f1=1 | contour title= min1=%g max1=%g min2=%g max2=%g dc=500 unit1= label1= plotcol=1'%(dmin1,dmax1,dmin2,dmax2))
#Flow('dixmapd',[timeimg,'dixcoord'],
#   'window min1=%g | put o1=0 | interp=spline nw=8 coord=${SOURCES[1]} | put label1=Depth unit1=m '%min1)
to = 4
larner = 35
titlsz = 6
for i1 in range(len(imglst)):
    img = imglst[i1]
    title = titlst[i1]
    Flow(img+'-depth',[img,'finalcoord'],
       'window min1=%g |inttest2 interp=spline nw=8 coord=${SOURCES[1]}| put d1=%g label1=Depth unit1=km'%(to,dz/1000))
    # depth window
    Result(name+'-'+img+'-depth',img+'-depth',
       '''
       window min1=%g max1=%g|
       grey title="%s" 
       min2=%g max2=%g 
       screenratio=%g screenht=%g  
       color=%s allpos=%s bias=%g %s
       larnersz=%g titlesz=%g
       '''%(dmin1,dmax1,title,dmin2,dmax2,
       screenratio,screenht,colorlst[i1],
       allposlst[i1],biaslst[i1],extra[i1],larner,titlsz))

    Plot(name+'-'+img+'-depth',img+'-depth',
       '''
       window min1=%g max1=%g|
       grey title= "%s" 
       min2=%g max2=%g 
       screenratio=%g screenht=%g  
       color=%s allpos=%s bias=%g %s
       larnersz=%g titlesz=%g
       '''%(dmin1,dmax1,title,dmin2,dmax2,
       screenratio,screenht,colorlst[i1],
       allposlst[i1],biaslst[i1],extra[i1],larner,titlsz))
    Plot(name+'-'+img+'-depth-w',img+'-depth',
       '''
       window min1=%g max1=%g j2=3 |
       wiggle title= "%s" transp=y
       min2=%g max2=%g min1=%g max1=%g
       screenratio=%g screenht=%g  plotcol=0
       color=%s allpos=%s bias=%g pclip=%g n2tic=0 grid1=n grid2=n n1tic=0
       larnersz=%g titlesz=%g label2= unit2=
       '''%(dmin1,dmax1,title,dmin2,dmax2,dmax1,dmin1,
       screenratio,screenht,colorlst[i1],
       allposlst[i1],biaslst[i1],pcliplst[i1],larner,titlsz))
    Plot(name+'-'+img+'-depth-w2',img+'-depth',
       '''
       window min1=%g max1=%g j2=3 |
       wiggle title="%s" transp=y
       min2=%g max2=%g min1=%g max1=%g
       screenratio=%g screenht=%g  plotcol=0
       color=%s allpos=%s bias=%g pclip=%g 
       crowd1=.742 grid1=n 
       larnersz=%g titlesz=%g grid2=n
       '''%(dmin1,dmax1,title,dmin2,dmax2,dmax1,dmin1,
       screenratio,screenht2,colorlst[i1],
       allposlst[i1],biaslst[i1],pcliplst[i1],larner,titlsz))
    Plot(name+'-'+img+'-depth-bar',img+'-depth',
             '''
             clip2 lower=%g |
             scale dscale=.001 |
             window min1=%g max1=%g|
             grey title= "%s" 
             min2=%g max2=%g 
             screenratio=%g screenht=%g  
             color=%s allpos=%s bias=%g %s crowd1=.75
             scalebar=y barwidth=0.335 barlabel="Dix Velocity" barunit="km/s" 
             larnersz=%g titlesz=%g n2tic=0  n1tic=0 label2= unit2= 
             '''%(biaslst[i1],dmin1,dmax1,title,dmin2,dmax2,
             screenratio/1.11,screenht2,colorlst[i1],
             allposlst[i1],biaslst[i1]*.001,extra[i1],larner,titlsz))
    Result(name+'-'+img+'0'+'-depth',img+'-depth',
       '''
       window min1=%g max1=%g |
       grey title= "%s" 
       min2=%g max2=%g 
       color=%s allpos=%s bias=%g %s
       larnersz=%g titlesz=%g
       '''%(dmin1,dmax1,title,dmin2,dmax2,
       colorlst[i1],allposlst[i1],biaslst[i1],
       extra[i1],larner,titlsz))
    Plot(name+'-'+img+'0'+'-depth',img+'-depth',
          '''
          window min1=%g max1=%g |
          grey title= "%s" 
          min2=%g max2=%g 
          color=%s allpos=%s bias=%g %s
          larnersz=%g titlesz=%g
          '''%(dmin1,dmax1,title,dmin2,dmax2,
          colorlst[i1],allposlst[i1],biaslst[i1],
          extra[i1],larner,titlsz))
    Result(name+'-'+img+'0'+'-depth-cont',[name+'-'+img+'0'+'-depth','finalcoord-t','finalcoord-x'],'Overlay')
#    Result(name+'-'+img+'0'+'-depth-cont0',[name+'-'+img+'0'+'-depth','finalcoord-t','finalcoord-x'],'Overlay')


Result(name+'-'+timeimg+'depthwiggles',[name+'-'+'v-dix2'+'-depth',name+'-'+timeimg+'-depth-w'],'Overlay')
Result(name+'-'+difimg+'depthwiggles',[name+'-'+'v-dix2'+'-depth',name+'-'+difimg+'-depth-w'],'Overlay')
Result(name+'-'+timeimg+'depthwiggles-bar',[name+'-'+'v-dix2'+'-depth-bar',name+'-'+timeimg+'-depth-w2'],'Overlay')
Result(name+'-'+difimg+'depthwiggles-bar',[name+'-'+'v-dix2'+'-depth-bar',name+'-'+difimg+'-depth-w2'],'Overlay')



End()