from rsf.proj import *
from m2d import *

# --------------------------------------------------------------------
# download
Fetch(['bpaitvel.hh','bpden.hh'],'bpait')

#Fetch('density_z6.25m_x12.5m.segy.gz',dir='2004_BP_Vel_Benchmark', 
#          server='ftp://software.seg.org',top='pub/datasets/2D')
#Flow('density_z6.25m_x12.5m.segy','density_z6.25m_x12.5m.segy.gz',
#          'zcat $SOURCE',stdin=0)

Flow('bpaitvel','bpaitvel.hh', '''dd form=native | put
     label1=Depth\ Z label2=Distance\ X unit1=m unit2=m''')

Flow('bpden','bpden.hh','dd form=native')

#Flow('bpden','density_z6.25m_x12.5m.segy',
#     '''
#     segyread | scale dscale=1000 | 
#     put label1=Depth\ Z label2=Distance\ X unit1=m unit2=m 
#     o1=0 o2=0 d1=6.25 d2=12.5
#     ''')

Flow('vel','bpaitvel',
     '''
     window j1=2 j2=1 | 
     window n2=1001 f2=2400 |
     put label1=Depth unit1=m label2=Distance unit2=m 
         d3=0 n3=1 o3=0 o2=0.0
     ''')

Flow('den','bpden',
     '''
     window j1=2 j2=1 | 
     window n2=1001 f2=2400  |
     put label1=Depth unit1=m label2=Distance unit2=m
         d3=0 n3=1 o3=0 o2=0.0
     ''')

# ------------------------------------------------------------------------------
model = {
    'X':    12500,
    'Z':    11937.5,
    'T':    1.5,
    'dt':   0.0015,
    'SelT':  1.4,
    'dx':   12.5,
    'dz':   12.5,
    'snpintvl': 1.0,
    'size'    : 8,   # FD order
    'frqcut'  : 0.9,
    'pml'     : 30,
    'vel'     : '4000',
    'den'     :  '2100'
}

# source & receiver
srp = {
     'bgn'   : 0.1,     # s, time of maximum ricker
     'frq'   : 20.0,     # source domain frequence
     'srcmms'  : 'y',      # point source
     'inject': 'y',      # if y, inject; if n, Initiate conditon
     'slx'   : 6250,   # source location (x), meter
     'slz'   : 5975,   # source location (x), meter
     'gdep'  : 1800      # receiver location (z), meter
     }


par = setpar(model, srp)

for m in ['den','vel']:
    pml  = m+'_pml'
    iwname = 'iw'+m
    Flow(pml,m,
         '''
         expand left=%(bd)d right=%(bd)d 
                top=%(bd)d  bottom=%(bd)d
         '''%par)
    Flow(iwname, m,'math output="input/1000" ')

Result('den',
       '''
       math output="input/1000" |
       put unit1="km" unit2="km" d1=%g d2=%g |
       grey color=j scalebar=y allpos=y transp=y screenht=7 screenratio=1
       wanttitle=n barlabel='Density' barunit='g/cm\^3\_'
       bartype=h   labelsz=6 titlesz=6 
       '''%(par['dz']/1000,par['dx']/1000))
Result('vel',
       '''
       math output="input/1000" |
       put unit1="km" unit2="km" d1=%g d2=%g |
       grey color=j scalebar=y allpos=y transp=y screenht=7 screenratio=1
       wanttitle=n labelsz=6 titlesz=6 
       bartype=h    barlabel='Velocity' barunit='km/s'
       '''%(par['dz']/1000,par['dx']/1000))


buildmodel(par, 'den0', 'vel0', 'iwden0', 'iwvel0', par['den'], par['vel'])
# ------------------------------------------------------------------------------
# build source 

Flow('srcp', None,
         '''
         spike n1=%(nt)d d1=%(dt)g k1=%(bgnp)g |
         ricker1 frequency=%(frq)g | 
         scale axis=1 
         '''%par)


# -----------------------------------------------------------------------------
# lr
lrw = 'lrwav'
lrr = 'lrrec'
lrs = 'lrsnap'
lrtr= 'lrtr'
lrtrfft = 'lrtrfft'
sglr2(lrw, lrr, 'srcp', 'vel', 'den', par, '', 0)
Flow(lrs, lrw, 'window n3=1 f3=%(snt)d ' %par)
Result(lrs,
       '''
       put unit1="km" unit2="km" d1=%g d2=%g |
       grey  screenht=7 screenratio=1
       wanttitle=y  title="snapshot by SGL method"
       labelsz=6 titlesz=6 
       bartype=h    barlabel='Velocity' barunit='m/s'
       '''%(par['dz']/1000,par['dx']/1000))

# lfd
lfdw = 'lfdwav'
lfdr = 'lfdrec'
lfds = 'lfdsnap'
lfdtr= 'lfdtr'
lfdtrfft = 'lfdtrfft'
sglfd2(lfdw, lfdr, 'srcp', 'vel_pml', 'den_pml', par, '', 0)
Flow(lfds, lfdw, 'window n3=1 f3=%(snt)d ' %par)
Result(lfds,
       '''
       put unit1="km" unit2="km" d1=%g d2=%g |
       grey  screenht=7 screenratio=1
       wanttitle=y title="snapshot by SGLFD method"
       labelsz=6 titlesz=6 
       bartype=h    barlabel='Velocity' barunit='m/s'
       '''%(par['dz']/1000,par['dx']/1000))
# iwave
iww = 'iwwav'
iwr = 'iwrec'
iws = 'iwsnap'
iwtr= 'iwtr'
iwtrfft = 'iwtrfft'
#iw(iww, iwr, 'srcp', 'iwvel', 'iwden', par, '', 0)
#Flow(iws, 'movie_p', 'window n3=1 f3=%(snt)d ' %par)
#Result(iws,
#       '''
#       put unit1="km" unit2="km" d1=%g d2=%g |
#       grey  screenht=7 screenratio=1
#       wanttitle=y title="snapshot by IWAVE" labelsz=6 titlesz=6 
#       bartype=h    barlabel='Velocity' barunit='m/s'
#       '''%(par['dz']/1000,par['dx']/1000))

End()