from rsf.proj import *
from rsf.recipes.beg import server

Fetch('seismic3.su','cwp',server)
title = {'seismic3':'Stacked data'}

Flow('seismic3','seismic3.su',
     '''
     suread what=d suxdr=y |
     costaper nw2=150 |
     put d2=0.0125 label2=Midpoint unit2=km
     ''')

# Cosine Fourier transform
Flow('cosft','seismic3','cosft sign2=1')

# RMS velocity
Flow('vrms.asc',None,
     '''
     echo
     0.0  1.5
     1.0  2.0
     2.5  3.16
     3.45 3.21
     4.36 3.36
     5.1  3.408
     5.45 3.6
     5.95 3.65
     n1=2 n2=8 in=$TARGET
     data_format=ascii_float
     ''')
Flow('vrms','vrms.asc seismic3',
     '''
     dd form=native |
     spline pattern=${SOURCES[1]}
     ''')

# Interval velocity
Flow('vint','vrms','dix rect1=50')

Result('vint','vrms vint',
       '''
       cat axis=2 ${SOURCES[1]} |
       graph wanttitle=n dash=1,0
       label1=Time unit1=s
       label2=Velocity unit2=km/s
       transp=y yreverse=y
       ''')

# Stolt migration
Flow('stoltm','cosft vint',
     '''
     stoltstretch velocity=${SOURCES[1]} nstretch=5 vel=1.5 pad=4000 |
     stolt vel=1.5 minstr=0.5
     ''',split=[2,2142])
Flow('stolt','stoltm vint',
     '''
     stoltstretch velocity=${SOURCES[1]} nstretch=5 vel=1.5 inv=y |
     cosft sign2=-1
     ''')
title['stolt'] = 'Stolt time section'

# Phase-shift migration
Flow('phasem','cosft vint','gazdag velocity=${SOURCES[1]} verb=y',split=[2,2142])
Flow('phase','phasem','cosft sign2=-1')
title['phase'] = 'Phase-shift time section'

# Convert time to depth

for case in ('stolt','phase'):
    Flow(case+'-depth',[case,'vint'],
         '''
         time2depth  velocity=${SOURCES[1]} intime=y nz=1500 dz=0.003 |
         put label1=Depth unit1=km
         ''')
    title[case+'-depth'] = title[case].replace('time','depth')

# Display
for case in title.keys():
    Result(case,'grey title="%s"' % title[case])

End()