from rsf.proj import *
import math
from rsf.recipes.tpx import FPX

### 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')

####### 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')

# find a term
Flow('recvsc',['arms','extindrecv',sc,'recvmodel'],
     '''
     conjgrad ./${SOURCES[2]} index=${SOURCES[1]} 
     mod=${SOURCES[3]} niter=150
     ''')

# project to a data space
Flow('recvscarms',['recvsc','extindrecv',sc],
     './${SOURCES[2]} 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.5
max1=8.0

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')

rect1=30
rect2=20
rect3=30

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

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

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

Flow('pwd-vc','vc15 dip-vc','pwd dip=${SOURCES[1]}')

#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',
        '''
        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))

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

#?# -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)

Flow('txp','fpx','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('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 = 60
dv = 20
v0 = 1400

Flow('ovc-1st-step','fkp','ovczo 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','ovczo 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=100 inv=y |
        t2warp inv=y
        ''',split=[4,np])

# 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('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

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 
        '''
# 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','semb-mute',picksemb)

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

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'))

### 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'))

### 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

#Result('tpx15','txvp',
#       '''
#       window min1=%g max1=%g n3=1 min3=1500 min2=%g max2=%g |
#       transp plane=23 memsize=1000 |
#       byte |
#       grey3 frame1=%d frame2=%d frame3=%d label1=Time flat=n unit2="s\^2\_/km"
#       '''%(ovcmin1,ovcmax1,ovcmin2,ovcmax2,frameovc1,frameovc2,frameovc3))

#Result('tpx25','txvp',
#       '''
#       window min1=%g max1=%g n3=1 min3=2500 min2=%g max2=%g |
#       transp plane=23 memsize=1000 |
#       byte |
#       grey3 frame1=%d frame2=%d frame3=%d label1=Time flat=n unit2="s\^2\_/km"
#       '''%(ovcmin1,ovcmax1,ovcmin2,ovcmax2,frameovc1,frameovc2,frameovc3))

# 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 

dx=200 #m

for case in (1,2,3):

        if case == 1:
                name = 'left'
                x = x1
        if case == 2:
                name = 'center-s'
                x = x2
        if case == 3:
                name = 'center-w'
                x = x3

        locxmin = x-dx
        locxmax = x+dx

        ### Gather with picked velocity

        Flow('vpick-semb-slice','vpick-semb','window n2=1 min2=%g | spray axis=2 n=%d d=0.000186667 o=%g'%(x,np,p0))

        Flow('txp-slice-'+name,'txvp vpick-semb-slice',
                '''
                window n2=1 min2=%g |
                slice pick=${SOURCES[1]}
                '''%x)

        Result('txp14-slice-'+name,'txvp',
                '''
                window n2=1 min2=%g n3=1 min3=1400 |
                grey min1=%g max1=%g min2=%g max2=%g unit2="s\^2\_/m" title="Slope Gather for v = 1.4 km/s"
                '''%(x,ovcmin1,ovcmax1,pmin,pmax))

        Result('txp25-slice-'+name,'txvp',
                '''
                window n2=1 min2=%g n3=1 min3=2500 |
                grey min1=%g max1=%g min2=%g max2=%g unit2="s\^2\_/m" title="Slope Gather for v = 2.5 km/s"
                '''%(x,ovcmin1,ovcmax1,pmin,pmax))

        Result('txp-slice-'+name,
                '''
                grey min1=%g max1=%g min2=%g max2=%g unit2="s\^2\_/m" title="Slope Gather for Migration Velocity"
                '''%(ovcmin1,ovcmax1,pmin,pmax))

        Result('image-'+name,'vc-slice-semb',section('Diffraction Image','Time','s',6.0,7.0,'min2=%g max2=%g'%(locxmin,locxmax)))

Plot('txp14-slice','txvp','window n2=1 min2=2700 n3=1 min3=1400 | grey min1=%g max1=%g min2=%g max2=%g wanttitle=n'%(ovcmin1,ovcmax1,pmin,pmax))

Plot('txp25-slice','txvp','window n2=1 min2=2700 n3=1 min3=2500 | grey min1=%g max1=%g min2=%g max2=%g wanttitle=n'%(ovcmin1,ovcmax1,pmin,pmax))

#Plot('txp-slice','grey min1=%g max1=%g min2=%g max2=%g wanttitle=n'%(ovcmin1,ovcmax1,pmin,pmax))

#Result('cigs','txp14-slice txp25-slice txp-slice','SideBySideAniso')

###

# check illumination

for cignum in (2700, 4500, 6500):

        plotillum14 = 'window n2=1 min2=%g n3=1 min3=%g | grey min1=%g max1=%g wanttitle=n'%(cignum,1400,ovcmin1,ovcmax1)

        plotillum25 = 'window n2=1 min2=%g n3=1 min3=%g | grey min1=%g max1=%g wanttitle=n'%(cignum,2500,ovcmin1,ovcmax1)

        plotillum2 = 'grey min1=%g max1=%g wanttitle=n'%(ovcmin1,ovcmax1)

        Plot('txp14-slice-%g'%cignum,'txvp',plotillum14)

        Plot('txp25-slice-%g'%cignum,'txvp',plotillum25)

        Flow('txp-slice-%g'%cignum,'txvp vpick-semb-slice',
                '''
                window n2=1 min2=%g |
                slice pick=${SOURCES[1]}
                '''%cignum)

        Plot('txp-slice-%g'%cignum,plotillum2)

#Result('cigs-2700','txp14-slice-2700 txp25-slice-2700 txp-slice-2700','SideBySideAniso')

#Result('cigs-4500','txp14-slice-4500 txp25-slice-4500 txp-slice-4500','SideBySideAniso')

#Result('cigs-6500','txp14-slice-6500 txp25-slice-6500 txp-slice-6500','SideBySideAniso')

# 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))

#Result('txp-slice',
#       '''
#       window min1=%g max1=%g min2=%g max2=%g |
#       transp plane=23 memsize=1000 |
#       byte |
#       grey3 frame1=%d frame2=%d frame3=%d label1=Time flat=n unit2="s\^2\_/km"
#       '''%(ovcmin1,ovcmax1,ovcmin2,ovcmax2,frameovc1,frameovc2,frameovc3))

End()