from rsf.proj import *

############################
# Data fetching, SEG-Y -> RSF
#############################

rawsegy=['L23535','L23536','L23537']

for file in rawsegy  :
    Fetch(file+'.SGY',
          server='http://certmapper.cr.usgs.gov',
          top='nersl/NPRA/SEISMIC/1981/31_81',
          dir='DMUX')
    Flow([file, file+'.bin',  file+'.asc', 't'+file],
         file+'.SGY',
         '''
         segyread bfile=${TARGETS[1]} hfile=${TARGETS[2]} tfile=${TARGETS[3]}
         ''')

# concatinate the input files

Flow('line',rawsegy,'cat axis=2 ${SOURCES[1:3]}')
Flow('tline',map(lambda x: 't'+x,rawsegy),'cat axis=2 ${SOURCES[1:3]}')

# display the first 3000 traces

Plot('first','line',
     'window n2=3000 | grey pclip=90 title="First 3,000 Traces" ')

Result('zoomfirst','line',
       '''
       window min2=1000 max2=1250 min1=0 max1=3 |
       grey pclip=90 title=Zoom
       ''')

Result('firstrec24','line',
       '''
       window min2=2324 max2=2424 min1=0 max1=6 |
       grey pclip=90 title="First Record"
       ''')

# Window out misfire
Flow('shotmask','tline',
     '''
     window n1=1 f1=2 | mask min=157 max=157 |
     add add=-1 | add scale=-1
     ''')

# Convert to shots
maskshot = '''
headerwindow mask=${SOURCES[1]} |
put n2=101 n3=67 |
window n2=96 f3=10 |
put o3=100 d3=1 label3=Shot
o2=-5225 d2=110 label2=Offset unit2=ft
'''

Flow('shots',  'line shotmask',maskshot)
Flow('tshots','tline shotmask',maskshot)

Plot('shots','grey title=Shots gainpanel=all pclip=90',view=1)

Result('shots',
       '''
       byte gainpanel=all pclip=90 |
       grey3 frame1=1000 frame2=50 frame3=30 title=Shots 
       flat=n point1=0.7 point2=0.7
       ''')

###########################
# Previous Stack
###########################
Fetch('31_81_PR.SGY',
      server='http://certmapper.cr.usgs.gov',
      top='nersl/NPRA/SEISMIC/1981/31_81',
      dir='PROCESSED')

file='checkstack'
Flow([file, file+'.bin',  file+'.asc', 't'+file],
     '31_81_PR.SGY',
     '''
     segyread bfile=${TARGETS[1]} hfile=${TARGETS[2]} tfile=${TARGETS[3]}
     ''')

Result('checkstack','grey title="Check Stack" pclip=90')
Result('zoomcheckstack','checkstack',
       '''
       window max1=3 |
       grey title="Zoom Check Stack" pclip=90
       ''')

# Local slope
Flow('dip','checkstack','fdip rect1=50 rect2=10 verb=y')
Result('dip','grey color=j scalebar=y title=Slope')

Flow('pwd','checkstack dip','pwd dip=${SOURCES[1]}')
Result('pwd','grey title=PWD')
Result('zoompwd','pwd',
       '''
       window max1=3 |
       grey title="Zoom PWD" pclip=90
       ''')

#############################
# Elevation profile, datuming
#############################

Flow('spelev','../line31-81/spnElev.txt',
     '''
     echo in=$SOURCE data_format=ascii_float n1=2 n2=33 |
     dd form=native 
     ''',stdin=0)
Plot('spelev',
     '''
     dd type=complex |
     window |
     graph wanttitle=n wantaxis=n symbol=o plotcol=5
     min2=75 max2=125 min1=85 max1=172
     ''')

Flow('elev','spelev','transp | linear o1=85 d1=1 n1=88 rect=2 niter=100')
Plot('elev',
     '''
     graph title=Elevation label1=Shot label2=Elevation unit2=ft
     min2=75 max2=125 min1=85 max1=172
     ''')
Result('elev','elev spelev','Overlay')

# Shot elevation
Flow('spoint','tshots','window n1=1 f1=2 | dd type=float') # fldr
Flow('selev','elev spoint','inttest1 coord=${SOURCES[1]} interp=lag nw=2')
Result('selev','window n1=1 | graph title="Shot Elevation" label2=Elevation unit2=ft')

# observer log says shot holes 82.5 ft before spn 123, then changes to 67.5 ft.
Flow('sdepth','spoint',
     'mask max=123 | dd type=float | math output="82.5*input+67.5*(1-input)" ')

# Receiver elevation
Flow('gpoint','tshots','window n1=1 f1=2 | dd type=float | math output="input+x1/440" ')
Flow('gelev','elev gpoint','inttest1 coord=${SOURCES[1]} interp=lag nw=2')

vnear = 10000 # replacement velocity (ft/s)

# Elevation statics
Flow('estat','selev sdepth gelev','add scale=1,-1,1 ${SOURCES[1:3]} | scale dscale=%g' % (1.0/vnear))
Result('estat','grey color=j mean=y scalebar=y barlabel=Time barunit=s title="Elevation Statics" ')

Flow('eshots','shots estat','datstretch datum=${SOURCES[1]}')

Result('eshots',
       '''
       byte gainpanel=all pclip=90 |
       grey3 frame1=1000 frame2=50 frame3=30 title="Elevation Statics"
       flat=n point1=0.7 point2=0.7
       ''')

############################
# Gain, mute, and dip filter
############################

Flow('gshots','eshots',
     '''
     pow pow1=1 | put d3= |
     mutter half=n v0=%g tp=0 | put d3=1 |
     shapeagc rect1=250
     ''' % (5280/0.488))

Result('gshots',
       '''
       byte gainpanel=all |
       grey3 frame1=1000 frame2=50 frame3=30 title="Gain & Mute"
       flat=n point1=0.7 point2=0.7
       ''')

# dip filter
Flow('fft','gshots','fft1 | fft3')

Result('fft',
       '''
       window max1=100 |
       math output="abs(input)" | real |
       byte gainpanel=all allpos=y |
       grey3 frame1=200 frame2=100 frame3=30 title="2-D Fourier"
       flat=n point1=0.7 point2=0.7 color=j
       ''')

Flow('filt','fft','dipfilter v1=-167 v2=-250 v3=84 v4=72 taper=2 pass=0')

Result('filt',
       '''
       window max1=100 |
       math output="abs(input)" | real |
       byte gainpanel=all allpos=y |
       grey3 frame1=200 frame2=100 frame3=30 title="2-D Fourier"
       flat=n point1=0.7 point2=0.7 color=j
       ''')


#sfdipfilter < in.rsf > out.rsf dim=2 angle=n v=-1. ang1=-50. ang2=-45. ang3=45. ang4=50. v1=0. v2=0.1 v3=99999. v4=999999. pass=y

#       bias=.003 \
#        dx=1 \
#        slopes=-.006,-.004,.012,.014 \
#        amps=0.0,1.0,1.0,0.0 

End()