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

Fetch('chev.HH','chevron',beg.server)

Flow('chev','chev.HH',
     '''
     dd form=native |
     put o2=7923.3979 d2=12.5 o3=10298.643 d3=12.5
     label1=Time unit1=s label2=Inline unit2=m label3=Crossline unit3=m |
     scale dscale=0.01 
     ''')
Flow('spec','chev','pad n1=256 | spectra all=y')
def cubeplot(title,clip='',extra=''):
    return '''
    window n1=34 min1=0.402 |  
    transp plane=13 memsize=1000 |
    byte gainpanel=all %s |   
    grey3 frame1=225 frame2=225 frame3=25 flat=y  point1=0.8 point2=0.8 
    label3=Time unit3=s label2=Inline unit2=m label1=Crossline unit1=m
    title="%s" %s screenratio=1 labelsz=5 
    ''' % (clip,title,extra)

def Cubeplot(name,frame=20):
    t=0.402+0.006*frame
    Plot(name+'0',name,
         '''
         window n1=1 f1=%g | transp | grey label2="Inline (m) a" unit2= label1=Crossline unit1=m wanttitle=n
         grid=y g1num0=12298.4 g1num=10000 g2num0=13423.6 g2num=10000 gridcol=6 gridfat=5 labelsz=5
         screenratio=1 scalebar=y barlabel=Amplitude parallel2=n
         ''' % frame)

    Plot(name+'1',name,
         '''
         window n3=1 f3=225 | grey label2= unit2= label1=Time unit1=s labelsz=5
         grid=y g1num0=12298.4 g1num=10000 g2num0=%g g2num=10000 gridcol=6 gridfat=5 screenratio=1
         scalebar=y barlabel=Amplitude parallel2=n format2="%%3.1f" o2num=0.4 d2num=0.1 formatbar="%%g"
         title="Inline (m) b" titlesz=5 wheretitle=t dbarnum=1 nbartic=10
         ''' % t)
    Plot(name+'2',name,
         '''
         window n2=1 f2=225 | grey label2= unit2= label1=Time unit1=s labelsz=5
         grid=y g1num0=13423.6 g1num=10000 g2num0=%g g2num=10000 gridcol=6 gridfat=5 screenratio=1
         scalebar=y barlabel=Amplitude parallel2=n format2="%%3.1f" o2num=0.4 d2num=0.1 formatbar="%%g"
         title="Crossline (m) c" titlesz=5 wheretitle=t dbarnum=1 nbartic=10
         ''' % t)
    Plot(name+'12',[name+'1',name+'2'],'OverUnderAniso',vppen='txscale=2')
    Result(name,[name+'0',name+'12'],'OverUnderIso')
    #Plot(name,[name+'0',name+'12'],'OverUnderIso')

Cubeplot('chev')


#####################################
#Result('chev',cubeplot(' '))

#cubeplot('chev')

Plot('chev',cubeplot(' '))

Flow('mask','chev','math output=1 | pad beg1=25 end1=25')

Flow('dips','chev mask','pad beg1=25 end1=25 | dip rect1=5 rect2=15 rect3=15 mask=${SOURCES[1]}')
Flow('dipx','dips','window n4=1')
Flow('dipy','dips','window f4=1')

Result('dipx',
       cubeplot('Inline Dip','bar=bar1.rsf',
                '''
                color=j scalebar=y bar=bar1.rsf wanttitle=n
                barlabel="Inline Dip" bartype=h
                '''))
Result('dipy',
       cubeplot('Crossline Dip','bar=bar2.rsf',
                '''
                color=j scalebar=y bar=bar2.rsf wanttitle=n
                barlabel="Crossline Dip" bartype=h
                '''))

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

Flow('last1','chev','window f2=449 squeeze=n')
Flow('last2','chev','window f3=449  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 chev','add mode=p ${SOURCES[1]} | stack axis=1 norm=n')
Flow('corr2','ref2 chev','add mode=p ${SOURCES[1]} | stack axis=1 norm=n')

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

Flow('cos1','corr1 chev2 ref1s',
     '''
     math s1=${SOURCES[1]} s2=${SOURCES[2]} output="(s1*s2)/(input*input)"
     ''')
