from rsf.proj import *
from math import *
xmax = 6.0
zmax = 3.0
layers = ((0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00),
(0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3),
(0.5,0.5,0.5,0.5,0.65,0.9,0.65,0.5),
(0.85,0.85,0.85,0.95,1.2,1.3,1.3,1.3),
(1.55,1.55,1.55,1.55,1.55,1.55,1.55,1.55),
(2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0))
flatlayers = ((0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00),
(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0),
(2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0),
(3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0),
(4.0,4.0,4.0,4.0,4.0,4.0,4.0,4.0),
(5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0))
def arr2str(array,sep=' '):
return sep.join(map(str,array))
n1 = len(layers[0])
n2 = len(layers)
Flow('layers.asc',None,
'''
echo %s
n1=%d n2=%d o1=0 d1=%g
data_format=ascii_float in=$TARGET
''' % (' '.join(map(arr2str,layers)),
n1,n2,xmax/(n1-1)))
Flow('layers','layers.asc','dd form=native')
d = 0.0101
Flow('refl1','layers',''' window n2=1 f2=0| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('refl2','layers',''' window n2=1 f2=1| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('refl3','layers',''' window n2=1 f2=2| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('refl4','layers',''' window n2=1 f2=3| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('refl5','layers',''' window n2=1 f2=4| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('refl6','layers',''' window n2=1 f2=5| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('ref','refl1 refl2 refl3 refl4 refl5 refl6', '''cat axis=2 ${SOURCES[1:6]} | math output="input*3/2"''')
Plot('ref','''graph yreverse=y wanttitle=n label1=Distance unit1=km plotfat=6
label2=Depth unit2=km min2=0 max2=3 min1=0 max1=6
screenht=5.0 screenratio=0.333 yll=3.5 xll=1.5
axisfat=3 titlefat=3 titlesz=10 labelfat=3 labelsz=6''')
Flow('flat.asc',None,
'''
echo %s
n1=%d n2=%d o1=0 d1=%g
data_format=ascii_float in=$TARGET
''' % (' '.join(map(arr2str,flatlayers)),
n1,n2,xmax/(n1-1)))
Flow('flats','flat.asc','dd form=native')
Flow('reflf1','flats',''' window n2=1 f2=0| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('reflf2','flats',''' window n2=1 f2=1| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('reflf3','flats',''' window n2=1 f2=2| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('reflf4','flats',''' window n2=1 f2=3| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('reflf5','flats',''' window n2=1 f2=4| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('reflf6','flats',''' window n2=1 f2=5| spline o1=0 d1=%g n1=%d '''% (d,int(1.5+xmax/d)))
Flow('reff','reflf1 reflf2 reflf3 reflf4 reflf5 reflf6', '''cat axis=2 ${SOURCES[1:6]} ''')
v1 = 2.5
v2 = 3.2
v3 = 3.4
v4 = 3.7
v5 = 3.86
updown=[1,2,3,4,5,4,3,2,1,0]
N = len(updown)-1
vstat = 1
scale=0.3
gx = [-0.06*scale,0.15*scale,0.2*scale,0.3*scale,0.35*scale]
gz = [0.0,0.0,0.0,0.0,0.0]
v = [v1,v2,v3,v4,v5]
vcol = [1.508,1.581,1.69,1.826,2,2.2]
xref = [3.0,3.0,3.0,3.0,3.0]
zref = [0.0,0.375,0.6,0.7,1.45]
xrefstr = ','.join(map(str,xref))
zrefstr = ','.join(map(str,zref))
gxstr = ','.join(map(str,gx))
gzstr = ','.join(map(str,gz))
updownstr = ','.join(map(str,updown))
vstr = ','.join(map(str,v))
vcolstr = ','.join(map(str,vcol))
Flow('refcut','ref','window n2=6 | put o2=0')
Flow('xrefl','refcut',
'''isaac2 niter=500 number=%d vstatus=%d debug=n xs=2 xr=4 velocity=%s layer=%s
xref=%s zref=%s xgradient=%s zgradient=%s tol=0.0001''' % (N,vstat,vstr,updownstr,xrefstr,zrefstr,gxstr,gzstr))
Plot('xrefl',
'''
dd type=complex | window |
graph wanttitle=n wantaxis=n yreverse=y min2=0 max2=3 min1=0 max1=6 plotcol=7 plotfat=6 screenht=5.0 screenratio=0.333 yll=3.5 xll=1.5
axisfat=3 titlefat=3 titlesz=10 labelfat=3 labelsz=6
''')
Plot('xrefl-points','xrefl',
'''
dd type=complex | window |
graph wanttitle=n wantaxis=n yreverse=y min2=0 max2=3 min1=0 max1=6 plotcol=7 plotfat=6 symbol=* symbolsz=10 screenht=5.0 screenratio=0.333 yll=3.5 xll=1.5
axisfat=3 titlefat=3 titlesz=10 labelfat=3 labelsz=6
''')
Result('v','xrefl xrefl-points ref','Overlay')
Flow('mod1','ref',
'''
window n2=5 f2=1 |
unif2 d1=%g n1=%d v00=%s dvdx=%s x0=%s
''' % (d,int(1.5+zmax/d),vstr,gxstr,xrefstr))
Flow('modcol','ref',
'''
window n2=5 f2=1 |
unif2 d1=%g n1=%d v00=%s
''' % (d,int(1.5+zmax/d),vcolstr))
Flow('velocity','reff',
'''
window n2=5 f2=1 |
unif2 d1=%g n1=%d v00=%s dvdx=%s x0=%s |
transp
''' % (1,5,vstr,gxstr,xrefstr))
Result('mod1',
'''
grey color=j
screenratio=%g screenht=4 wanttitle=n
labelfat=3 labelsz=6 scalebar=y
maxval=5.2 minval=2.2 allpos=y bias=2.2 clip=3
label1="Depth (km)"
label2="Distance (km)"
barlabel="Velocity (km/s)"
''' % (zmax/xmax))
Flow('lineshot','modcol','window n2=1 f2=426 | math output=x1')
Flow('lineshotc','lineshot','math output="426*0.0101" | cat axis=2 $SOURCE | transp | dd type=complex | window ')
Plot('lineshotc','graph screenratio=%g yreverse=y screenht=4 min1=0 max1=6 wanttitle=n wantaxis=n plotfat=4 plotcol=7' % (zmax/xmax))
Plot('modcol',
'''
grey color=j
screenratio=%g screenht=4 wanttitle=n
mean=y labelfat=2 labelsz=4
label1="Depth (km)"
label2="Distance (km)"
''' % (zmax/xmax))
Result('modcol','modcol lineshotc','Overlay')
Flow('refs','ref','window n2=5 f2=1 | put o2=0')
Flow('dips','refs','deriv scale=y')
Flow('diffrac1','refs','math output="1.0" | cut n2=4 | math output="1.0 - input" ')
Flow('diffrac2','refs','math output="1.0" | cut f1=420 n1=11 | math output="1.0 - input" ')
Flow('diffractivity','diffrac1 diffrac2','add ${SOURCES[1]} | mask min=1 | dd type=float ')
shotscmp = (2.5,3.0,3.5)
plotscmp = []
Flow('zo','refs dips diffractivity',
'''
kirmod_newton nt=751 dt=0.004 freq=15 cmp=y
ns=201 s0=3.2825 ds=0.0101 nh=1 dh=0.04 h0=0 verb=y
debug=n fwdxini=y
vstatus=%d velocity=%s debug=n
xref=%s zref=%s xgradient=%s zgradient=%s dip=${SOURCES[1]} refl=${SOURCES[2]}
''' %(vstat,vstr,xrefstr,zrefstr,gxstr,gzstr))
Flow('zotaper','zo','window | pow pow1=1 | costaper nw1=25 | costaper nw2=25')
Result('zotaper',
'''
window | grey transp=y yreverse=y poly=y
title="Zero-offset section" label2=Midpoint unit2=km screenratio=1.3
axisfat=3 titlefat=3 titlesz=18 labelfat=3 labelsz=14 max1=2
wherexlabel=top
''')
Flow('depth','ref refs','window n2=5 | math s=${SOURCES[1]} output="s-input" | put o2=0')
Flow('t0','depth velocity','math v=${SOURCES[1]} output="input/v"')
Flow('v2t0','t0 velocity','math v=${SOURCES[1]} output="input*v^2" ')
Flow('t0sum','t0','transp | causint | transp')
Flow('t0sumext','t0sum','spray axis=3 n=201 d=0.0101 o=-1.01')
Flow('vnmosq','v2t0 t0sum','transp | causint | transp | math t0=${SOURCES[1]} output="input/t0" ')
Flow('vnmosqext','vnmosq','spray axis=3 n=201 d=0.0101 o=-1.01')
Flow('hypertime','vnmosqext t0sumext','math t0=${SOURCES[1]} output="2*sqrt(t0^2 + (x3)^2/input)"')
Flow('offset',None,'spike n1=201 d1=0.0101 o1=-1.01 | math output=x1')
Flow('slow','velocity',' math output="1/input" ')
Flow('vnmosqbylayer','velocity',' math output="input^2" ')
Flow('vnmosqhet','refs vnmosqbylayer slow t0sum',
'''
sffermatrecursion vnmosq=${SOURCES[1]} slow=${SOURCES[2]} t0sum=${SOURCES[3]}
''')
Flow('vnmosqhetext','vnmosqhet','spray axis=3 n=201 d=0.0101 o=-1.01')
Flow('hyperhettime','vnmosqhetext t0sumext','math t0=${SOURCES[1]} output="2*sqrt(t0^2 + (x3)^2/input)"')
s3=426
for j in [s3]:
sh = 3
shot = str(sh)
for i in range(4,5):
num = str(i+1)
Flow('zo'+shot,'zotaper',' window ')
Plot('zo'+shot,
'''
grey transp=y yreverse=y poly=y pclip=95
title="Zero-offset diffraction at %s km" min1=0 max1=2 label2=Midpoint unit2=km screenratio=1.5
axisfat=3 titlefat=3 titlesz=6 labelfat=3 labelsz=6
wherexlabel=top
''' % (str(j*0.0101)) )
Plot('hyperto'+num+'-'+shot,'hypertime offset',
'''
window n1=1 f1=%d n2=1 f2=%d | cat axis=2 ${SOURCES[1]} order=2,1 | transp | dd type=complex | window |
graph yreverse=y plotcol=%d plotfat=4 min2=0 max2=2 screenratio=1.5
axisfat=3 wanttitle=n wantaxis=n
''' % (j,i,i+1) )
Plot('hyperhetto'+num+'-'+shot,'hyperhettime offset',
'''
window n1=1 f1=%d n2=1 f2=%d | cat axis=2 ${SOURCES[1]} order=2,1 | transp | dd type=complex | window |
graph yreverse=y plotcol=%d plotfat=4 min2=0 max2=2 screenratio=1.5 dash=1
axisfat=3 wanttitle=n wantaxis=n
''' % (j,i,i+1) )
Result('hypercomparediff3', ''' zo3 hyperto5-3
hyperhetto5-3
''','Overlay')
for j in [s3]:
sh = 3
shot = str(sh)
Flow('hypertimeshift'+shot,'vnmosqext t0sumext',
''' math t0=${SOURCES[1]} output="2*sqrt(t0^2 + (x3)^2/input) - 2*t0"|
window n1=1 f1=%d n2=1 f2=4
''' % j)
Flow('hyperhettimeshift'+shot,'vnmosqhetext t0sumext',
''' math t0=${SOURCES[1]} output="2*sqrt(t0^2 + (x3)^2/input) - 2*t0"|
window n1=1 f1=%d n2=1 f2=4
''' % j)
Flow('warped'+shot,'zo'+shot+' hypertimeshift'+shot,'datstretch datum=${SOURCES[1]}')
Flow('warpedhet'+shot,'zo'+shot+' hyperhettimeshift'+shot,'datstretch datum=${SOURCES[1]}')
Plot('warped'+shot,
'''
grey transp=y yreverse=y poly=y
title="Flattened at %s km (1D)" min1=1.65 max1=1.95 label2=Midpoint unit2=km screenratio=0.5
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=10
wherexlabel=top
''' % str(j*0.0101) )
Plot('warpedhet'+shot,
'''
grey transp=y yreverse=y poly=y
title="Flattened at %s km (Proposed)" min1=1.65 max1=1.95 label2=Midpoint unit2=km screenratio=0.5
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=10
wherexlabel=top
''' % str(j*0.0101) )
Plot('t0to-'+shot,'t0sumext offset',
'''
window n1=1 f1=%d n2=1 f2=4 | scale dscale=2 | cat axis=2 ${SOURCES[1]} order=2,1 | transp | dd type=complex | window |
graph yreverse=y plotcol=5 plotfat=16 min2=1.65 max2=1.95 screenratio=0.5 dash=1
axisfat=3 wanttitle=n wantaxis=n
''' % j )
Plot('warpedn'+shot,'warped'+shot+' t0to-'+shot,'Overlay')
Plot('warpedhetn'+shot,'warpedhet'+shot+' t0to-'+shot,'Overlay')
Result('warpedhypercompare3','warpedn3 warpedhetn3','TwoRows')
Flow('vold','vnmosq','math output="sqrt(input)" | window n1=1 f1=426')
vold = '2.477,2.973,3.12,3.222,3.382'
Flow('vnew','vnmosqhet','math output="sqrt(input)" | window n1=1 f1=426')
vnew = '2.477,2.974,3.165,3.336,3.693'
Flow('refdiff','t0sum','window n1=1 f1=425 | add scale=2 | spray axis=2 n=201 o=3.2825 d=0.0101 | transp')
Flow('vdiffmodelold','refdiff','unif2 d1=0.004 n1=751 v00=%s | smooth rect1=10' %vold)
Flow('vdiffmodelnew','refdiff','unif2 d1=0.004 n1=751 v00=%s | smooth rect1=10' %vnew)
Flow('oldmig','zotaper vdiffmodelold',' window | mig2 antialias=1.0 vel=${SOURCES[1]}')
Flow('newmig','zotaper vdiffmodelnew',' window | mig2 antialias=1.0 vel=${SOURCES[1]}')
Plot('oldmig',
'''
grey transp=y yreverse=y poly=y pclip=99
title="Migrated diffraction (1D)" min1=0 max1=2 label2=Midpoint unit2=km screenratio=0.75
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=12
wherexlabel=top
''' )
Plot('newmig',
'''
grey transp=y yreverse=y poly=y pclip=99
title="Migrated diffraction (Proposed)" min1=0 max1=2 label2=Midpoint unit2=km screenratio=0.75
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=12
wherexlabel=top
''' )
Result('migcompare','oldmig newmig','SideBySideAniso')
Flow('oldmigfocus','oldmig','window min1=1.65 max1=1.9 | focus dim=2 rect1=10 rect2=3')
Flow('newmigfocus','newmig','window min1=1.65 max1=1.9 | focus dim=2 rect1=10 rect2=3')
Plot('oldmigfocus',
'''
grey transp=y yreverse=y poly=y pclip=99
title="Focus of migrated diffraction (1D)" min1=1.65 max1=1.9 label2=Midpoint unit2=km screenratio=0.55
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=10
wherexlabel=top color=j allpos=y
''' )
Plot('newmigfocus',
'''
grey transp=y yreverse=y poly=y pclip=99
title="Focus of migrated diffraction (Proposed)" min1=1.65 max1=1.9 label2=Midpoint unit2=km screenratio=0.55
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=10
wherexlabel=top color=j allpos=y
''' )
Result('focuscompare','oldmigfocus newmigfocus','TwoRows')
Flow('oldmigenv','oldmig','window min1=1.65 max1=1.9 | envelope')
Flow('newmigenv','newmig','window min1=1.65 max1=1.9 | envelope')
Plot('oldmigenv',
'''
grey transp=y yreverse=y poly=y
title="Focus of migrated diffraction (1D)" min1=1.65 max1=1.9 label2=Midpoint unit2=km screenratio=0.55
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=10
wherexlabel=top color=j allpos=y scalebar=y maxval=0.077 clip=0.077
''' )
Plot('newmigenv',
'''
grey transp=y yreverse=y poly=y
title="Focus of migrated diffraction (Proposed)" min1=1.65 max1=1.9 label2=Midpoint unit2=km screenratio=0.55
axisfat=3 titlefat=6 titlesz=14 labelfat=6 labelsz=10
wherexlabel=top color=j allpos=y scalebar=y maxval=0.077 clip=0.077
''' )
Result('envcompare','oldmigenv newmigenv','TwoRows')
updown2=[1,2,3,4,5]
N2 = len(updown2)-1
xinitial = [4.,4.,4.,4.]
updown2str = ','.join(map(str,updown2))
xinistr = ','.join(map(str,xinitial))
Flow('ray-2','refcut',
'''isaac2 niter=10 number=%d vstatus=%d debug=n xs=4.3026 xr=4.2426 velocity=%s layer=%s
xref=%s zref=%s xgradient=%s zgradient=%s xinitial=%s tol=1e-5''' % (N2,vstat,vstr,updown2str,xrefstr,zrefstr,gxstr,gzstr,xinistr))
Flow('ray-1','refcut',
'''isaac2 niter=10 number=%d vstatus=%d debug=n xs=4.3026 xr=4.2626 velocity=%s layer=%s
xref=%s zref=%s xgradient=%s zgradient=%s xinitial=%s tol=1e-5''' % (N2,vstat,vstr,updown2str,xrefstr,zrefstr,gxstr,gzstr,xinistr))
Flow('ray0','refcut',
'''isaac2 niter=10 number=%d vstatus=%d debug=n xs=4.3026 xr=4.2826 velocity=%s layer=%s
xref=%s zref=%s xgradient=%s zgradient=%s xinitial=%s tol=1e-5''' % (N2,vstat,vstr,updown2str,xrefstr,zrefstr,gxstr,gzstr,xinistr))
Flow('ray+1','refcut',
'''isaac2 niter=10 number=%d vstatus=%d debug=n xs=4.3026 xr=4.3026 velocity=%s layer=%s
xref=%s zref=%s xgradient=%s zgradient=%s xinitial=%s tol=1e-5''' % (N2,vstat,vstr,updown2str,xrefstr,zrefstr,gxstr,gzstr,xinistr))
Flow('ray+2','refcut',
'''isaac2 niter=10 number=%d vstatus=%d debug=n xs=4.3026 xr=4.3226 velocity=%s layer=%s
xref=%s zref=%s xgradient=%s zgradient=%s xinitial=%s tol=1e-5''' % (N2,vstat,vstr,updown2str,xrefstr,zrefstr,gxstr,gzstr,xinistr))
for k in ['-2','-1','0','+1','+2']:
Plot('ray'+k,
'''
dd type=complex | window |
graph wanttitle=n wantaxis=n yreverse=y min2=0 max2=3 min1=0 max1=6 plotcol=7 plotfat=6 screenht=5.0 screenratio=0.333 yll=3.5 xll=1.5
axisfat=3 titlefat=3 titlesz=10 labelfat=3 labelsz=6
''')
Result('dt2dx2','ray-2 ray-1 ray0 ray+1 ray+2 ref','Overlay')
End()
|