from rsf.proj import *
from rsf.prog import RSFROOT

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

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

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

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

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

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

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


def Grey1(data,data0,other):
    Result(data,data0,'''
    put d2=1 o2=1|grey 
    label1=Time unit1=s color=g label2=Distance unit2=m clip=0.1 wanttitle=n title= screenratio=1.2 min1=5.5 max1=8.5 label2=Trace unit2= max2=300
    %s ''' % (other))

def Grey(data,other):
    Result(data,'''
    put d2=1 o2=1|grey
    label1=Time unit1=s color=g label2=Distance unit2=m clip=0.1 wanttitle=n title= screenratio=1.2 min1=5.5 max1=8.5 label2=Trace unit2= max2=300
    %s ''' % (other))

def Greyz(data,other):
    Result(data,'''
    grey
    label1=Time unit1=s color=g label2=Distance unit2=m clip=0.1 wanttitle=n title= screenratio=1.2 label2=Trace unit2= screenratio=1.0
    %s ''' % (other))

def Greyv(data,data0,other):
    Result(data,data0,
        '''
        grey 
        color=j scalebar=y bias=%g barlabel=Velocity 
        barreverse=y %s
        ''' % (v0+0.5*nv*dv,other))

#Flow('k-post','vc15','cp|scale axis=2')
Flow('k-post','vc15','cp|scale axis=2')

Grey('k-post','')

Flow('k-dip','k-post','dip rect1=10 rect2=10')
Flow('k-pwd-n','k-post k-dip','pwd dip=${SOURCES[1]}')
Flow('k-pwd-s','k-post k-pwd-n','add scale=1,1 ${SOURCES[1]}')

Grey('k-pwd-s','')
Grey('k-pwd-n','clip=0.05')

########################################################################
# INITIALIZATION
########################################################################
matlab         = WhereIs('matlab')
matROOT = '../matfun/'
matfun = 'LDRR'
matlabpath = os.environ.get('MATLABPATH',os.path.join(RSFROOT,'lib'))

if not matlab:
    sys.stderr.write('\nCannot find Matlab.\n')
    sys.exit(1)

n1=1750
n2=401
n3=1
d1=0.004
d2=16.667
d3=1
o1=4
o2=900
o3=0
lf=0
hf=120
N=5
verb=0

n1win=100
n2win=20
n3win=1
r1=0.5
r2=0.5
r3=0.5

put='n1=%d n2=%d n3=%d d1=%g d2=%g d3=%g o1=%g o2=%g o3=%g'%(n1,n2,n3,d1,d2,d3,o1,o2,o3)