Flow('cos2','corr2 chev2 ref2s',
     '''
     math s1=${SOURCES[1]} s2=${SOURCES[2]} output="(s1*s2)/(input*input)"
     ''')
Flow('cos','cos1 cos2',
     'cat axis=3 ${SOURCES[1]} ')

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

Flow('seed','dips','window n2=1 n3=1 n4=1 | math output=x1')
Flow('pick','dips seed time',
     '''
     pwpaint2 seed=${SOURCES[1]} cost=${SOURCES[2]} eps=1 |
     clip2 lower=0.25 upper=0.75
     ''')

Result('pick',cubeplot('','allpos=y','color=j wanttitle=n') + ' color=j')
##############################################################################################
#Single refrence

#Flow('flat','chev-pad pick','iwarp warp=${SOURCES[1]} eps=1 n1=34 o1=0.402')
#Result('flat',cubeplot(''))
# Contour3

#Flow('wcont','chev',
     #'''
     #window n2=1 n3=1 f2=350 f3=100 | pad beg1=25 end1=35 | envelope | 
    # max1 | window n1=25 | real
    # ''')
#Plot('wcont','pick wcont',
     #'''
     #window n1=34 min1=0.402 |     
    # contour3 frame1=18 frame2=350 frame3=100 flat=y point2=0.8 point1=0.8
    # wanttitle=n wantaxis=n screenratio=1 plotfat=5 cfile=${SOURCES[1]}
    # ''')
#Result('wcont','chev wcont','Overlay')

#################################################################################################
# multiple references

picks = []
for ref in ((150,150),(200,100),(250,150),(250,250),(350,300),(300,150),(300,300)):
        pick = 'pick%d-%d' % ref
        picks.append(pick)

        Flow(pick,'dips seed cos','pwpaint3 seed=${SOURCES[1]} cost=${SOURCES[2]} ref2=%d ref3=%d'% ref)

np = len(picks)
Flow('picks',picks,
     'add ${SOURCES[1:%d]} | scale dscale=%g' % (np,1.0/np))
# Contour3

Flow('wcont','chev',
     '''
     window n2=1 n3=1 f2=225 f3=225 | pad beg1=25 end1=35 | envelope | 
     max1 | window n1=25 | real
     ''')
Plot('wcont','pick wcont',
     '''
     window n1=34 min1=0.402 | transp plane=13 |   
     contour3 frame1=225 frame2=225 frame3=20 flat=y point2=0.8 point1=0.8
     wanttitle=n wantaxis=n screenratio=1 plotfat=5 cfile=${SOURCES[1]}
     ''')
Result('wcont','chev wcont','Overlay')
#Spectral Decomposition in regular coordinate system

Flow('timefreq','chev','timefreq rect=50 verb=y nw=250 dw=0.244141')

Result('tf-curved2','timefreq',
       'window n1=1 f1=27 n4=1 f4=50 | grey allpos=y color=j title=Curved')
Result('slice','chev','window n1=1 f1=24 | grey transp=n title=Amplitude unit1=m unit2=m screenratio=1 scalebar=y parallel2=n labelsz=6')

for freq in (10,15,20,25,30,35,40,45,50,60):
    slice = 'slice%d' % freq
    Flow(slice,'timefreq','window n1=1 f1=25 n2=1 min2=%d' % freq)
    Result(slice,' grey title="%g Hz" bias=-0.01 transp=n unit1=m unit2=m allpos=y scalebar=y parallel2=n labelsz=6' % freq)


######################################################################################################### t0

# Convert time to pseudo-depth

dt=0.006
t0=0.252
dz=20
dzdt = dz/dt
dtdz=dt/dz
z0 = t0*dzdt

#Flow('pickz','pick','put d1=%g o1=%g | scale dscale=%g | pad beg1=100' % (dz,z0,dzdt))
Flow('pickz','pick','put d1=%g o1=%g | scale dscale=%g' % (dz,z0,dzdt))
#Flow('slowreal1','pickz','lineiko what=s | window n1=84 min1=%g' % z0)
Flow('slowreal1','pickz','lineiko what=s')
Result('slowreal1','put d1=%g o1=%g | ' % (dt,t0) + cubeplot('','allpos=y','color=j'))

Flow('t0real1','slowreal1','eikonal vel=n plane3=y plane2=y zshot=%g | clip2 upper=5000' % z0)

