from rsf.proj import *
from rsf.recipes.beg import server

def Grey(data,other):
        Result(data,
           '''
                        grey transp=y color=j 
                        unit1=s unit2= allpos=y scalebar=y transp=n
                        parallel2=n labelsz=6 clip=0.5 label1=Inline label2=Crossline %s'''%(other))#screenratio=1

def Grey33f(data,other):
        Result(data,
           '''
                        put d1=0.002 o1=0 | scale axis=3| byte allpos=y bar=bar.rsf| grey3 flat=n transp=y allpos=y
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  
                        parallel2=n labelsz=6 screenratio=1 point1=0.8 point2=0.8 minval=0 maxval=1 %s'''%(other))

sgy='BOONS.SGY'
#datapath: /media/CHENYK2/Seismicdata/boonsville
#Fetch(sgy,'boonsegy',server)

#Flow('boons tboons boons.asc boons.bin',
#     sgy,
#     'segyread tfile=${TARGETS[1]} hfile=${TARGETS[2]} bfile=${TARGETS[3]}')

#C 3 Line number in bytes 9  - 12      105 - 201 (fldr)
#C 4 Trace number in bytes 21 - 24      74 - 206 (cdp)                       

# 110 foot stacking bins

Flow('boon3','../boonsville/boon3.rsf','cp')

Result('boon3','byte gainpanel=all | grey3 frame1=1000 frame2=76 frame3=45 title="Boonsville Data" flat=y point1=0.7 point2=0.7')

Flow('boont','boon3','window j1=2')

Flow('dip3','boont','fdip rect1=20 rect2=5 rect3=5 verb=y')

Flow('dip1','dip3','window n4=1')
Flow('dip2','dip3','window f4=1')

Result('dip1','byte gainpanel=all bar=bar1.rsf | grey3 frame1=500 frame2=76 frame3=45 title="Inline Dip" color=j scalebar=y flat=y point1=0.7 point2=0.7')
Result('dip2','byte gainpanel=all bar=bar2.rsf | grey3 frame1=500 frame2=76 frame3=45 title="Crossline Dip" color=j scalebar=y flat=y point1=0.7 point2=0.7')

Flow('seed','dip3','window n2=1 n3=1 n4=1 | math output=x1')

Flow('shift1','boont','window f2=1')
Flow('shift2','boont','window f3=1')

Flow('last1','boont','window f2=-1 squeeze=n')
Flow('last2','boont','window f3=-1 squeeze=n')

Flow('ref1','shift1 last1','cat axis=2 ${SOURCES[1]}')
Flow('ref2','shift2 last2','cat axis=3 ${SOURCES[1]}')

Flow('ref1s','ref1','add mode=p $SOURCE | stack axis=1 norm=n')
Flow('ref2s','ref2','add mode=p $SOURCE | stack axis=1 norm=n')

Flow('corr1','ref1 boont','add mode=p ${SOURCES[1]} | stack axis=1 norm=n')
Flow('corr2','ref2 boont','add mode=p ${SOURCES[1]} | stack axis=1 norm=n')

Flow('boon2','boont','add mode=p $SOURCE | stack axis=1 norm=n')

Flow('wcos1','corr1 boon2 ref1s',
     '''
     math s1=${SOURCES[1]} s2=${SOURCES[2]} output="(s1*s2)/(input*input)"
     ''')
Flow('wcos2','corr2 boon2 ref2s',
     '''
     math s1=${SOURCES[1]} s2=${SOURCES[2]} output="(s1*s2)/(input*input)"
     ''')
Flow('wcos','wcos1 wcos2',
     '''
     cat axis=3 ${SOURCES[1]} |
     smooth rect1=40 rect2=40 
     ''')

Flow('time','wcos',
     '''
     mul $SOURCE | stack axis=3 norm=n |
     put o1=0 d1=1 o2=0 d2=1 o3=0 d3=1 |
     eikonal vel=n zshot=76 yshot=45
     ''')

