from rsf.proj import *
import fdmod
## 
 # Sigsbee 2A
##
data='sigsbee2a_stratigraphy.sgy'
Fetch(data,'sigsbee')
Flow(['vstr2A','vstr2Ahead'],data,
         '''
         segyread tape=$SOURCE tfile=${TARGETS[1]} |
         put
         o1=0  d1=25     label1=z unit1=ft
         o2=10000 d2=25  label2=x unit2=ft
         ''',stdin=0)
# ------------------------------------------------------------
par = {
    'nt':7000, 'dt':0.001,'ot':0,    'lt':'t','ut':'s',
    'kt':100,    # wavelet delay
    'nx':3201, 'ox':0, 'dx':0.00762, 'lx':'x','ux':'km',
    'nz':1201, 'oz':0, 'dz':0.00762, 'lz':'z','uz':'km',
    }
# add F-D modeling parameters
fdmod.param(par)
# ------------------------------------------------------------
# wavelet
Flow('wav',None,
         '''
         spike nsp=1 mag=1 n1=%(nt)d d1=%(dt)g o1=%(ot)g k1=%(kt)d |
         ricker1 frequency=15 |
         scale axis=123 |
         put label1=t label2=x label3=y |
         transp
         ''' % par)
Result('wav',
       'transp | window n1=200 | graph title="" label1="t" label2= unit2=')
# ------------------------------------------------------------
# experiment setup
Flow('r_',None,'math n1=%(nx)d d1=%(dx)g o1=%(ox)g output=0' % par)
Flow('s_',None,'math n1=1      d1=0      o1=0      output=0' % par)
# receiver positions
Flow('zr','r_','math output="0" ')
Flow('xr','r_','math output="x1"')
Flow('rr',['xr','zr'],
     '''
     cat axis=2 space=n
     ${SOURCES[0]} ${SOURCES[1]} | transp
     ''', stdin=0)
Plot('rr',fdmod.rrplot('',par))
# source positions
Flow('zs','s_','math output=.01')
Flow('xs','s_','math output=10.0')
Flow('rs','s_','math output=1')
Flow('ss',['xs','zs','rs'],
     '''
     cat axis=2 space=n
     ${SOURCES[0]} ${SOURCES[1]} ${SOURCES[2]} | transp
     ''', stdin=0)
Plot('ss',fdmod.ssplot('',par))

# ------------------------------------------------------------

# Velocity
Flow('vel','vstr2A',
     '''
     scale rscale=.0003048 |
     put o1=%(oz)g d1=%(dz)g  o2=%(oz)g d2=%(dz)g
     ''' % par)

Plot('vel',fdmod.cgrey('''
allpos=y bias=1.5 pclip=100 color=j title=Survey\ Design labelsz=4 titlesz=6 wheretitle=t
''',par))
Result('vel',['vel','rr','ss'],'Overlay')

# ------------------------------------------------------------

# density
Flow('den','vel','math output=1')

# ------------------------------------------------------------
# finite-differences modeling
fdmod.awefd('dat','wfl','wav','vel','den','ss','rr','free=y dens=y',par)

Plot('wfl',fdmod.wgrey('pclip=99 title=Wavefield\ Movie labelsz=4 titlesz=6 wheretitle=t',par),view=1)

times=['1','2','3','4']
cntr=0
for item in ['9','19','29','39']:
    Result('time'+item,'wfl',
           '''
           window f3=%s n3=1 min1=0 min2=0 | grey gainpanel=a
           pclip=99 wantframenum=y title=Wavefield\ at\ %s\ s labelsz=4
           titlesz=6 screenratio=.375 screenht=2 wheretitle=t
           label1=z label2=x unit1=kft unit2=kft
           ''' % (item,times[cntr]))
    cntr = cntr + 1
Result('dat','window j2=4 j1=2 | transp |' + fdmod.dgrey('pclip=99',par))
End()