Result('t0real1','put d1=%g o1=%g | ' % (dt,t0) + cubeplot('','allpos=y','color=j'))

Plot('t0real1',
      '''
      put d1=%g o1=%g | window n1=34 min1=0.402 | transp plane=13 memsize=1000 |
      contour3 frame1=225 frame2=225 frame3=25 flat=y point2=0.8 point1=0.8
      wanttitle=n wantaxis=n screenratio=1 plotfat=5 plotcol=4 
      ''' % (dt,t0))

Result('chev-t0real1','chev t0real1','Overlay')


############################################################################################################ x0

# Inline direction
Flow('t0real','t0real1','put d1=%g o1=%g | scale dscale=%g' % (dt,t0,dtdz))
Result('t0real',cubeplot('','allpos=y','color=j'))

Flow('distreal1','t0real1','math output=x2')

Flow('zeroreal1','t0real1','math output=0')

Flow('x0real1','t0real1 distreal1 zeroreal1','lineiko what=i time=${SOURCES[1]} slow=${SOURCES[2]}')

Result('x0real1',cubeplot('','allpos=y','color=j'))

Plot('x0real1',
      '''
      put d1=%g o1=%g | window n1=34 min1=0.402 | transp plane=13 memsize=1000 |      
      contour3 frame1=225 frame2=225 frame3=20 flat=y point2=0.8 point1=0.8
      wanttitle=n wantaxis=n screenratio=1 plotfat=5 plotcol=5 
      ''' % (dt,t0))

Result('x0t01','chev t0real1 x0real1','Overlay')

############################################################################################################# y0

# Crossline direction

Flow('distreal2','t0real1','math output=x3')

Flow('y0real2','t0real1 distreal2 zeroreal1','lineiko what=i time=${SOURCES[1]} slow=${SOURCES[2]}')

Result('y0real2','put d1=%g o1=%g | ' % (dt,t0) + cubeplot('','allpos=n mean=y','color=j'))
Plot('y0real2',
      '''
      put d1=%g o1=%g | window n1=34 min1=0.402 | transp plane=13 memsize=1000 |      
      contour3 frame1=225 frame2=225 frame3=20 flat=y point2=0.8 point1=0.8
      wanttitle=n wantaxis=n screenratio=1 plotfat=5 plotcol=6 
      ''' % (dt,t0))

Result('coord','chev t0real1 x0real1 y0real2','Overlay')

def Cubeplotx0(name,frame=20):
    t=0.402+0.006*frame
    Plot(name+'0',name,
         '''
         put d1=%g o1=%g | window n1=1 f1=%g | transp | contour label2=Inline unit2=m label1=Crossline unit1=m wanttitle=n
         grid=y g1num0=12298.4 g1num=10000 g2num0=13423.6 g2num=10000 gridcol=6 gridfat=5 labelsz=5
         screenratio=1 scalebar=y barlabel=Amplitude parallel2=n plotcol=4 
         ''' % (dt,t0,frame))

    Plot(name+'1',name,
         '''
         put d1=%g o1=%g | window n3=1 f3=225 | contour label2= unit2= label1=Time unit1=s labelsz=5
         grid=y g1num0=12298.4 g1num=10000 g2num0=%g g2num=10000 gridcol=6 gridfat=5 screenratio=1
         scalebar=y barlabel=Amplitude parallel2=n format2="%%3.1f" o2num=0.4 d2num=0.1 formatbar="%%g"
         title="Inline (m)" titlesz=5 wheretitle=t dbarnum=1 nbartic=10 plotcol=4 
         ''' % (dt,t0,t))
    Plot(name+'2',name,
         '''
         put d1=%g o1=%g | window n2=1 f2=225 | contour label2= unit2= label1=Time unit1=s labelsz=5
         grid=y g1num0=13423.6 g1num=10000 g2num0=%g g2num=10000 gridcol=6 gridfat=5 screenratio=1
         scalebar=y barlabel=Amplitude parallel2=n format2="%%3.1f" o2num=0.4 d2num=0.1 formatbar="%%g"
         title="Crossline (m)" titlesz=5 wheretitle=t dbarnum=1 nbartic=10 plotcol=4 
         ''' % (dt,t0,t))
    Plot(name+'12',[name+'1',name+'2'],'OverUnderAniso',vppen='txscale=2')
    Result(name,[name+'0',name+'12'],'OverUnderIso')
    Plot(name,[name+'0',name+'12'],'OverUnderIso')

