sconstruct

国外免费地震资料处理软件包

from rsfproj import *
from galilee import *

lag = [1,2] + range(98,103) + range(198,203)
ns = len(lag)

Flow('lag.asc',None,
     '''
     echo %s n1=%d n=100,100 in=$TARGET
     data_format=ascii_int
     ''' % (string.join(map(str,lag),' '),ns))
Flow('lag','lag.asc','dd form=native')

tension = (0,0.3,0.7,1)
nt = len(tension)

g = 'grey pclip=100 title="%s" wantaxis=n crowd=.85 gpow=.7'

Galilee('mesh',nx=560,ny=880)

Plot('mesh',
     '''
     grey allpos=y pclip=94 wantaxis=y yreverse=n transp=n 
     title="Binned" label2=South-North unit2=km label1=West-East unit1=km
     ''')
Result('mesh','Overlay',vppen='xsize=3 ysize=4')

tens = []
fils = []
cross = []
for tt in range(nt):
    t = tension[tt]
    s0 = 18*(114-73*t)
    s = [96*(8*t-11),
          84*(1-t),
          16*(9-8*t),
          112*(2*t-3)]
    s.extend([s[0],s[3],s[2],5*t-6])
    s.extend([s[2],s[1],s[2],s[7]])
    
    ss = 'ss%d' % tt
    Flow(ss+'.asc','lag',
         '''
         echo %s n1=%d lag=$SOURCE in=$TARGET
         data_format=ascii_float a0=%g
         ''' % (string.join(map(str,s),' '),ns,s0))
    Flow(ss,ss+'.asc','dd form=native')
    filt = 'flt%d' % tt
    slag = 'lag%d' % tt
    Flow([filt,ss+'0',slag,'s'+slag],ss,
         '''
         wilson eps=5.e-3 lagout=${TARGETS[2]} > ${TARGETS[1]} &&
         wilson eps=5.e-3 < ${SOURCES[0]} lagout=${TARGETS[3]} lagin=${TARGETS[2]} 
         ''')

    inp = 'inp%d' % tt
    Flow('s'+inp,filt,
         '''
         spike n1=40 n2=40 nsp=2 k1=11,16 k2=8,3 mag=1,-1 |
         helicon filt=$SOURCE
         ''',stdin=0)
    Flow(inp,'s'+inp,
         '''
         spike n1=40 n2=40 nsp=2 k1=31,28 k2=24,16 mag=1,-1 |
         add $SOURCE
         ''',stdin=0)
    div = 'div%d' % tt
    adj = 'adj%d' % tt
    Flow(div,[inp,filt],'helicon filt=${SOURCES[1]} div=y')
    Flow(adj,[div,filt],'helicon filt=${SOURCES[1]} div=y adj=y')

    Plot(inp,g % ' tension=%g' % t)
    Plot(div,g % 'input/filter' + ' clip=1.5')
    Plot(adj,g % "(input/filter)/filter'" + ' clip=1.5')

    ten = 'ten%d' % tt
    Plot(ten,[inp,div,adj],'SideBySideAniso')
    tens.append(ten)

    fill = 'fill%d' % tt
    Flow(fill,['mesh',filt],'miss filt=${SOURCES[1]} eps=0.1 exact=n niter=20')
    Plot(fill,
         '''
         igrad | grey pclip=94  yreverse=n transp=n title="tension=%g"
         ''' % t)
    fils.append(fill)

    cros = 'cros%d' % tt
    Flow(cros,fill,'window n2=1 f2=300 n1=440 f1=80')
    cross.append(cros)

Result('splin',tens[1:],'OverUnderAniso',vppen='txscale=4')

Plot('top',fils[:2],'SideBySideAniso')
Plot('bot',fils[2:],'SideBySideAniso')
Result('gal','top bot','OverUnderAniso',vppen='txscale=2 xsize=6 ysize=8')

Result('cross',cross,
       '''
       cat axis=2 ${SOURCES[1:%d]} |
       scale dscale=-1 |
       dots gaineach=0 overlap=1 dots=0 connect=4
       labels=%s label1=West-East unit1=km yreverse=y
       ''' % (nt,string.join(map(lambda x: 'tension=%g' % x,tension),':')))
    
End()