```from rsf.proj import * import sys sys.path.append('Python') import fdmod,stiff,pot,fdd,pot,spk,pplot from math import * def arr2str(array,sep=' '): return string.join(map(str,array),sep) par = { 'nx':600, 'ox':0, 'dx':0.001, 'lx':'x', 'ux':'km', 'nz':600, 'oz':0, 'dz':0.001, 'lz':'z', 'uz':'km', 'nt':1601,'ot':0, 'dt':0.0001, 'lt':'t', 'ut':'s', 'kt':150, 'jsnap':200, 'height':10, 'nb':0, 'frq':100, 'ratio':1,'wweight':100 } fdmod.param(par) par['labelattr']=par['labelattr']+''' titlesz=20 labelsz=10 bartype=h font=4 barlabelsz=10 n2tic=6 o2num=0 d2num=.1 ''' xmax = par['dx']*(par['nx']-1) zmax = par['dz']*(par['nz']-1) x0=xmax/2 # -C*(x-x0)^2+z0 C=1.0 z0=0.2 z1=z0-C*(0*xmax/3-x0)*(0*xmax/3-x0) z2=z0-C*(1*xmax/3-x0)*(1*xmax/3-x0) z3=z0-C*(2*xmax/3-x0)*(2*xmax/3-x0) z4=z0-C*(3*xmax/3-x0)*(3*xmax/3-x0) dz1=.1 dz2=.2 dz3=.4 #z1=0.015 #z2=.011 #z3=-.001 #z4=-.021 # #dz1=.01 #dz2=.03 #dz3=.04 layers = ((z1, z2, z3, z4), (z1+dz1,z2+dz1,z3+dz1,z4+dz1), (z1+dz2,z2+dz2,z3+dz2,z4+dz2), (z1+dz3,z2+dz3,z3+dz3,z4+dz3)) velocities = (2.00, 2.20, 2.35, 2.50, 3.25) #vs=velocities/2 nu0=2*C*x0 nu=(nu0,nu0,nu0,nu0,nu0) n1 = len(layers[0]) n2 = len(layers) Flow('layers',None, ''' echo %s n1=%d n2=%d o1=0 d1=%g data_format=ascii_float in=\$TARGET ''' % (string.join(map(arr2str,layers),' '), n1,n2,xmax/(n1-1))) Flow('vp','layers', ''' spline o1=0 d1=%g n1=%d | unif2 d1=%g n1=%d v00=%s| dd type=float form=native ''' % (par['dz'],int(1.5+xmax/par['dz']), par['dz'],int(1.5+zmax/par['dz']), arr2str(velocities,','),)) Flow('zero','','spike o1=%(oz)g o2=%(ox)g d1=%(dz)g d2=%(dx)g n1=%(nz)d n2=%(nx)d|scale rscale=0'%par) #Result('vp', # ''' # grey color=j title="Model 1" # screenratio=%g wantscalebar=y # allpos=y titlesz=8 labelsz=6 # label1="Depth (km)" # label2="Distance (km)" # ''' % (zmax/xmax)) rate=-2*C Flow('nu','layers', ''' spline o1=0 d1=%g n1=%d | unif2 dvdx=%g d1=%g n1=%d v00=%s | dd type=float form=native | math output="atan(input)*180/3.14159265" ''' % (par['dz'],int(1.5+xmax/par['dz']),rate, par['dz'],int(1.5+zmax/par['dz']), arr2str(nu,','), )) #Flow('nu','','spike o1=%(oz)g o2=%(ox)g d1=%(dz)g d2=%(dx)g n1=%(nz)d n2=%(nx)d|scale rscale=45'%par) #Result('nu', # ''' # grey color=j title="Model 1" # screenratio=%g wantscalebar=y # titlesz=8 labelsz=6 bias=0 # label1="Depth (km)" # label2="Distance (km)" # ''' % (zmax/xmax)) Flow('vs', 'vp','math output="input/2" ') Flow('ro', 'vp','math output="input-0.6" ') Flow('epsilon','vp','math output="input/4-0.4" ') Flow('delta', 'vp','math output="input/8-0.2" ') #Flow('nu', 'vp','math output=45 ') #Plot('vs',fdmod.cgrey('bias=1 allpos=y',par)) #Plot('ro',fdmod.cgrey('bias=1 allpos=y',par)) #Plot('epsilon',fdmod.cgrey('',par)) #Plot('delta',fdmod.cgrey('bias=-0.29',par)) # ------------------------------------------------------------ # source/receiver coordinates # ------------------------------------------------------------ fdmod.point('ss', par['ox']+(par['nx']/2*par['dx']), par['oz']+(par['nz']/6*par['dz']),par) fdmod.horizontal('rr',0,par) Plot('rr',fdmod.rrplot('screenratio=1',par)) Plot('ss',fdmod.ssplot('yll=2.15 screenht=10.15 screenratio=.925',par)) # ------------------------------------------------------------ barattr=' xll=2 yll=1.3 wherebartic=top wherebarlabel=top' Plot('vp',fdmod.cgrey('color=j wantscalebar=y allpos=y barlabel="V\_P0\^ (km/s)" '+barattr,par)) Plot('vs',fdmod.cgrey('color=j wantscalebar=y allpos=y barlabel="V\_S0\^ (km/s)"'+barattr ,par)) Plot('ro',fdmod.cgrey('color=j wantscalebar=y allpos=y barlabel="\F10 r \F3 (g/cm\^3\_)" '+barattr,par)) Plot('epsilon',fdmod.cgrey('color=j wantscalebar=y allpos=y barlabel="\s140 \F10 e" formatbar=%4.2f'+barattr,par)) Plot('delta',fdmod.cgrey('color=j wantscalebar=y allpos=y barlabel="\s140 \F10 d" formatbar=%4.2f'+barattr,par)) Plot('nu',fdmod.cgrey('color=e wantscalebar=y formatbar=%2.0f barlabel="\s140 \F10 n \s100 (\^o\_)" '+barattr,par)) for k in (['vp', 'vs', 'ro', 'epsilon', 'delta', 'nu']): Result(k,[k],'Overlay') # ------------------------------------------------------------ # elastic source # ------------------------------------------------------------ fdmod.wavelet('wav_',par['frq'],par) Flow('ver','wav_','math output="1*input"') Flow('hor','wav_','math output="1*input"') Flow('wave',['ver','hor'], ''' cat axis=2 space=n \${SOURCES[1:2]} | transp plane=12 | transp plane=23 | transp plane=12 ''') fdmod.ewavelet('wave','',par) stiff.tti2d('cA','vp','vs','ro','epsilon','delta','nu',par) Flow( ['dAu','uA'],['wave','cA','ro','ss','rr'], ''' sfewefd2dtti ompchunk=1 ompnth=0 verb=y free=n snap=y jsnap=200 nb=0 ssou=n opot=n nbell=5 anitype=t ccc=\${SOURCES[1]} den=\${SOURCES[2]} sou=\${SOURCES[3]} rec=\${SOURCES[4]} wfl=\${TARGETS[1]} ''' % par) pot.displacementsC('uA','uA','uAz','uAx',7,'',par) pot.cliptogether( 'uA','uAz','uAx','"\F5 W\_z"','"\F5 W\_x"',1,'pclip=98 labelsz=12 yll=1.7 xll=2.3',par) nframe=(par['nt']-1)//par['jsnap']+1 Flow('uAc','uA','window n1=%(nz)d n2=%(nx)d f1=%(nbell)d f2=%(nbell)d |scale axis=123|scale rscale=10'%par) Flow('uAwom',['vp','uAc'],'add add=-2.4 | spray axis=3 n=2|spray axis=4 n=%d | scale axis=123|add \${SOURCES[1]}'%nframe) fdmod.emovie('uAmovie','uAwom',nframe,'pclip=98',2,par) fdd.derivatives(par) order=8 #separate in X domain spk.delt('spk',64,64,par) fdd.separatorD('dzX','dxX','spk','cA','n','x','gaussian',1.0,8,27,27,par) #fdd.separator('dzX','dxX','spk','cA','nu','n','x',27,27,par) fdd.oneirST('xop','dzX-tmp','dxX-tmp',27,27,'color=j pclip=100',par) pot.potentials( 'pA','uAz','uAx','dzX','dxX','n','','q',par) pot.cliptogether( 'pA','pAp','pAs','"\F5 qP"','"\F5 qS"',1,'pclip=98 labelsz=12 yll=1.7 xll=2.3',par) Flow('uAx-n','uAx','scale axis=123| noise range=.2') Flow('uAz-n','uAz','scale axis=123| noise range=.1') #pot.cliptogether( 'uA-n','uAz-n','uAx-n','"u\_z"','"u\_x"',1,'pclip=98',par) #pot.potentials( 'pA-n','uAz-n','uAx-n','dzX','dxX','n','','q',par) #pot.cliptogether( 'pA-n','pA-np','pA-ns',"qP","qS",1,'pclip=98',par) #fdd.oneirNS('aop','dzX','dxX',7,7,'wheretitle=top pclip=100 wantaxis=n titlesz=30 color=e',par) fdd.oneirNS_rot('rop','dzX','dxX','nu',7,7,'wheretitle=top pclip=100 wantaxis=n wanttitle=y titlesz=30 color=e ',par) fdmod.boxarray('aoppos',3,(par['nx']/4+1)*par['dx'],par['nx']/4*par['dx'], 3,(par['nz']/4+1)*par['dz'],par['nz']/4*par['dz'],par) Plot('aoppos',fdmod.qqplot('''symbol=o plotcol=6 plotfat=10 wantscalebar=y allpos=y barlabel="V\_P0\^ (km/s)" xll=2 yll=1.3''',par)) Result('aoppos',['nu','aoppos'],'Overlay') allplots=['rop00','rop01','rop02', 'rop10','rop11','rop12', 'rop20','rop21','rop22'] #allplots=['aop00','aop01','aop02', # 'aop10','aop11','aop12', # 'aop20','aop21','aop22'] pplot.multip('rop',allplots,3,3,0.25,0.25,-10,-16) #separate with isotropic operators stiff.tti2d('cI','vp','vs','ro','zero','zero','zero',par) fdd.separatorD('dzI','dxI','spk','cI','y','x','gaussian',1.0,8,27,27,par) #fdd.separator('dzI','dxI','spk','cI','zero','y','x',27,27,par) pot.potentials( 'pI','uAz','uAx','dzI','dxI','y','','q',par) pot.cliptogether( 'pI','pIp','pIs','"\F5 qP"','"\F5 qS"',1,'pclip=98 labelsz=12 yll=1.7 xll=2.3',par) #separate with VTI operators stiff.tti2d('cV','vp','vs','ro','epsilon','delta','zero',par) fdd.separatorD('dzV','dxV','spk','cV','n','x','gaussian',1.0,8,27,27,par) #fdd.separator('dzV','dxV','spk','cV','zero','n','x',27,27,par) pot.potentials( 'pV','uAz','uAx','dzV','dxV','n','','q',par) pot.cliptogether( 'pV','pVp','pVs','"\F5 qP"','"\F5 qS"',1,'pclip=98 labelsz=12 yll=1.7 xll=2.3',par) End()```