#Cubeplotx0('t0real1')
#Cubeplotx0('x0real1')
#Cubeplotx0('y0real2')
################################################################################mapping from (t,x,y) to (t0,x0,y0)
Flow('x0real','x0real1','put d1=%g o1=%g | scale dscale=%g' % (dt,t0,dtdz))
Flow('y0real','y0real2','put d1=%g o1=%g | scale dscale=%g' % (dt,t0,dtdz))
Result('x0real',cubeplot('','allpos=y','color=j'))
Result('y0real',cubeplot('','allpos=y','color=j'))
Flow('warpreal1','t0real1 x0real1 y0real2',
      '''
      add add=%g |
      cat axis=4 ${SOURCES[1:3]} 
      ''' % z0)

Flow('chevp','chev','pad beg1=25 end1=25')
Flow('chevvp','chevp warpreal','iwarp3 warp=${SOURCES[1]} eps=1')

Flow('warpreal','warpreal1','put d1=%g o1=%g | scale dscale=%g' % (dt,t0,dtdz))
Flow('chevz','chev','pad beg1=25 end1=25 | put d1=%g o1=%g | scale dscale=%g ' % (dz,z0,dzdt))

Flow('chevvz','chevz warpreal1','iwarp3 warp=${SOURCES[1]} eps=1')

Flow('chevv','chevvz','put d1=%g o1=%g | scale dscale=%g' % (dt,t0,dtdz))
#Flow('chevvs','chevv','window n1=34 min1=0.402')
Result('chevv',cubeplot('Warped'))
#Result('chevv',cubeplot('','allpos=y','color=j'))


Flow('flats','chev pick','pad beg1=25 end1=25 | iwarp warp=${SOURCES[1]} eps=1 n1=34 o1=0.402')
Flow('chev-hor','flats pick','iwarp warp=${SOURCES[1]} eps=1 n1=34 o1=0.40 inv=n')
Result('chev-hor',cubeplot('Warped'))
Result('flats',cubeplot('Warped'))
#Spectral decomposition in horizon coordinates

Flow('timefreq_2','flats','timefreq rect=5 nw=24 dw=5')
Flow('timefreq2','flats','timefreq rect=50 verb=y nw=250 dw=0.244141')

Result('tf-curved1','timefreq2',
       'window n1=1 f1=27 n4=1 f4=50 | grey allpos=y color=j title=Curved')
#for freq in (10,15,20,25,30,35,40,45,50,60):
   # slice = 'slice-4%d' % freq
    #Flow(slice,'timefreq_2','window n1=1 f1=25 n2=1 min2=%d  | sharpen perc= 40' % freq)
    #Result(slice,' grey title="%g Hz" bias=0.1 transp=n unit1=ft unit2=ft allpos=y scalebar=y parallel2=n labelsz=6' % freq)

#for freq in (10,15,20,25,30,35,40,45,50,60):
   # slice = 'slice-6%d' % freq
    #Flow(slice,'timefreq_2 pick',
       #  '''
        # window n2=1 min2=%d  | 
        # iwarp warp=${SOURCES[1]} eps=1 n1=34 o1=0.40 inv=n 
        # ''' % freq)
    #Result(slice,' window n1=1 f1=45 | grey title="%g Hz" bias=-0.01 transp=n unit1=ft unit2=ft allpos=y scalebar=y parallel2=n labelsz=6' % freq)


Flow('spec_chev','chev','pad n1=256 | spectra all=y')

Flow('spec_flat','flats','pad n1=256 | spectra all=y')
Flow('spec_strat','chevv','pad n1=256 | spectra all=y')

Result('spec','spec_strat spec_chev','cat axis=2 ${SOURCES[1:2]} | scale axis=1 | graph dash=2,0 title=Spectrum label2= unit2=')
#####################################

#Spectral Decomposition in stratigraphic coordinate system

Flow('timefreq_1','chevv','timefreq rect=5 nw=24 dw=5')

