```# Low frequency scalar wave scattering attenuation and penetration in 3-D isotropic von Karman medium # # December 2007 # # Thomas Jules Browaeys # Bureau of Economic Geology # University of Texas at Austin # mailto:jules.browaeys@beg.utexas.edu from rsf.proj import * from math import pi # ---------------------------------------------- # FREQUENCY 1-D GRID # ---------------------------------------------- # nx = points number # dx = frequency steps pgrid = {'nx':300, 'ox':1., 'dx':1.} # Create frequency grid Flow('fgrd',None,'spike nsp=1 mag=1 n1=%(nx)d d1=%(dx)g o1=%(ox)g | put label1=f unit1=Hz' % pgrid) # Logarithmic scale Flow('lfgrd','fgrd','math output="log(x1)/log(10.)"') # ---------------------------------------------- # STOCHASTIC MEDIUM # ---------------------------------------------- # hu = H exponent > -0.5 # cu = H spectrum constant # Exponents and constants pgrids = {'M025':{'hu':-0.25, 'cu':1.311, 'shu':'\F15 8\F2 0.25'}, '025':{'hu':0.25, 'cu':0.5991, 'shu':'0.25'}, '05':{'hu':0.5, 'cu':1.0, 'shu':'0.5'}, '075':{'hu':0.75, 'cu':1.3317, 'shu':'0.75'}, '1':{'hu':1.0, 'cu':1.57, 'shu':'1.0'} } # Correlation lengths (m) bvalues = {'b25':2.5, 'b05':5.0, 'b10':10.0} bvalstr = {'b25':'2.5', 'b05':'5', 'b10':'10'} # VP VS variations and velocity parq = { 'sp':0.3, # P wave standard deviation 'ss':0.3, # S wave standard deviation 'cp':2700., # P wave velocity 'cs':1230. # S wave velocity } # Variances parq['vrs2'] = parq['ss']*parq['ss'] parq['vrp2'] = parq['sp']*parq['sp'] for lb in bvalues.keys(): parq['b'] = bvalues[lb] bstring = bvalstr[lb] ksb = 'ksb' + lb ksb2 = 'ksb2' + lb kpb = 'kpb' + lb kpb2 = 'kpb2' + lb # Wave vector k*b for VS Flow(ksb,'fgrd','math output="(%g)*x1/(%g)"' % (2.*pi*parq['b'],parq['cs'])) # Squared (k*b) for VS Flow(ksb2,ksb,'math output="input*input"') # Wave vector k*b for VP Flow(kpb,'fgrd','math output="(%g)*x1/(%g)"' % (2.*pi*parq['b'],parq['cp'])) # Squared (k*b) for VP Flow(kpb2,kpb,'math output="input*input"') # ----------------------- # Exponential correlation # ----------------------- # S waves dvse = 'dvse' + lb qse = 'qse' + lb depthse = 'depthse' + lb # Dispersion velocity Flow(dvse,ksb2,'math k2=\${SOURCES[0]} output="1./(1.+(%g)*(0.5+4.*k2)/(1.+4.*k2))"' % (parq['vrs2'])) # Attenuation 1/Q Flow(qse,[ksb2, ksb],'math k2=\${SOURCES[0]} k=\${SOURCES[1]} output="(%g)*(k2*k)/(1.+4.*k2)"' % (8.*parq['vrs2'])) # Depth penetration Flow(depthse,ksb2,'math k2=\${SOURCES[0]} output="(%g)/(%g)*(1.+4.*k2)/(k2*k2)"' % (parq['b'],4.*parq['vrs2'])) # P waves dvpe = 'dvpe' + lb qpe = 'qpe' + lb depthpe = 'depthpe' + lb # Dispersion velocity Flow(dvpe,kpb2,'math k2=\${SOURCES[0]} output="1./(1.+(%g)*(0.5+4.*k2)/(1.+4.*k2))"' % (parq['vrp2'])) # Attenuation 1/Q Flow(qpe,[kpb2, kpb],'math k2=\${SOURCES[0]} k=\${SOURCES[1]} output="(%g)*(k2*k)/(1.+4.*k2)"' % (8.*parq['vrp2'])) # Depth penetration Flow(depthpe,kpb2,'math k2=\${SOURCES[0]} output="(%g)/(%g)*(1.+4.*k2)/(k2*k2)"' % (parq['b'],4.*parq['vrp2'])) lfdvse = 'lfdvse' + lb lfdvpe = 'lfdvpe' + lb lqse = 'lqse' + lb lqpe = 'lqpe' + lb lfqse = 'lfqse' + lb lfqpe = 'lfqpe' + lb lldve = 'lldve' + lb llqe = 'llqe' + lb depthle = 'depthle' + lb # Plots Flow(lfdvse,['lfgrd', dvse],'cmplx \${SOURCES[0:2]}',stdin=0) Flow(lfdvpe,['lfgrd', dvpe],'cmplx \${SOURCES[0:2]}',stdin=0) Flow(lqse,qse,'math output="log(input)/log(10.)"') Flow(lqpe,qpe,'math output="log(input)/log(10.)"') Flow(lfqse,['lfgrd', lqse],'cmplx \${SOURCES[0:2]}',stdin=0) Flow(lfqpe,['lfgrd', lqpe],'cmplx \${SOURCES[0:2]}',stdin=0) # Plot velocity dispersion Result(lldve,[lfdvpe, lfdvse], ''' cat \${SOURCES[0:2]} axis=2 | put label1='log[f]' unit1='Hz' label2='c/c\_0\^' unit2= | graph dash=0,3 title="Dispersion velocity - H=0.5" min1=0. max1=2.5 min2=0.9 max2=1. ''',stdin=0) # Plot attenuation Result(llqe,[lfqpe, lfqse], ''' cat \${SOURCES[0:2]} axis=2 | put label1='log[f]' unit1='Hz' label2='log[1/Q]' unit2= | graph dash=0,3 title="Attenuation - H=0.5" min1=0. max1=2.5 max2=0. min2=-6. ''',stdin=0) # Plot depth penetration Result(depthle,[depthpe, depthse], ''' cat \${SOURCES[0:2]} axis=2 | put label2='d' unit2='m' | graph dash=0,3 title="Penetration depth - H=0.5" min1=0. max1=100. min2=0. max2=2500. ''',stdin=0) # ------------------- # Fractal correlation # ------------------- for fexp in ('M025','025','075'): parq.update(pgrids[fexp]) # S waves hksb4 = 'hksb4' + lb + fexp qsf = 'qsf' + lb + fexp depthsf = 'depthsf' + lb + fexp Flow(hksb4,[ksb2, ksb],'math k2=\${SOURCES[0]} k=\${SOURCES[1]} output="k-k/(1.+4.*k2)^(%g)"' % (0.5+parq['hu'])) # Attenuation 1/Q Flow(qsf,hksb4,'math hk4=\${SOURCES[0]} output="2.*(%g)*hk4"' % (parq['vrs2']*parq['cu'])) # Depth penetration Flow(depthsf,[hksb4, ksb],'math hk4=\${SOURCES[0]} k=\${SOURCES[1]} output="(%g)/(%g)/(k*hk4)"' % (parq['b'],parq['vrs2']*parq['cu'])) # P waves hkpb4 = 'hkpb4' + lb + fexp qpf = 'qpf' + lb + fexp depthpf = 'depthpf' + lb + fexp Flow(hkpb4,[kpb2, kpb],'math k2=\${SOURCES[0]} k=\${SOURCES[1]} output="k-k/(1.+4.*k2)^(%g)"' % (0.5+parq['hu'])) # Attenuation 1/Q Flow(qpf,hkpb4,'math hk4=\${SOURCES[0]} output="2.*(%g)*hk4"' % (parq['vrp2']*parq['cu'])) # Depth penetration Flow(depthpf,[hkpb4, kpb],'math hk4=\${SOURCES[0]} k=\${SOURCES[1]} output="(%g)/(%g)/(k*hk4)"' % (parq['b'],parq['vrp2']*parq['cu'])) lqsf = 'lqsf' + lb + fexp lqpf = 'lqpf' + lb + fexp lfqsf = 'lfqsf' + lb + fexp lfqpf = 'lfqpf' + lb + fexp llqf = 'llqf' + lb + fexp depthlf = 'depthlf' + lb + fexp # Plots Flow(lqsf,qsf,'math output="log(input)/log(10.)"') Flow(lqpf,qpf,'math output="log(input)/log(10.)"') Flow(lfqsf,['lfgrd', lqsf],'cmplx \${SOURCES[0:2]}',stdin=0) Flow(lfqpf,['lfgrd', lqpf],'cmplx \${SOURCES[0:2]}',stdin=0) # Plot attenuation titleatt = 'title="Attenuation \k100 \k100 \F5 b\F2 \k70 =\k70 ' + bstring + '\k10 m\g \G \F5 H\F2 \k60 =\k60 ' + parq['shu'] + '"' Result(llqf,[lfqpf, lfqsf], ''' cat \${SOURCES[0:2]} axis=2 | put label1='log[f]' unit1='Hz' label2='log[1/Q]' unit2= | graph dash=0,3 min1=0. max1=2.5 max2=0. min2=-6. font=2 parallel2=n plotfat=4,4 n2tic=7 o2num=-6.0 d2num=1.0 %s ''' % (titleatt),stdin=0) # Plot depth penetration titledep = 'title="Penetration depth \k100 \k100 \F5 b\F2 \k70 =\k70 ' + bstring + '\k10 m\g \G \F5 H\F2 \k60 =\k60 ' + parq['shu'] + '"' Result(depthlf,[depthpf, depthsf], ''' cat \${SOURCES[0:2]} axis=2 | put label2='d' unit2='m' | graph dash=0,3 min1=0. max1=100. min2=0. max2=2500. font=2 parallel2=n plotfat=4,4 n2tic=6 o2num=0.0 d2num=500.0 %s ''' % (titledep),stdin=0) End()```