sconstruct

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

from rsfproj import *

yshot=0.2
zshot=0.2

contour = '''
contour c0=0. dc=0.08 nc=15 max1=1. max2=1. wanttitle=n
label1=Depth label2=Lateral screenratio=1. labelsz=12 screenht=14
'''

samples = (1,2,4,8,10,15,20,30,40,50,80,100)

Flow('const',None,
     'spike n1=2000 n2=2000 d1=0.0005 d2=0.0005 unit1=km unit2=km mag=1')
     
Flow('exact','const','math output="sqrt((x1-%g)^2+(x2-%g)^2)" ' % (zshot,yshot))
Plot('exact','window j1=5 j2=5 | ' + contour)

for order in (1,2):
    eiko = 'eiko%d' % order
    const = 'const%d' % order
    Plot(eiko,'const',
         '''
         window j1=100 j2=100 |
         eikonal order=%d zshot=%g yshot=%g br1=0.05 br2=0.05 |
         %s dash=3 wantaxis=n
         ''' % (order,zshot,yshot,contour))
    Plot(const,['exact',eiko],'Overlay')

Result('circles','const1 const2','SideBySideIso')

err = [[],[]]
tim = [[],[]]
for sample in samples:
    vel = 'vel%d' % sample
    ext = 'ext%d' % sample
    spc = 'spc%d' % sample
    siz = 'siz%d' % sample
    Flow(vel,'const','window j1=%d j2=%d' % (sample,sample))
    Flow(ext,'exact','window j1=%d j2=%d' % (sample,sample))
    Flow(spc,vel,
         'window n2=1 | math output=1 | stack axis=1 norm=n | math output=1000/input')
    Flow(siz,vel,
         'math output=1 | stack axis=2 norm=n | stack axis=1 norm=n | scale dscale=0.000001')
         
    for order in (1,2):
        eik = 'eik%d-%d' % (order,sample)
        dif = 'dif%d-%d' % (order,sample)
        Flow(eik,[vel,ext],
             '''
             eikonal order=%d zshot=%g yshot=%g br1=0.1 br2=0.1 |
             math exact=${SOURCES[1]} output="(input-exact)^2" |
             stack axis=2 | stack axis=1 | math output="1000*sqrt(input)" 
             ''' % (order,zshot,yshot))
        Flow(dif,[spc,eik],'cmplx ${SOURCES[:2]}',stdin=0)
        err[order-1].append(dif)

        time = 'tim%d-%d' % (order,sample)
        Flow(time,[vel,siz],
             '''
             (/usr/bin/time -f "%%S %%U" %s < ${SOURCES[0]} order=%d zshot=%g yshot=%g br1=0.1 br2=0.1 > /dev/null ) >& time.out &&
             (tail -1 time.out; echo in=time0.rsf n1=2 data_format=ascii_float) > time0.rsf &&
             dd form=native < time0.rsf | stack axis=1 norm=n > time1.rsf &&
             cmplx ${SOURCES[1]} time1.rsf > $TARGET &&
             /bin/rm time.out time0.rsf &&
             rm time1.rsf
             ''' % (WhereIs('sfeikonal'),order,zshot,yshot),
             stdin=0,stdout=-1)
        tim[order-1].append(time)
        
for order in (1,2):
    Flow('err%d' % order,err[order-1],'cat axis=1 ${SOURCES[1:%d]}' % len(samples))
    Flow('tim%d' % order,tim[order-1],'cat axis=1 ${SOURCES[1:%d]}' % len(samples))
    
Flow('err','err1 err2','cat axis=2 ${SOURCES[1]}')
Flow('tim','tim1 tim2','cat axis=2 ${SOURCES[1]}')

Flow('logerr','err',
     '''
     real < $SOURCE | math output="log(input)" > real.rsf &&
     imag < $SOURCE | math output="log(input)" > imag.rsf &&
     cmplx real.rsf imag.rsf > $TARGET &&
     rm real.rsf imag.rsf
     ''',stdin=0,stdout=-1)

Flow('normtim','tim',
     '''
     real < $SOURCE > size.rsf &&
     imag < $SOURCE | add mode=d size.rsf > imag.rsf &&
     cmplx size.rsf imag.rsf > $TARGET &&
     rm size.rsf imag.rsf
     ''',stdin=0,stdout=-1)

Plot('err',
     '''
     graph wanttitle=n
     label2="rms error" unit2=ms label1="grid spacing" unit1=m
     dash=0,3 screenratio=1. labelsz=12 screenht=14
     ''')
Plot('logerr',
     '''
     graph wanttitle=n label2="log(rms error)" label1="log(grid spacing)"
     unit2=" " unit1=" " dash=0,3 screenratio=1. labelsz=12 screenht=14
     ''')
Result('error','err logerr','SideBySideIso')

Plot('tim',
     '''
     graph wanttitle=n dash=0,3 screenratio=1. labelsz=9
     label2="CPU time" unit2=s label1="N/1000000" unit1=
     ''')
Plot('normtim',
     '''
     graph wanttitle=n dash=0,3 screenratio=1. labelsz=9
     label2="1000000 CPU time/N" unit2=s label1="N/1000000" unit1=
     min2=0 max2=2
     ''')
Result('times','tim normtim','SideBySideIso')

End()