for freq in (10,15,20,25,30,35,40,45,50,60):
    slice = 'slice-1%d' % freq
    Flow(slice,'timefreq_1','window n1=1 f1=45 n2=1 min2=%d  | sharpen perc= 60' % freq)
    Result(slice,' grey title="%g Hz" transp=n unit1=m unit2=m allpos=y parallel2=n labelsz=6 labelfat=4 font=2 scalebar=y barlabel=Amplitude' % freq)


#############################################
#Spectral Decomposition in stratigraphic coordinate system shown in cartesian
#for freq in (10,15,20,25,30,35,40,45,50,60):
   # slice = 'slice-2%d' % freq
   # Flow(slice,'timefreq_1 warpreal',
        # '''
        # window n2=1 min2=%d |
        # iwarp3 warp=${SOURCES[1]} inv=n |
         #sharpen perc=40
        # ''' % freq)
    #Result(slice,' window n1=1 f1=45 | grey title="%g Hz" bias=0 transp=n unit1=ft unit2=ft allpos=y scalebar=y parallel2=n labelsz=6' % freq)


#Flow('timefreqz','timefreq','pad beg1=25 end1=25 | put d1=%g o1=%g | scale dscale=%g' % (dz,z0,dzdt))

#Spectral decomposition in Cartesian Coordinates shown in stratigraphic coordinates
#Flow('timefreqd','timefreq','window n2=1 f2=30 | pad beg1=25 end1=25 | put d1=%g o1=%g | scale dscale=%g ' % (dz,z0,dzdt))


for freq in (10,15,20,25,30,35,40,45,50,60):
    slice = 'slice-3%d' % freq
    timefreqd = 'timefreqd%d' %freq
    Flow(timefreqd,'timefreq','window n2=1 min2=%d | pad beg1=25 end1=25 | put d1=%g o1=%g | scale dscale=%g ' % (freq,dz,z0,dzdt))
    Flow(slice,[timefreqd ,'warpreal1'],'iwarp3 warp=${SOURCES[1]} eps=1')
    Result(slice,
           '''
            put d1=%g o1=%g | scale dscale=%g |
            window n1=1 f1=45 | 
            grey title="%g Hz" bias=-0.1 transp=n unit1=m unit2=m allpos=y scalebar=y parallel2=n labelsz=6 labelfat=4 font=2 barlabel=Amplitude
           ''' % (dt,t0,dtdz,freq))



###########################################################################Mapping back from (t0,x0,y0) to (t,x,y)

Flow ('chevvvz','chevvz warpreal1','iwarp3 warp=${SOURCES[1]} inv=n')
Flow('chevvv','chevvvz','put d1=%g o1=%g | scale dscale=%g' % (dt,t0,dtdz))
Result('chevvv',cubeplot('Warped back'))



#Flow('den2','timefreq','stack')
#Flow('dom2','timefreq den','math output="x2*input" | stack | div ${SOURCES[1]}')
#Flow('den','timefreq_2','stack')
#Flow('dom','timefreq_2 den','math output="x2*input" | stack | div ${SOURCES[1]}')
#Flow('den1','timefreq_1','stack')
#Flow('dom1','timefreq_1 den1','math output="x2*input" | stack | div ${SOURCES[1]}')
#Result('dom',cubeplot(' '))
##Result('domm1','dom1','window n2=1 f2=225 | grey color=j scalebar=y allpos=y clip=70 bias=5 minval=5 maxval=180 title="Dominant Frequency" barlabel=Frequency barunit=Hz')
#Result('domm2','dom2','window n2=1 f2=225 | grey color=j scalebar=y allpos=y clip=70 bias=5 minval=5 maxval=180 title="Dominant Frequency" barlabel=Frequency barunit=Hz')
#Result('ddom','dom2','window n1=1 f1=25 | grey color=j scalebar=y allpos=y clip=70 bias=40 minval=5 maxval=180 title="Dominant Frequency" barlabel=Frequency barunit=Hz')
#Result('dom1',cubeplot(' '))
#Result('dom11','dom1 warpreal',
      # '''
       #iwarp3 warp=${SOURCES[1]} inv=n  | 
      # window n1=1 f1=50 |
       #grey color=j scalebar=y allpos=y clip=70 
       #bias=40 minval=5 maxval=180 title="Dominant Frequency" barlabel=Frequency barunit=Hz
      # ''')
End()