from rsf.proj import *
import rsf.recipes.stack as stack
from rsf.prog import RSFROOT

name='bei'
v0=1.4
nv=51
dv=0.025
x0=7.705
dx=0.0335
nx=250
nt=751
padx=521
an=2
tmin=0
tmax=3.996
rect1=40
rect2=15
srect1=1
srect2=5
frect1=40
frect2=10
units='m'
beg1=0
f1=15
j3=1
vslope=0.67
vx0=1.5
nout=2048


def Grey(data,data0,other):
    Result(data,data0,'''
    grey %s
    label1=Time unit1=s label2=Distance clip=10000000
    ''' % (other))
def Greyv(data,data0,other):
    Result(data,data0,
        '''
        grey 
        color=j scalebar=y label1=Time unit1=s label2=Distance unit2=km bias=%g barlabel=Velocity barunit="km/s"
        barreverse=y %s
        ''' % (v0+0.5*nv*dv,other))

Fetch('midpts.hh','midpts')

Flow('bei','midpts.hh','dd form=native | put d2=0.067 o2=0.132')

scn = name+'-scn'
vel = name+'-vel'

Flow(scn,name,
         'mutter v0=%g | vscan semblance=y v0=%g nv=%d dv=%g' % (v0,v0,nv,dv))

if vslope:
        pick = 'mutter x0=%g v0=%g half=n | ' % (vx0,vslope)
else:
        pick = ''

pick = pick + 'pick rect1=%d rect2=%d | window' % (rect1,rect2)

def grey(title):
        return '''
        window n1=%d |
        grey title="%s" 
        label1=Time unit1=s label2=Distance unit2="%s"
        ''' % (nt,title,units)

def velgrey(title):
        return grey(title) + '''
        color=j scalebar=y mean=y barlabel=Velocity barunit="%s/s"
        barreverse=y
        ''' % units

Flow(vel,scn,pick)
Result(vel,velgrey('RMS Velocity'))

nmo = name+'-nmo'
stk = name+'-stk'

Flow(nmo,[name,vel],'mutter v0=%g | nmo velocity=${SOURCES[1]}' % v0)
Flow(stk,nmo,'stack')
Result(stk,grey('NMO Stack'))

stk2=stk+'2'
Flow(stk2,nmo,
         '''
         window f1=%d | logstretch nout=%d |
         fft1 | transp plane=13 memsize=500 |
         finstack | transp memsize=500 |
         fft1 inv=y | window n1=%d |
         logstretch inv=y | pad beg1=%d |
         transp | bandpass fhi=%g | transp
         ''' % (f1,nout,nout,f1,0.75*0.5/dx))
Result(stk2,grey('DMO Stack'))


if frect1==0:
        frect1 = 2*rect1

if frect2==0:
        frect2 = 2*rect2/j3

def grey(title):
        return '''
        window n1=%d |
        grey title="%s" 
        label1=Time unit1=s label2=Distance unit2="%s"
        ''' % (nt,title,units)

def velgrey(title):
        return grey(title) + '''
        color=j scalebar=y bias=%g barlabel=Velocity barunit="%s/s"
        barreverse=y
        ''' % (v0+0.5*nv*dv,units)

dip = 'bei-dip'
Flow('bei-dip','bei-stk2','dip rect1=%d rect2=%d' % (rect1,rect2))
Result('bei-dip',grey('Dominant Slope') + \
       ' color=j scalebar=y barlabel=Slope barunit=samples ')


Flow('bei-pwk','bei-dip',
     'noise rep=y seed=2007 | pwdsmooth2 dip=$SOURCE rect1=3 rect2=%d' %  srect2)
Result('bei-pwk',grey('Pattern'))

Flow('bei-pwd',['bei-stk2',dip],'pwd dip=${SOURCES[1]}')
# Flow('bei-pwd.bin','bei-pwd','rsf2bin bfile=$TARGET')
Result('bei-pwd',grey('Separated Diffractions'))


Flow('bei-ref',['bei-stk2',dip],
     'pwdsmooth2 dip=${SOURCES[1]} rect1=%d rect2=%d' %  (srect1,srect2))
