sconstruct

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

from rsfproj import *

# Make a hyperbolic reflector
#############################

a = 1
b = 0
# b = 1 # !!! Uncomment this


Flow('hyper',None,
     'math n1=5001 o1=-2.5 d1=0.001 output="sqrt(%g+%g*x1*x1)" ' % (a*a,b*b))

# Reflector dip

Flow('dip','hyper','math output="%g*x1/input" ' % (b*b))


Result('hyper',
       '''
       unif2 n1=201 d1=0.01 v00=1,2 allpos=y |
       grey color=j screenratio=0.4 screenht=5
       title=Model
       label1=Depth unit1=km
       label2=Lateral unit2=km
       ''')

# Kirchhoff modeling
####################

ns = 3  # number of shot gathers
s0 = -1 # first shot 
ds = 1  # shot increment

vel = 1 # velocity

Flow('model','hyper dip',
     '''
     kirmod ns=%d ds=%g s0=%g nh=201 dh=0.01 h0=-1
     vel=%g freq=10 nt=875 dt=0.004 dip=${SOURCES[1]}
     ''' % (ns,ds,s0,vel))

shots = []
for s in range(3):
    # Plot shot gather
    shot = 'shot%d' % s
    Flow(shot,'model','window n3=1 f3=%d' % s)
    Plot(shot+'g',shot,
         '''
         grey title="Shot at %g"
         label1=Time unit1=s
         label2=Offset unit2=km
         labelsz=10 titlesz=15
         ''' % (s0+s*ds))

    # !!! Plot analytical formula - change for nonzero b
    time = 'time%d' % s
    Plot(time,'model',
         '''
         window n3=1 n1=1 | math output="sqrt(x1*x1+%g)/%g" |
         graph yreverse=y min2=0 max2=3.5 pad=n
         wanttitle=n plotcol=2 plotfat=5 wantaxis=n
          ''' % ((4*a*a),vel))

    Plot(shot,[shot+'g',time],'Overlay')    
    shots.append(shot)

Result('shots',shots,'SideBySideAniso')

# Wave-equation migration
#########################

# Spike wavelet

fft = 'fft1 | window max1=30'

Flow('wav',None,'spike n1=875 k1=1 | ' + fft)

# Slowness on a grid

Flow('slow',None,
     '''
     spike mag=%g
     n3=201 d3=0.01
     n2=1 d2=0.01
     n1=500 d1=0.01
     ''' % (1./vel))

# Shot and receiver wavefields

Flow('rwfl swfl','model wav',
     fft + ''' |
     spray axis=3 n=1 o=0 d=1 |
     spray axis=5 n=1 o=0 d=1 |
     srsyn nx=500 dx=0.01 ox=-2.5 wav=${SOURCES[1]} swf=swf.rsf |
     transp plane=12 | transp plane=23 > ${TARGETS[0]} &&
     transp plane=12 < swf.rsf | transp plane=23 > ${TARGETS[1]} &&
     rm swf.rsf
     ''',stdout=-1)

# Migrate and create common-image gathers

Flow('image cigs','swfl rwfl slow',
     '''
     srmig2 slo=${SOURCES[2]} rwf=${SOURCES[1]} cig=${TARGETS[1]}
     itype=x nhx=51 nhz=1 hsym=y jcx=50 readwrite=y verb=y 
     ''')

# Display an image

Result('image',
       '''
       window | grey pclip=100 screenratio=0.4 screenht=5
       title=Image transp=n
       label2=Depth unit2=km
       label1=Lateral unit1=km
       ''')

# Display common-image gathers

cigs = []
for s in range(3):
    cig = 'cig%d' % s
    x = -2.5+(s+4)*0.5
    Flow(cig,'cigs',
         '''
         window n1=1 f1=%d | slant adj=y p0=-1.2 np=61 dp=0.04 |
         window min1=0.5 max1=1.5
         ''' % (s+4))
    Plot(cig,
         '''
         grey pclip=100
         label1=Depth unit1=km
         label2="Tan(theta)"
         title="Angle Gather at x=%g" grid=y gridcol=5
         labelsz=10 titlesz=15
         ''' % x)
    cigs.append(cig)

Result('cigs',cigs,'SideBySideAniso')

End()