from rsf.proj import *

choose = {
'dn': '1',
'vp': '2'
}

label = {
'dn': 'Density',
'vp': 'Velocity'
}

unit = {
'dn': 'g/cm\^3\_',
'vp': 'km/s'
}

for m in ('dn','vp'):
    # Get model
    mod = m + '0.rsf'
    # Fetch(mod,'iwave')

    Flow([mod, mod +'@'], None,
         '''
         standardmodel 
         model=8 
         hfile=
         '''
         + mod +
         '''
         choose=
         '''
         + choose[m] +
         '''
         d1=2.5 d2=2.5 d3=1.0 
         n1=721 n2=3121 n3=1 
         f3=3300
         label1=Depth unit1=m
         label2=Distance unit2=m
         label=
         '''
         + label[m] +
         '''
         unit=
         '''
         + unit[m] +
         '''
         datapath=.
         ''',stdin=0,stdout=0)

    j = 2 # subsampling
    for r in range(1,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','g/cm\^3\_')[m=='dn']))

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

        j *= 2

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')

# 2-4 results
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,'hdr.su','demo.par'],
         '''
         /bin/cat ${SOURCES[3]} > $TARGET &&
         echo
         hdrfile  = ${SOURCES[2]}
         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
         ''',stdin=0)

    Result(data,'grey title="%g-m data, 2-4 scheme" ' % (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=4,2,6,7 wanttitle=n label2=Pressure unit2=MPa
       ''')

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

# 2-8 results
j = 1
traces8k = []
for r in range(4):

    dn = 'dn%d' % r
    vp = 'vp%d' % r

    par = 'iwave8k%d.par' % r
    data = 'data8k%d' % r
    sudata = data + '.su'

    Flow(par,[dn,vp,'hdr.su','demo.par'],
         '''
         /bin/cat ${SOURCES[3]} | sed s/order=2/order=4/ > $TARGET &&
         echo
         hdrfile  = ${SOURCES[2]}
         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
         ''',stdin=0)

    trace = 'trace' + str(r+4)
    Flow(trace,data,'window n2=1 f2=100')
    traces8k.append(trace)

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

    j *= 2

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



End()