Result('bei-ref',grey('Separated Reflections'))


Flow('bei-shp',['bei-stk2','bei-ref'],'add ${SOURCES[1]} scale=1,-1')
Flow('bei-shp.bin','bei-shp','rsf2bin bfile=$TARGET')

Result('bei-shp',grey('Diffractions'))
# 
# velcon = '''
# 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)
# 
# vlf='bei-vlf'
# Flow('bei-vlf','bei-pwd',velcon)
# Flow('bei-vlfq','bei-pwd',
#      '''
#      math output="input*input" | %s | clip2 lower=0
#      ''' % velcon)
# 
# if j3 > 1:
#     focus = 'window j3=%d | ' % j3
# else:
#     focus = ''
#     
# focus = focus + '''    
# focus rect1=%d rect3=%d |
# math output="1/abs(input)" |
# cut max1=%g | cut min1=%g 
# ''' % (frect1,frect2,tmin,tmax)
# 
# foc='bei-foc'
# Flow(foc,'bei-vlf',focus)
# 
# sem='bei-sem'
# Flow('bei-sem',['bei-vlf','bei-vlfq'],
#      '''
#      mul $SOURCE |
#      divn den=${SOURCES[1]} rect1=%d rect3=%d
#      ''' % (rect1,rect2))
# 
# pik='bei-pik'
# 
# if vslope:
#     pick = 'mutter x0=%g v0=%g half=n | ' % (vx0,vslope)
# else:
#     pick = ''
# 
# pick = pick + 'scale axis=2 | pick an=%g rect1=%d rect2=%d | window' % (an,rect1,rect2)
# 
# 
# if j3 > 1:
#     pick2 = pick + ''' |
#     transp |
#     spline n1=%d d1=%g o1=%g |
#     transp
#     ''' % (nx,dx,x0)
# else:
#     pick2 = pick
#     
# Flow('bei-pik','bei-sem',pick2)
# Result('bei-pik',velgrey('Migration Velocity'))
# 
# slc='bei-slc'
# Flow('bei-slc',['bei-vlf','bei-pik'],'slice pick=${SOURCES[1]}')
# 
# 
# Flow('bei-vlf2','bei-stk2',velcon)
# Flow('bei-slc2',['bei-vlf2','bei-pik'],'slice pick=${SOURCES[1]}')
# 
# 
# 
# Flow('bei-slc1','bei-slc2','shapeagc rect1=200')
# 
# 
# Flow('bei-vlf0','bei-ref',velcon)
# Flow('bei-slc0',['bei-vlf0','bei-pik'],'slice pick=${SOURCES[1]}')
# 
# Result('bei-slc',grey('Migrated Diffractions'))
# Result('bei-slc2',grey('Migrated Stack'))
# Result('bei-slc1',grey('Migrated Stack'))
# Result('bei-slc0',grey('Migrated Reflections'))
# 
# Flow('bei-slc3','bei-slc0','agc rect1=200')
# 
# Result('bei-slc3',grey('Migrated Reflections'))
#  
#     
# panels = []
# for cig in (41,83,110,137,165):
#     panel = 'cig%d' % cig
#     panels.append(panel)
# 
#     Plot(panel+'g','bei-sem',
#          '''
#          window min1=1 max1=3 n3=1 f3=%d | grey allpos=y
#          label1=Time unit1=s label2=Velocity unit2=km/s
#          labelsz=15 titlesz=18 color=j title="IG at %d km"
#          ''' % (cig,int(7.705+0.0335*cig)))
# 
#     pik = 'pik%d' % cig
#     Plot(pik,'bei-pik',
#          '''
#          window min1=1 max1=3 n2=1 f2=%d |
#          graph transp=y yreverse=y pad=n
#          wanttitle=n labelsz=15 titlesz=18 wantaxis=n
#          min2=1.425 max2=2.6 plotcol=7 plotfat=10
#          ''' % cig)
# 
#     Plot(panel,[panel+'g',pik],'Overlay')
#     
# Result('panel',panels,'SideBySideAniso',vppen='txscale=1.2')
# 
# Flow('pred','bei-dip','pwpaint i0=125 eps=0.1')
# Plot('stk','bei-stk2','grey title=Intepretation label2=Distance unit2=km label1=Time unit1=s')
# Plot('pred','contour wanttitle=n wantaxis=n plotfat=5')
# Result('pred','stk pred','Overlay')




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

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

