from rsf.proj import *

for m in ('dn','vp'):
    # Get model
    mod = m + '2d_2.5m.hh'
    Fetch(mod,'iwave')

    j = 1 # subsampling
    for r in range(4):
        modrsf = '%s%d' % (m,r)
        Flow(modrsf,mod,
             '''
             dd form=native | window j1=%d j2=%d |
             put label1=Depth unit1=m label2=Distance unit2=m
             label=%s unit="%s"
             ''' % (j,j,('Velocity','Density')[m=='dn'],('km/s','kg/m\^3\_')[m=='dn']))

        Result(modrsf,'grey color=j mean=y title="%g-m model" scalebar=y barreverse=y' % (2.5*j))

        j *= 2

Flow('wavelet',None,
     '''
     spike n1=401 o1=-0.4 d1=0.002 k1=201 |
     trapez frequency=0.,2.5,15.0,20.0 |
     sfscale dscale=1.e6 
     ''')
Flow('thead','wavelet','segyheader')
Flow('wavelet.su','wavelet thead','suwrite tfile=${SOURCES[1]} endian=0')

Flow('data',None,'spike n1=1501 n2=301 d1=0.002 d2=1 o2=0')
Flow('slice','data','window n1=1')

hkeys = dict(sx=3300,gx='100+20*x1',delrt=0,selev=-40,gelev=-20)
hks = hkeys.keys()

i=0
hstr = ''
for hk in hks:
    Flow(hk,'slice','math output=%s | dd type=int' % str(hkeys[hk]))
    i += 1
    hstr += '%s=${SOURCES[%d]} ' % (hk,i)

Flow('tdata',['data']+hks,'segyheader ' + hstr)
Flow('hdr.su','data tdata','suwrite tfile=${SOURCES[1]} endian=0')

j = 1
traces = []
for r in range(4):
    dn = 'dn%d' % r
    vp = 'vp%d' % r

    par = 'iwave%d.par' % r
    data = 'data%d' % r
    sudata = data + '.su'

    Flow(par,[dn,vp,'wavelet.su','hdr.su','demo.par'],
         '''
         /bin/cat ${SOURCES[4]} > $TARGET &&
         echo
         source = ${SOURCES[2]}
         hdrfile  = ${SOURCES[3]}
         datafile  = %s
         velocity = ${SOURCES[1]}
         density = ${SOURCES[0]} >> $TARGET
         ''' % sudata,stdin=0,stdout=-1)

    Flow([data,sudata],par,
         '''
         asg par=$SOURCE &&
         suread < ${TARGETS[1]} read=data endian=0
         ''' % data,stdin=0)

    Result(data,'grey title="%g-m data" ' % (2.5*j))

    j *= 2

    trace = 'trace%d' % r
    Flow(trace,data,'window n2=1 f2=100')
    traces.append(trace)

Result('trace',traces,
       'cat axis=2 ${SOURCES[1:4]} | graph plotcol=7,6,2,4 wanttitle=n label2=Pressure unit2=MPa')

Result('wtrace',traces,
       '''
       cat axis=2 ${SOURCES[1:4]} | window min1=1.8 max1=2.5 |
       graph plotcol=7,6,2,4 wanttitle=n label2=Pressure unit2=MPa
       ''')

End()