from rsf.proj import *

# Download velocity model from the data server
##############################################
vstr = 'sigsbee2a_stratigraphy.sgy'
Fetch(vstr,'sigsbee')
Flow('zvstr',vstr,'segyread read=data')

Flow('vel','zvstr',
     '''
     put d1=0.00762 o2=3.048 d2=0.00762
     label1=Depth unit1=km label2=Distance unit2=km
     label=Velocity unit=km/s |
     scale dscale=0.0003048
     ''')
Result('vel',
       '''
       grey wanttitle=n color=j allpos=y 
       screenratio=0.3125 screenht=4 labelsz=4
       scalebar=y barreverse=y
       ''')

# Window a portion
Flow('vel2','vel','window max2=9.5')

dt = 0.002
nt = 5001

# Convert to RMS
Flow('voft','vel2',
     'depth2time velocity=$SOURCE dt=%g nt=%d' % (dt,nt))
Flow('vrms','voft',
     '''
     add mode=p $SOURCE | causint | 
     math output="sqrt(input*%g/(x1+%g))" |
     smooth rect2=10 | window j2=2
     ''' % (dt,dt))

# Download zero-offset from the data server
###########################################
Fetch('sigexp_ns.rsf','sigsbee')
Flow('data','sigexp_ns.rsf',
     '''
     dd form=native |
     bandpass flo=2 fhi=60 |
     window max2=9.5 j1=2 j2=2 n1=%d |
     costaper nw1=50 nw2=50
     ''' % nt)

Result('data',
       'window min1=3 | grey title="Zero-Offset Data" ')

# Slope estimation
Flow('dip','data','fdip rect1=100 rect2=10')
Result('dip',
       '''
       grey color=j scalebar=y
       title="Dominant Slope"
       barlabel=Slope barunit=samples
       ''')

# Plane-wave destruction
Flow('dif','data dip','pwd dip=${SOURCES[1]}')
Result('dif',
       '''
       window min1=3 | 
       grey title="Separated Diffractions" 
       ''')

# Velocity continuation
Flow('fourier','dif','pad n2=1025 | cosft sign2=1')
Flow('velconf','fourier',
     '''
     put o3=0 | stolt vel=1.4 | 
     spray axis=2 n=1 o=0 d=1 | 
     fourvc pad2=8192 nv=61 dv=0.02 v0=1.4 verb=y
     ''',split=[2,1025],reduce='cat axis=3')
Flow('velcon','velconf',
     '''
     transp plane=23 memsize=1000 | 
     cosft sign2=-1  | window n2=424 | 
     transp plane=23 memsize=1000
     ''')

# Picking a slice
#################
Flow('dimage','velcon vrms',
     'slice pick=${SOURCES[1]}')
Result('dimage',
       '''
       window min1=3 |
       grey title="Imaged Diffractions"
       ''')

# Angle-gather migration
########################
proj = Project()
prog = proj.Program('anglemig.c')

Flow('anglemig','dif %s' % prog[0],
     './${SOURCES[1]} vel=2 na=90 a0=-45 da=1')

Result('anglemig',
       '''
       window min2=2 |
       transp | transp plane=23 memsize=1000 | 
       byte gainpanel=all | grey3
       frame1=1000 frame2=200 frame3=60 unit3="\^o\_"
       title="Dip Angle Gathers" point1=0.7 point2=0.7
       ''')

End()