from rsf.proj import *


# Download Hess VTI model
#########################
zcat = WhereIs('gzcat') or WhereIs('zcat')
for case in ('vp','epsilon'):
    sgy = 'timodel_%s.segy' % case
    sgyz = sgy + '.gz'
    Fetch(sgyz,dir='Hess_VTI',
          server='ftp://software.seg.org',
          top='pub/datasets/2D')
    # Uncompress
    Flow(sgy,sgyz,zcat + ' $SOURCE',stdin=0)
    # Convert to RSF format
    Flow(case,sgy,
         '''
         segyread read=data | 
         window j1=2 j2=2 | put d1=40 d2=40 
         unit1=ft label1=Depth unit2=ft label2=Distance 
         ''')

# Horizontal velocity
Flow('vx','vp epsilon',
     'math e=${SOURCES[1]} output="input*sqrt(1+2*e)"')

for case in ('vp','vx'):
    Result(case,
           '''
           grey color=j pclip=100 allpos=y bias=5000 
           scalebar=y barreverse=y wanttitle=n 
           barlabel=Velocity barunit=ft/s 
           screenht=5 screenwd=12 labelsz=6
           ''')

Flow('wav','vp vx','sfTestaniso vx=${SOURCES[1]} nb=30 fm=10 nt=1000 dt=0.001 verb=y')
Result('snapaniso','wav',
        '''
        window n3=1 min3=0.9 |grey title="Hess VTI:t=0.9s" label1=Depth unit1=km 
        label2=Lateral unit2=km screenht=5 screenwd=10 labelsz=6
        ''')



End()