Flow('pick','dip3 seed time',
     'pwpaint2 seed=${SOURCES[1]} cost=${SOURCES[2]} order=1')
Result('pick','byte gainpanel=all allpos=y | grey3 frame1=500 frame2=76 frame3=45 title="Relative Age" color=j flat=y point1=0.7 point2=0.7')

Flow('flat','boont pick','iwarp warp=${SOURCES[1]}')
Result('boont','byte gainpanel=all | grey3 frame1=500 frame2=76 frame3=45 wanttitle=n flat=y')
Result('flat','byte gainpanel=all | grey3 frame1=500 frame2=76 frame3=45 wanttitle=n flat=y')

Flow('wboon','boont','window min1=0 max1=1.598')
Flow('wflat','flat', 'window min1=0 max1=1.598')

Result('wboon','byte gainpanel=all | grey3 frame1=119 frame2=76 frame3=45 wanttitle=n flat=y color=G')
Result('wflat','byte gainpanel=all | grey3 frame1=119 frame2=76 frame3=45 wanttitle=n flat=y color=G')

# Contour3

Flow('cont','boont',
     '''
     window n2=1 n3=1 f2=76 f3=45 | envelope | 
     max1 | window n1=10 | real
     ''')
Plot('wboon','put o2=0 d2=1 o3=0 d3=1 label2=Inline label3=Crossline |scale axis=3 | byte gainpanel=all bar=bar.rsf | grey3 flat=y transp=y unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=0.6 title="Strata Interpretation" minval=-1 maxval=1 scalebar=y')
Plot('wboon-2','wboon','math output="abs(input)" | put o2=0 d2=1 o3=0 d3=1 label2=Inline label3=Crossline | scale axis=3  | byte gainpanel=all  allpos=y bar=bar.rsf| grey3 flat=y transp=y unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=0.6 title="Strata Interpretation" color=j minval=0 maxval=1 scalebar=y')

Plot('cont','pick cont',
     '''
     window min1=0 max1=1.598 |
     contour3 flat=y transp=y plotfat=5
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  
                        parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=0.6 cfile=${SOURCES[1]}
     ''')
Result('cont','wboon cont','Overlay')

# Strongest horizon

Flow('cont0','boont',
     '''
     window n2=1 n3=1 f2=76 f3=45 | envelope | 
     max1 | window n1=1 | real | math output="1.27" 
     ''')
Flow('cont-c','boont',
     '''
     window n2=1 n3=1 f2=76 f3=45 | envelope | 
     max1 | window n1=1 | real | math output="0.86" 
     ''')
Flow('cont-v','boont',
     '''
     window n2=1 n3=1 f2=76 f3=45 | envelope | 
     max1 | window n1=1 | real | math output="1.05" 
     ''')
Plot('cont0','pick cont0',
     '''
     put o2=0 d2=1 o3=0 d3=1 label2=Inline label3=Crossline |window min1=0 max1=1.598 |
     contour3 flat=y transp=y plotfat=5
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  
                        parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=0.6 cfile=${SOURCES[1]} scalebar=y title= 
     ''')
Plot('cont-c','pick cont-c',
     '''
     put o2=0 d2=1 o3=0 d3=1 label2=Inline label3=Crossline |window min1=0 max1=1.598 |
     contour3 flat=y transp=y plotfat=5
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  plotcol=4 
                        parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=0.6 cfile=${SOURCES[1]} scalebar=y title= 
     ''')
Plot('cont-v','pick cont-v',
     '''
     put o2=0 d2=1 o3=0 d3=1 label2=Inline label3=Crossline |window min1=0 max1=1.598 |
     contour3 flat=y transp=y plotfat=1
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  plotcol=3
                        parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=0.6 cfile=${SOURCES[1]} scalebar=y title= 
     ''')
Result('stra-int','wboon cont0 cont-c cont-v','Overlay')


Result('stra-int-2','wboon-2 cont0 cont-c cont-v','Overlay')

