sconstruct

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

from rsfproj import *
import math, string

def sag(nt,nx,alpha):
    k1 = []
    k2 = []
    for ix in range(8):
        z = 0.25 + ix*0.5
        if alpha > 0:
            t = 2 * math.log(1 + alpha*z/1.5) / alpha
        else:
            t = 2 * z/1.5
        it = int(1.5+t*nt/4)
        ix = int(1.5+z*nx/4)
        if it <= nt:
            k1.append(str(it))
            k2.append(str(ix))
    return '''
    spike n1=%d n2=%d d1=%g d2=%g nsp=%d k1=%s k2=%s
    ''' % (nt,nx,4.0/nt,4.0/nx,len(k1),string.join(k1,','),string.join(k2,','))

Flow('sagmod',None,sag(64,64,0.5))
Result('sagmod','grey pclip=100 wanttitle=n')

Flow('sagvel','sagmod','voft v0=1.5 alpha=0.5')
Flow('sagdat','sagmod sagvel',
     'cosft sign2=1 | gazdag velocity=${SOURCES[1]} inv=1 | cosft sign2=-1')
Result('sagdat','grey pclip=100 wanttitle=n')

for case in range(2):
    alpha = case*0.5
    mod = 'mod%d' % case
    vel = 'vel%d' % case
    Flow(mod,None,sag(64,64,alpha) + ' | smooth rect1=3')
    Flow(vel,mod,'voft v0=1.5 alpha=%g' % alpha)
    Plot(mod,[mod,vel],
         '''
         cosft sign2=1 | gazdag velocity=${SOURCES[1]} inv=1 | 
         tpow tpow=0.5 | gazdag velocity=${SOURCES[1]} inv=0 |
         cosft sign2=-1 | grey pclip=100 wanttitle=n
         ''')

Result('sagres','mod0 mod1','SideBySideAniso',vppen='txscale=1.5')

for file in ('beistack','beivrms'):
    Fetch(file+'.HH','midpts')
    Flow(file,file+'.HH','dd form=native')

Flow('wgphase',['beistack','beivrms'],
     'cosft sign2=1 | gazdag velocity=${SOURCES[1]} | cosft sign2=-1')
Result('wgphase','wgphase','grey pclip=98 wanttitle=n')

End()