from rsf.proj import *

#==============================================================
# job running parameters

xv0=12.192
zv0=0

par = dict(
    ns=1376,
    dt=0.0008,

    isep=0,
    ihomo=0,
    itaper=1,
#    tapertype="DaveHale",
#    tapertype="Sinc",
#    tapertype="Cosine",
    tapertype="TTI",
#    tapertype="None",
    nstep=2,
    nxv=1400, ox=xv0, lx='x',  ux='km',
    nzv=1320, oz=zv0, lz='z',  uz='km',
    )
# =============================================================
# download Hess VTI models from 'ftp://software.seg.org'
# Depth grid: 1501 
# Horizontal grid: 3617
# Vertical and horizontal spacing are both 20ft.

from rsf.gallery import hessvti

#hessvti.get_smallmodel('vp vs epsilon delta')
hessvti.get_model('vp vs epsilon delta')

ft2km = 0.0003048
dx = 20*ft2km
dz = 20*ft2km

for mod in Split('vp vs epsilon delta'):
        Flow(mod+'1',mod,
             '''
             put d1=%g d2=%g unit1=km label1=Depth unit2=km label2=Distance unit=m/s |
             window n1=%d min1=%g n2=%d min2=%g 
             ''' % (dx, dz,
                    par['nzv'],par['oz'],
                    par['nxv'],par['ox']))

Flow('vs_hess','vp_hess',
     '''
     math output="input*0.5"
     ''')

# =============================================================
# produce vp0 & vs0 model by scaling vp model 
Flow('vp0','vp1',
     '''
     math output="input*1000" | smooth rect1=5 rect2=5
     ''')

Flow('vs0','vs1',
     '''
     math output="input*1000" | smooth rect1=5 rect2=5
     ''')

Flow('epsi','epsilon1',
     '''
     smooth rect1=3 rect2=3
     ''')

Flow('del','delta1',
     '''
     smooth rect1=3 rect2=3
     ''')

# =============================================================

# produce vp0 & vs0 model by scaling vp model 
Flow('hessvp0','vp1',
     '''
     math output="input*1000" 
     ''')

Flow('hessvs0','vs1',
     '''
     math output="input*1000" 
     ''')

Flow('hessepsilon','epsilon1',
     '''
     math output="input" 
     ''')

Flow('hessdelta','delta1',
     '''
     math output="input" 
     ''')

name0='''
hessvp0 hessvs0
'''

name00='''
hessepsilon hessdelta
'''

Result('vp0',
        '''
        grey color=g scalebar=y screenht=4 screenwd=10 barreverse=y wanttitle=n}
        ''')

for ff in Split(name0):
        Result(ff,
        '''
        grey color=g scalebar=y screenht=4 screenwd=10 barreverse=y wanttitle=n}
        ''')

for gg in Split(name00):
        Result(gg,
        '''
        grey color=g scalebar=y screenht=4 screenwd=10 barreverse=y wanttitle=n}
        ''')
#       grey color=j scalebar=y screenratio=0.47 allpos=n barreverse=y wanttitle=n}
################################################################
# Stage I:Forward modeling with real SV0 
################################################################
name1='''
PseudoPurePx PseudoPurePz PseudoPureP
'''

Flow(['PseudoPurePx',  'PseudoPurePz',  'PseudoPureP'],
         'vp0  vs0  epsi del',
         '''
         vti2dpseudop
         vp0=${SOURCES[0]}
         vs0=${SOURCES[1]}
         epsi=${SOURCES[2]} 
         del=${SOURCES[3]}
         PseudoPurePz=${TARGETS[1]}
         PseudoPureP=${TARGETS[2]}
         ns=%d 
         dt=%g
         isep=%d
         ihomo=%d
         itaper=%d
         tapertype=%s
         nstep=%d
         ''' % (par['ns'],par['dt'],par['isep'],par['ihomo'],par['itaper'],par['tapertype'],par['nstep'])
    )

for qq in Split(name1):
        Result(qq,
        '''
        grey polarity=n scalebar=n screenratio=1 wanttitle=n
        ''')

################################################################
# Stage II:Forward modeling with finite SV0 
################################################################
name2='''
PseudoPurePxFsv0 PseudoPurePzFsv0 PseudoPurePFsv0
'''

Flow(['PseudoPurePxFsv0',  'PseudoPurePzFsv0',  'PseudoPurePFsv0', 'Fvs0'],
         'vp0  epsi del',
         '''
         vti2dpseudopfvs0
         vp0=${SOURCES[0]}
         epsi=${SOURCES[1]} 
         del=${SOURCES[2]}
         PseudoPurePz=${TARGETS[1]}
         PseudoPureP=${TARGETS[2]}
         Fvs0=${TARGETS[3]}
         ns=%d 
         dt=%g
         isep=%d
         ihomo=%d
         itaper=%d
         tapertype=%s
         nstep=%d
         ''' % (par['ns'],par['dt'],par['isep'],par['ihomo'],par['itaper'],par['tapertype'],par['nstep'])
    )

for qq in Split(name2):
        Result(qq,
        '''
        grey polarity=n scalebar=n screenratio=1 wanttitle=n
        ''')
################################################################
# Stage III: Coverted the RTM (cjb-SU) data to RSF format and plot
################################################################
dx=6.096*2/1000
dz=6.096*2/1000
nx0=3617

# Not really reproducible, just display

from rsf.recipes.beg import server

Fetch('hessrtmillum.su','hess',server)

Flow('hessrtm','hessrtmillum.su',
     '''
     datasucjb2rsf2d fn=$SOURCE
     | window n2=1809 min2=0
     | put d1=0.012192 d2=0.012192 unit1=km label1=Depth unit2=km label2=Distance
     ''')

Result('hessrtm',
     '''
     grey pclip=96 color=e polarity=n scalebar=n screenht=6 screenwd=10 wanttitle=n
     ''')

End()