Flow('tfsswt2-2d3d','../boonsville/tfsswt2-2d3d.rsf','cp')
Grey33f('tfsswt2-2d3d','title="Strata Interpretation" flat=y frame3=75 frame1=525 frame2=80 color=j screenratio=0.5 label1=Time unit1=s label2=Inline label3=Crossline point1=0.5 point2=0.5 screenratio=0.6 scalebar=y')

Result('stra-tf','Fig/tfsswt2-2d3d.vpl cont0 cont-c cont-v','Overlay')


Flow('stra1','tfsswt2-2d3d','scale axis=3 | window n1=1 f1=430')
Flow('stra2','tfsswt2-2d3d','scale axis=3 | window n1=1 f1=480')
Flow('stra3','tfsswt2-2d3d','scale axis=3 | window n1=1 f1=525')
Flow('stra4','tfsswt2-2d3d','scale axis=3 | window n1=1 f1=565')

Grey('stra1','title="Top of Caddo" unit1= minval=0 maxval=0.85')
Grey('stra2','title="Strata between Caddo and Vineyard" unit1= minval=0 maxval=0.85')
Grey('stra3','title="Top of Vineyard" unit1= minval=0 maxval=0.85')
Grey('stra4','title="Strata between Vineyard and Base" unit1= minval=0 maxval=0.85')




# cont1 and cont2 -> mask

t0=1.129
dt=0.1

t1=t0-dt
t2=t0+dt

t1=0
t2=1.598

Flow('mask','pick','mask min=%g max=%g | dd type=float | smooth rect1=10' % (t1,t2))
Flow('boonm','boont mask','mul ${SOURCES[1]}')

Plot('wboonm','boonm',
     '''
     window min1=0 max1=1.598 |
     byte gainpanel=all | 
     grey3 flat=y transp=y
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  
                        parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=1
     ''')
Plot('contm','pick cont0',
     '''
     window min1=0 max1=1.598 |
     contour3 flat=y transp=y plotfat=5
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  
                        parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=1 cfile=${SOURCES[1]}
     ''')

Plot('boonm',
     '''
     window min1=0 max1=1.598 |
     byte gainpanel=all |
     grey3 flat=y transp=y
                        unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  
                        parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=1
     ''' )

Result('boonm','wboonm contm','Overlay')

# previous method - similarity with envelope
#Flow('rot','boonm','simenv verb=n rect=1000',split=[3,97],reduce='cat axis=4')

# new method - local skewness
# Flow('sqr','boonm','localskew verb=n rect=1000',split=[3,97],reduce='cat axis=4')
# Flow('kur','boonm','localskew verb=n rect=1000 const=y inv=n',split=[3,97],reduce='cat axis=4')
# Flow('skew','sqr kur','mul ${SOURCES[1]}')
# 
# Flow('mrot','sqr','stack axis=1 | transp',split=[4,97],reduce='cat axis=3')
# 
# Flow('mpik','mrot','pick rect1=10 rect2=10 | add add=90')
# 
# Result('mpik',
#        '''
#        grey color=j scalebar=y title="Estimated Phase Rotation"
#        unit1= unit2= transp=n yreverse=n 
#        ''')

# Flow('ampl','flat','window n1=1 min1=1.129')
# Result('ampl','grey title="Amplitude Before Rotation" unit1= unit2= transp=n yreverse=n color=G')
# 
# Flow('npik','mpik','spray axis=1 n=1001 d=0.002')
# 
# Flow('boonr','boonm npik','phaserot phase=${SOURCES[1]}')
# 
# Result('boonr',
#        '''
#        window min1=0 max1=1.598 |
#        byte gainpanel=all |
#        grey3 flat=y transp=y
#                       unit1=s unit2= allpos=y  frame3=75 frame1=525 frame2=80 unit1=s unit2= unit3= framelabelcol=VP_BLACK  
#                       parallel2=n labelsz=6 point1=0.5 point2=0.5 screenratio=1
#        ''' )


End()