############################################################
## with parameter
############################################################
Flow('bei-rr-t',[os.path.join(matROOT,matfun+'.m'),'bei-stk2'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}');quit"
     '''%vars(),stdin=0,stdout=-1)

# Flow('bei-rr','bei-rr.bin','bin2rsf bfile=${SOURCES[0]} n1=1000 n2=250 d1=0.004 d2=0.0334 o1=0 o2=7.705')
Flow('bei-rr','bei-rr-t','put n1=1000 n2=250 d1=0.004 d2=0.0334 o1=0 o2=7.705')
Flow('bei-rrf','bei-stk2 bei-rr','add scale=1,-1 ${SOURCES[1]}')

Grey('bei-rr','bei-rr','title="Separated reflections by LRRA"')
Grey('bei-rrf','bei-rrf','title="Separated diffractions by LRRA"')
Grey('bei-tra','bei-shp','title="Separated diffractions by PWD"')
Grey('bei-in','bei-stk2','title="Input seismic data"')



### Imaging for RR
Flow('bei-rr-vlf','bei-rrf',
        '''
        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('bei-rr-vlfq','bei-rrf',
     '''
    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('bei-rr-sem',['bei-rr-vlf','bei-rr-vlfq'],
     '''
     mul $SOURCE |
     divn den=${SOURCES[1]} rect1=%d rect3=%d
     ''' % (rect1,rect2))
Flow('bei-rr-pik','bei-rr-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('bei-rr-slc',['bei-rr-vlf','bei-rr-pik'],'slice pick=${SOURCES[1]}')
Grey('bei-rr-slc','bei-rr-slc','title="Migrated Diffractions LRRA"')
Greyv('bei-rr-pik','bei-rr-pik','title="Velocity of LRRA"')


### Imaging for PWD
Flow('bei-tra-vlf','bei-shp',
        '''
        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('bei-tra-vlfq','bei-shp',
     '''
    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('bei-tra-sem',['bei-tra-vlf','bei-tra-vlfq'],
     '''
     mul $SOURCE |
     divn den=${SOURCES[1]} rect1=%d rect3=%d
     ''' % (rect1,rect2))
Flow('bei-tra-pik','bei-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('bei-tra-slc',['bei-tra-vlf','bei-tra-pik'],'slice pick=${SOURCES[1]}')
Grey('bei-tra-slc','bei-tra-slc','title="Migrated Diffractions PWD"')
Greyv('bei-tra-pik','bei-tra-pik','title="Velocity of PWD"')


## framebox A
x=9.3
y=2.8
w=3.7
w1=1.0

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=7.705 max1=16.0465 min2=0 max2=4.0 pad=n plotfat=15 plotcol=4 
        wantaxis=n wanttitle=n yreverse=y scalebar=n
        ''')
## framebox B
x=10.0
y=1.25
w=2.5
w1=0.7

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=7.705 max1=16.0465 min2=0 max2=4.0 pad=n plotfat=15 plotcol=3 
        wantaxis=n wanttitle=n yreverse=y scalebar=n
        ''')

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

Plot('label2',None,
        '''
        box x0=10.2 y0=5.6 label="" xt=0.5 yt=0.5 length=0.75 
        ''')

Result('bei-tra-slc-0','Fig/bei-tra-slc.vpl frame1 frame2 label1 label2','Overlay')
Result('bei-rr-slc-0','Fig/bei-rr-slc.vpl frame1 frame2 label1 label2','Overlay')
Result('bei-tra-0','Fig/bei-tra.vpl frame1 frame2 label1 label2','Overlay')
Result('bei-rr-0','Fig/bei-rrf.vpl frame1 frame2 label1 label2','Overlay')



End()