N=5
NN=3
ifdamp=0
############################################################
## with parameter
############################################################
Flow('k-lrr0 k-lrra0',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}','${TARGETS[1]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(NN)d,%(n1win)d,%(n2win)d,%(n3win)d,%(r1)g,%(r2)g,%(r3)g,%(ifdamp)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-lrr','k-lrr0','put %s'%put)
Flow('k-lrra','k-lrra0','put %s'%put)
Flow('k-lrr-n','k-post k-lrr','add scale=1,-1 ${SOURCES[1]}')
Flow('k-lrra-n','k-post k-lrra','add scale=1,-1 ${SOURCES[1]}')


NN=3
ifdamp=1
Flow('k-ldrr0 k-ldrra0',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}','${TARGETS[1]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(NN)d,%(n1win)d,%(n2win)d,%(n3win)d,%(r1)g,%(r2)g,%(r3)g,%(ifdamp)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-ldrr','k-ldrr0','put %s'%put)
Flow('k-ldrra','k-ldrra0','put %s'%put)
Flow('k-ldrr-n','k-post k-ldrr','add scale=1,-1 ${SOURCES[1]}')
Flow('k-ldrra-n','k-post k-ldrra','add scale=1,-1 ${SOURCES[1]}')

Grey('k-lrr','title="LRR"')
Grey('k-lrra','title="LRRA"')
Grey('k-ldrr','title="LDRR"')
Grey('k-ldrra','title="LDRRA"')

Grey('k-lrr-n','title="LRR" clip=0.05')
Grey('k-lrra-n','title="LRRA" clip=0.05')
Grey('k-ldrr-n','title="LDRR" clip=0.05')
Grey('k-ldrra-n','title="LDRRA" clip=0.05')

# LRR with fixed N
N=1
ifdamp=0
Flow('k-lrr-N1-0 k-lrra0-tmp1',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}','${TARGETS[1]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(NN)d,%(n1win)d,%(n2win)d,%(n3win)d,%(r1)g,%(r2)g,%(r3)g,%(ifdamp)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-lrr-N1','k-lrr-N1-0','put %s'%put)
Flow('k-lrr-N1-n','k-post k-lrr-N1','add scale=1,-1 ${SOURCES[1]}')

N=2
ifdamp=0
Flow('k-lrr-N2-0 k-lrra0-tmp2',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}','${TARGETS[1]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(NN)d,%(n1win)d,%(n2win)d,%(n3win)d,%(r1)g,%(r2)g,%(r3)g,%(ifdamp)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-lrr-N2','k-lrr-N2-0','put %s'%put)
Flow('k-lrr-N2-n','k-post k-lrr-N2','add scale=1,-1 ${SOURCES[1]}')

N=3
ifdamp=0
Flow('k-lrr-N3-0 k-lrra0-tmp3',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}','${TARGETS[1]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(NN)d,%(n1win)d,%(n2win)d,%(n3win)d,%(r1)g,%(r2)g,%(r3)g,%(ifdamp)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-lrr-N3','k-lrr-N3-0','put %s'%put)
Flow('k-lrr-N3-n','k-post k-lrr-N3','add scale=1,-1 ${SOURCES[1]}')

N=4
ifdamp=0
Flow('k-lrr-N4-0 k-lrra0-tmp4',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}','${TARGETS[1]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(NN)d,%(n1win)d,%(n2win)d,%(n3win)d,%(r1)g,%(r2)g,%(r3)g,%(ifdamp)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-lrr-N4','k-lrr-N4-0','put %s'%put)
Flow('k-lrr-N4-n','k-post k-lrr-N4','add scale=1,-1 ${SOURCES[1]}')

Grey('k-lrr-N1-n','title="LRR (N=1)" clip=0.05')
Grey('k-lrr-N2-n','title="LRR (N=2)" clip=0.05')
Grey('k-lrr-N3-n','title="LRR (N=3)" clip=0.05')
Grey('k-lrr-N4-n','title="LRR (N=4)" clip=0.05')


#Global
N=30
matROOT = '../matfun/'
matfun = 'FXY_MSSA'
Flow('k-grr0-N30',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-grr-N30','k-grr0-N30','put %s'%put)
Flow('k-grr-N30-n','k-post k-grr-N30','add scale=1,-1 ${SOURCES[1]}')
Grey('k-grr-N30','title="GRR (N=30)"')
Grey('k-grr-N30-n','title="GRR (N=30)" clip=0.05')

N=20
matROOT = '../matfun/'
matfun = 'FXY_MSSA'
Flow('k-grr0-N20',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-grr-N20','k-grr0-N20','put %s'%put)
Flow('k-grr-N20-n','k-post k-grr-N20','add scale=1,-1 ${SOURCES[1]}')
Grey('k-grr-N20','title="GRR (N=20)"')
Grey('k-grr-N20-n','title="GRR (N=20)" clip=0.05')

N=10
matROOT = '../matfun/'
matfun = 'FXY_MSSA'
Flow('k-grr0-N10',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-grr-N10','k-grr0-N10','put %s'%put)
Flow('k-grr-N10-n','k-post k-grr-N10','add scale=1,-1 ${SOURCES[1]}')
Grey('k-grr-N10','title="GRR (N=10)"')
Grey('k-grr-N10-n','title="GRR (N=10)" clip=0.05')

N=5
matROOT = '../matfun/'
matfun = 'FXY_MSSA'
Flow('k-grr0-N5',[os.path.join(matROOT,matfun+'.m'),'k-post'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}',%(n1)d,%(n2)d,%(n3)d,%(d1)g,%(lf)g,%(hf)g,%(N)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('k-grr-N5','k-grr0-N5','put %s'%put)
Flow('k-grr-N5-n','k-post k-grr-N5','add scale=1,-1 ${SOURCES[1]}')
Grey('k-grr-N5','title="GRR (N=5)"')
Grey('k-grr-N5-n','title="GRR (N=5)" clip=0.05')



# parameterization

# 
# 
# padx=401
# beg1=0
# v0=1400
# nv=51
# dv=25
# nx=401
# beg1=0
# x0=900
# vx0=1500
# vslope=0.67
# an=2
# 
# Flow('k-tra-vlf','k-pwd-n',
#       '''
#       pad n2=%d beg1=%d | cosft sign2=1 | put o3=0 | 
#       stolt vel=%g | 
#       vczo nv=%d dv=%g v0=%g |
#       transp plane=23 |
#       cosft sign2=-1 | 
#       window n2=%d f1=%d |
#       put o2=%g |
#       transp plane=23'''%(padx,beg1,v0,nv,dv,v0,nx,beg1,x0))
# 
# Flow('k-tra-vlfq','k-pwd-n',
#      '''
#     math output="input*input" | 
#     pad n2=%d beg1=%d | cosft sign2=1 | put o3=0 | 
#       stolt vel=%g | 
#       vczo nv=%d dv=%g v0=%g |
#       transp plane=23 |
#       cosft sign2=-1 | 
#       window n2=%d f1=%d |
#       put o2=%g |
#       transp plane=23 | clip2 lower=0
#      ''' %(padx,beg1,v0,nv,dv,v0,nx,beg1,x0))
# Flow('k-tra-sem',['k-tra-vlf','k-tra-vlfq'],
#      '''
#      mul $SOURCE |
#      divn den=${SOURCES[1]} rect1=10 rect3=10
#      ''' )
# # Flow('k-tra-pik','k-tra-sem','mutter x0=%g v0=%g half=n | scale axis=2 | pick an=%g rect1=%d rect2=%d | window' % (vx0,vslope,an,30,30))
# Flow('k-tra-pik','k-tra-sem','scale axis=2 | pick an=%g rect1=%d rect2=%d | window' % (an,30,30))
# # Flow('k-tra-pik0','k-tra-sem','mutter x0=%g v0=%g half=n | scale axis=2 | window' % (vx0,vslope*1000.0))
# 
# Flow('k-tra-slc',['k-tra-vlf','k-tra-pik'],'slice pick=${SOURCES[1]}')
# Grey1('k-tra-slc','k-tra-slc','title="Migrated Diffractions PWD"')
# Greyv('k-tra-pik','k-tra-pik','title="Velocity of PWD"')
# 

Flow('k-dmo','k-post',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')

Flow('k-pwd-dif','k-pwd-n',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')
Flow('k-pwd-ref','k-pwd-s',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')
Flow('k-ldrra-dif','k-ldrra-n',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')
Flow('k-lrra-dif','k-lrra-n',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')
Flow('k-lrra-ref','k-lrra',
        '''
        cosft sign2=1 |
        vczo v0=1500.0 nv=1 dv=-1500 |
        window |
        cosft sign2=-1
        ''')
Grey('k-dmo','title="DMO" clip=0.05')
Grey('k-pwd-dif','title="Diffraction by PWD" clip=0.05')
Grey('k-pwd-ref','title="Reflection by PWD" clip=0.05')

Grey('k-ldrra-dif','title="Diffraction by LDRRA" clip=0.05')
Grey('k-lrra-dif','title="Diffraction by LRRA" clip=0.05')
Grey('k-lrra-ref','title="Reflection by LRRA" clip=0.05')

Plot('label01',None,
        '''
        box x0=1.55 y0=6.75 label="Reflection" xt=0.5 yt=0.5 length=1.5 
        ''')
Plot('label02',None,
        '''
        box x0=3.4 y0=3.95 label="Reflection" xt=0.5 yt=0.5 length=1.8 
        ''')
Plot('label03',None,
        '''
        box x0=5.0 y0=3.9 label="Reflection" xt=0.5 yt=0.5 length=0.75 
        ''')
Plot('label001',None,
        '''
        box x0=1.55 y0=6.75 label="" xt=0.5 yt=0.5 length=1.5 
        ''')
Plot('label002',None,
        '''
        box x0=3.4 y0=3.95 label="" xt=0.5 yt=0.5 length=1.8 
        ''')
Plot('label003',None,
        '''
        box x0=5.0 y0=3.9 label="" xt=0.5 yt=0.5 length=0.75 
        ''')
Result('k-pwd-dif0','Fig/k-pwd-dif.vpl label01 label02 label03','Overlay')
Result('k-lrra-dif0','Fig/k-lrra-dif.vpl label001 label002 label003','Overlay')


### zoom 
# compare k-lrra-n and k-pwd-n
# 
## framebox A
x=50
y=6.5
w=50
w1=0.5

Flow('frame1.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        string.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame1','frame1.asc',
        '''
        dd type=complex form=native |
        graph  min1=1 max1=300 min2=5.5 max2=8.5 pad=n plotfat=15 plotcol=4 
        wantaxis=n wanttitle=n yreverse=y scalebar=n screenratio=1.2
        ''')
Flow('k-pwd-n-z1','k-pwd-n','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-lrra-n-z1','k-lrra-n','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-pwd-s-z1','k-pwd-s','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-lrra-s-z1','k-lrra','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))

## framebox B
x=125
y=6.4
w=50
w1=0.3

Flow('frame2.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        string.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame2','frame2.asc',
        '''
        dd type=complex form=native |
        graph min1=1 max1=300 min2=5.5 max2=8.5 pad=n plotfat=15 plotcol=4 
        wantaxis=n wanttitle=n yreverse=y scalebar=n screenratio=1.2
        ''')
Flow('k-pwd-n-z2','k-pwd-n','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-lrra-n-z2','k-lrra-n','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-pwd-s-z2','k-pwd-s','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-lrra-s-z2','k-lrra','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))

## framebox C
x=200
y=6.25
w=40
w1=0.25

Flow('frame3.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        string.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame3','frame3.asc',
        '''
        dd type=complex form=native |
        graph min1=1 max1=300 min2=5.5 max2=8.5 pad=n plotfat=15 plotcol=4 
        wantaxis=n wanttitle=n yreverse=y scalebar=n screenratio=1.2
        ''')
Flow('k-pwd-n-z3','k-pwd-n','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-lrra-n-z3','k-lrra-n','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-pwd-s-z3','k-pwd-s','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))
Flow('k-lrra-s-z3','k-lrra','put o2=1 d2=1 | window min2=%g max2=%g min1=%g max1=%g'%(x,x+w,y,y+w1))

## Create label A
Plot('label1',None,
        '''
        box x0=3.0 y0=6.45 label="A" xt=0.5 yt=0.5 length=0.75 
        ''')

Plot('label2',None,
        '''
        box x0=4.6 y0=6.65 label="B" xt=0.5 yt=0.5 length=0.75 
        ''')

Plot('label3',None,
        '''
        box x0=6.0 y0=7.05 label="C" xt=0.5 yt=0.5 length=0.75 
        ''')

Result('k-pwd-n0','Fig/k-pwd-n.vpl frame1 frame2 frame3 label1 label2 label3','Overlay')
Result('k-lrra-n0','Fig/k-lrra-n.vpl frame1 frame2 frame3 label1 label2 label3','Overlay')
Result('k-pwd-s0','Fig/k-pwd-s.vpl frame1 frame2 frame3 label1 label2 label3','Overlay')
Result('k-lrra-s0','Fig/k-lrra.vpl frame1 frame2 frame3 label1 label2 label3','Overlay')


Greyz('k-pwd-n-z1','title="PWD" clip=0.05')
Greyz('k-pwd-n-z2','title="PWD" clip=0.05')
Greyz('k-pwd-n-z3','title="PWD" clip=0.05')
Greyz('k-pwd-s-z1','title="PWD"')
Greyz('k-pwd-s-z2','title="PWD"')
Greyz('k-pwd-s-z3','title="PWD"')

Greyz('k-lrra-n-z1','title="LRRA" clip=0.05')
Greyz('k-lrra-n-z2','title="LRRA" clip=0.05')
Greyz('k-lrra-n-z3','title="LRRA" clip=0.05')
Greyz('k-lrra-s-z1','title="LRRA"')
Greyz('k-lrra-s-z2','title="LRRA"')
Greyz('k-lrra-s-z3','title="LRRA"')


End()