```############################################################## #### Interferometry/Imaging/Exponential -- Random Model for validation #### of acoustic scattering reciprocity. #### Author: Ivan Vasconcelos ############################################################## from rsf.proj import * import fdmod,gfield,rules par = { 'nx':1001, 'ox':0, 'dx':0.002, 'lx':'x', 'ux':'km', 'nz':501, 'oz':0, 'dz':0.002, 'lz':'z', 'uz':'km', 'nt':5001, 'ot':0, 'dt':0.0002, 'lt':'t', 'ut':'s', 'kt':250, 'jsnap':200, 'nb':200, 'frq':30, 'ssou':'y', } fdmod.param(par) #################################################################### ### Setting up the model .... #################################################################### # coordinates for some model features par['kf']=0.624/par['dz'] + 1 # = 0.624 par['kg']=1.0/par['dz'] + 1 # = 0.824 par['kz']=0.770/par['dz'] + 1 # = 0.598 par['kh']=0.448/par['dz'] + 1 # = 0.448 par['oxf']=0.5/par['dx'] + 1 # = 0.500 #par['kf']=1.25/2.*par['nz'] #par['kg']=1.65/2.*par['nz'] #par['kz']=0.60*par['nz'] #par['kh']=0.90/2.*par['nz'] #par['oxf']=0.5*par['nx'] par['oxda']=0.742/par['dx'] + 1 par['kxda']=0.752/par['dx'] + 1 par['ozda']=0.492/par['dz'] + 1 par['kzda']=0.504/par['dz'] + 1 #par['oxda']=122 0.242/par['dx'] + 1 #par['kxda']=127 #par['ozda']=247 #par['kzda']=252 par['oxdb']=1.244/par['dx'] + 1 par['kxdb']=1.252/par['dx'] + 1 par['ozdb']=0.492/par['dz'] + 1 par['kzdb']=0.504/par['dz'] + 1 #par['oxdb']=373 #par['kxdb']=377 #par['ozdb']=247 #par['kzdb']=252 par['oxdc']=0.996/par['dx'] + 1 par['kxdc']=1.004/par['dx'] + 1 par['ozdc']=0.352/par['dx'] + 1 par['kzdc']=0.364/par['dx'] + 1 #par['oxdc']=249 #par['kxdc']=253 #par['ozdc']=177 #par['kzdc']=182 #par['frq']=45 par['ssou']='n' par['vscale']=0 par['hscale']=1 # ------------------------------------------------------------ # Make models using sfspike # P velocity (km/s) Flow('vpb',None, ''' spike nsp=6 mag=0.25,0.6,0.35,2,2,2 p2=-0.15,0.25,0.0,0.0,0.0,0.0 n1=%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d o1=%(oz)g,%(oz)d,%(oz)d,%(oz)d,%(oz)d,%(oz)d d1=%(dz)g,%(dz)d,%(dz)d,%(nz)d,%(nz)d,%(nz)d k1=%(kz)d,%(kf)d,%(kh)d,%(ozda)d,%(ozdb)d,%(ozdc)d l1=%(kg)d,%(nz)d,%(kf)d,%(kzda)d,%(kzdb)d,%(kzdc)d n2=%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d o2=%(ox)g,%(ox)d,%(ox)d,%(ox)d,%(ox)d,%(ox)d d2=%(dx)g,%(dx)d,%(dx)d,%(dx)d,%(dx)d,%(dx)d k2=%(ox)d,%(ox)d,%(ox)d,%(oxda)d,%(oxdb)d,%(oxdc)d l2=%(nx)d,%(nx)d,%(nx)d,%(kxda)d,%(kxdb)d,%(kxdc)d | add add=3.0 | put label1=z label2=x unit1=km unit2=km ''' % par) # Flow('vpbo',None, ''' spike nsp=6 mag=0,0,0,0,0,0 p2=-0.15,0.25,0.0,0.0,0.0,0.0 n1=%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d o1=%(oz)g,%(oz)d,%(oz)d,%(oz)d,%(oz)d,%(oz)d d1=%(dz)g,%(dz)d,%(dz)d,%(nz)d,%(nz)d,%(nz)d k1=%(kz)d,%(kf)d,%(kh)d,%(ozda)d,%(ozdb)d,%(ozdc)d l1=%(kg)d,%(nz)d,%(kf)d,%(kzda)d,%(kzdb)d,%(kzdc)d n2=%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d o2=%(ox)g,%(ox)d,%(ox)d,%(ox)d,%(ox)d,%(ox)d d2=%(dx)g,%(dx)d,%(dx)d,%(dx)d,%(dx)d,%(dx)d k2=%(ox)d,%(ox)d,%(ox)d,%(oxda)d,%(oxdb)d,%(oxdc)d l2=%(nx)d,%(nx)d,%(nx)d,%(kxda)d,%(kxdb)d,%(kxdc)d | add add=3.0 | put label1=z label2=x unit1=km unit2=km ''' % par) # S velocity (km/s) Flow('vsb',None, ''' spike nsp=6 mag=1.0,1.0,1.0,2.5,2.5,2.5 p2=-0.15,0.25,0.0,0.0,0.0,0.0 n1=%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d o1=%(oz)g,%(oz)d,%(oz)d,%(oz)d,%(oz)d,%(oz)d d1=%(dz)g,%(dz)d,%(dz)d,%(nz)d,%(nz)d,%(nz)d k1=%(kz)d,%(kf)d,%(kh)d,%(ozda)d,%(ozdb)d,%(ozdc)d l1=%(kg)d,%(nz)d,%(kf)d,%(kzda)d,%(kzdb)d,%(kzdc)d n2=%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d o2=%(ox)g,%(ox)d,%(ox)d,%(ox)d,%(ox)d,%(ox)d d2=%(dx)g,%(dx)d,%(dx)d,%(dx)d,%(dx)d,%(dx)d k2=%(ox)d,%(ox)d,%(ox)d,%(oxda)d,%(oxdb)d,%(oxdc)d l2=%(nx)d,%(nx)d,%(nx)d,%(kxda)d,%(kxdb)d,%(kxdc)d | add add=0.00001 | put label1=z label2=x unit1=km unit2=km ''' % par) # Flow('vsbo',None, ''' spike nsp=6 mag=0,0,0,0,0,0 p2=-0.15,0.25,0.0,0.0,0.0,0.0 n1=%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d o1=%(oz)g,%(oz)d,%(oz)d,%(oz)d,%(oz)d,%(oz)d d1=%(dz)g,%(dz)d,%(dz)d,%(nz)d,%(nz)d,%(nz)d k1=%(kz)d,%(kf)d,%(kh)d,%(ozda)d,%(ozdb)d,%(ozdc)d l1=%(kg)d,%(nz)d,%(kf)d,%(kzda)d,%(kzdb)d,%(kzdc)d n2=%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d o2=%(ox)g,%(ox)d,%(ox)d,%(ox)d,%(ox)d,%(ox)d d2=%(dx)g,%(dx)d,%(dx)d,%(dx)d,%(dx)d,%(dx)d k2=%(ox)d,%(ox)d,%(ox)d,%(oxda)d,%(oxdb)d,%(oxdc)d l2=%(nx)d,%(nx)d,%(nx)d,%(kxda)d,%(kxdb)d,%(kxdc)d | add add=0.00001 | put label1=z label2=x unit1=km unit2=km ''' % par) # Density (kg/km^3) Flow('rob',None, ''' spike nsp=6 mag=500000,500000,500000,1000000,1000000,1000000 p2=-0.15,0.25,0.0,0.0,0.0,0.0 n1=%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d o1=%(oz)g,%(oz)d,%(oz)d,%(oz)d,%(oz)d,%(oz)d d1=%(dz)g,%(dz)d,%(dz)d,%(nz)d,%(nz)d,%(nz)d k1=%(kz)d,%(kf)d,%(kh)d,%(ozda)d,%(ozdb)d,%(ozdc)d l1=%(kg)d,%(nz)d,%(kf)d,%(kzda)d,%(kzdb)d,%(kzdc)d n2=%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d o2=%(ox)g,%(ox)d,%(ox)d,%(ox)d,%(ox)d,%(ox)d d2=%(dx)g,%(dx)d,%(dx)d,%(dx)d,%(dx)d,%(dx)d k2=%(ox)d,%(ox)d,%(ox)d,%(oxda)d,%(oxdb)d,%(oxdc)d l2=%(nx)d,%(nx)d,%(nx)d,%(kxda)d,%(kxdb)d,%(kxdc)d | add add=2000000 | put label1=z label2=x unit1=km unit2=km ''' % par) # Flow('robo',None, ''' spike nsp=6 mag=0,0,0,0,0,0 p2=-0.15,0.25,0.0,0.0,0.0,0.0 n1=%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d,%(nz)d o1=%(oz)g,%(oz)d,%(oz)d,%(oz)d,%(oz)d,%(oz)d d1=%(dz)g,%(dz)d,%(dz)d,%(nz)d,%(nz)d,%(nz)d k1=%(kz)d,%(kf)d,%(kh)d,%(ozda)d,%(ozdb)d,%(ozdc)d l1=%(kg)d,%(nz)d,%(kf)d,%(kzda)d,%(kzdb)d,%(kzdc)d n2=%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d,%(nx)d o2=%(ox)g,%(ox)d,%(ox)d,%(ox)d,%(ox)d,%(ox)d d2=%(dx)g,%(dx)d,%(dx)d,%(dx)d,%(dx)d,%(dx)d k2=%(ox)d,%(ox)d,%(ox)d,%(oxda)d,%(oxdb)d,%(oxdc)d l2=%(nx)d,%(nx)d,%(nx)d,%(kxda)d,%(kxdb)d,%(kxdc)d | add add=2000000 | put label1=z label2=x unit1=km unit2=km ''' % par) # Thomsen anisotropy parameters (set to zero in these examples) # Epsilon Flow('epsilon',None, ''' spike nsp=1 mag=0 n1=%(nz)d o1=%(oz)g d1=%(dz)g n2=%(nx)d o2=%(ox)g d2=%(dx)g k1=%(kz)d l1=%(nz)d | add add=0 | put label1=z label2=x unit1=km unit2=km ''' % par) # Delta Flow('delta',None, ''' spike nsp=1 mag=0 n1=%(nz)d o1=%(oz)g d1=%(dz)g n2=%(nx)d o2=%(ox)g d2=%(dx)g k1=%(kz)d l1=%(nz)d | add add=0 | put label1=z label2=x unit1=km unit2=km ''' % par) # ------------------------------------------------------------ # random field - generating random model perturbation # The ACTUAL model for this exercise is made here! # background velocities # P-wave Flow('vpo',None, ''' spike nsp=1 mag=0 n1=%(nz)d o1=%(oz)g d1=%(dz)g n2=%(nx)d o2=%(ox)g d2=%(dx)g k1=%(kz)d l1=%(nz)d | add add=3.0 | put label1=z label2=x unit1=km unit2=km ''' % par) # S-wave Flow('vso',None, ''' spike nsp=1 mag=0 n1=%(nz)d o1=%(oz)g d1=%(dz)g n2=%(nx)d o2=%(ox)g d2=%(dx)g k1=%(kz)d l1=%(nz)d | add add=0.00001 | put label1=z label2=x unit1=km unit2=km ''' % par) # Density Flow('roo',None, ''' spike nsp=1 mag=0 n1=%(nz)d o1=%(oz)g d1=%(dz)g n2=%(nx)d o2=%(ox)g d2=%(dx)g k1=%(kz)d l1=%(nz)d | add add=2000000 | put label1=z label2=x unit1=km unit2=km ''' % par) # Epsilon Flow('epso',None, ''' spike nsp=1 mag=0 n1=%(nz)d o1=%(oz)g d1=%(dz)g n2=%(nx)d o2=%(ox)g d2=%(dx)g k1=%(kz)d l1=%(nz)d | add add=0.0 | put label1=z label2=x unit1=km unit2=km ''' % par) # Delta Flow('delo',None, ''' spike nsp=1 mag=0 n1=%(nz)d o1=%(oz)g d1=%(dz)g n2=%(nx)d o2=%(ox)g d2=%(dx)g k1=%(kz)d l1=%(nz)d | add add=0.0 | put label1=z label2=x unit1=km unit2=km ''' % par) # mask: selecting the perturbation region par['kzs']=0.53*par['nz'] + 1 par['lzs']=0.67/par['dx'] + 1 par['kxs']=0.7/par['dx'] + 1 par['lxs']=1.3/par['dx'] + 1 Flow('mask',None, ''' spike nsp=1 mag=1 n1=%(nz)d o1=%(oz)g d1=%(dz)g k1=%(kzs)d l1=%(lzs)d n2=%(nx)d o2=%(ox)g d2=%(dx)g k2=%(kxs)d l2=%(lxs)d | smooth rect1=20 rect2=30 repeat=3 | put label1=z label2=x unit1=km unit2=km ''' % par) Result('mask',fdmod.cgrey('title="mask" wantscalebar=y',par)) # Parameters for generating random perturbations par['ff']=9.8 par['aa']=1.0 #6.3 # distribution shape parameter 'alpha' par['ru']=0.036 par['rv']=0.009 # Use package gfield.py to generate perturbations gfield.execute('hh',122009,par) Result('hh',fdmod.cgrey('wantscalebar=y color=j',par)) Flow('gg','hh mask','add mode=p \${SOURCES[1]}') Result('gg',fdmod.cgrey('wantscalebar=y color=j',par)) Flow('kk','gg','add scale=0.75 | add add=1') Result('kk',fdmod.cgrey('wantscalebar=y bias=1.0 color=j',par)) #-------------------------------------------------------------- # Final velocity and density fields Flow('vp','vpbo kk','add mode=p \${SOURCES[1]}') Flow('vs','vsbo','cat') Flow('ro','rob kk','add mode=p \${SOURCES[1]}') Flow('vph','vpo kk','add mode=p \${SOURCES[1]}') Flow('vsh','vso','cat') Flow('roh','roo','cat') ### Set up ACQUISITION geometry # source/receiver positions #fdmod.point3('ss', # par['ox']+(par['nx']/2*par['dx']), # par['oz'],1,par) #fdmod.horizontal('rrt',par['dz'],par) #fdmod.horizontal('rrb',par['oz']+(par['nz']*par['dx']),par) #fdmod.vertical('rrl',par['dx'],par) #fdmod.vertical('rrr',par['ox']+(par['nx']*par['dx']),par) #Flow('rr','rrt rrr rrb rrl','cat \${SOURCES[1]} \${SOURCES[2]} \${SOURCES[3]}') #fdmod.boxarray('rr',par['nz'],par['oz'],par['dz'],par['nx'],par['ox'],par['dx'],par) # first receiver point par['sx'] = 0.9 par['sz'] = 0.15 fdmod.point3('ss', par['sx'], par['sz'],1,par) # second receiver point par['sx2'] = 1.0 par['sz2'] = 0.25 fdmod.point3('ss2', par['sx2'], par['sz2'],1,par) ## Integration/Source "surfaces" (they're lines in 2D!) # top par['rz1']=0.08 fdmod.horizontal('rr1',par['rz1'],par) # 5 different bottom integration/source lines par['rz2']=0.30 fdmod.horizontal('rr2',par['rz2'],par) par['rz3']=0.45 fdmod.horizontal('rr3',par['rz3'],par) par['rz4']=0.60 fdmod.horizontal('rr4',par['rz4'],par) par['rz5']=0.75 fdmod.horizontal('rr5',par['rz5'],par) par['rz6']=0.90 fdmod.horizontal('rr6',par['rz6'],par) Flow('rr','rr1 rr2 rr3 rr4 rr5 rr6','cat axis=2 space=n \${SOURCES[1]} \${SOURCES[2]} \${SOURCES[3]} \${SOURCES[4]} \${SOURCES[5]}') # ------------------------------------------------------------ # source/receiver coordinates Plot('rr1','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Plot('rr2','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Plot('rr3','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Plot('rr4','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Plot('rr5','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Plot('rr6','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Plot('rr','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Plot('ss','window n1=2 | dd type=complex | window | ' + fdmod.cgraph('wantscalebar=y symbol=o plotcol=1',par)) Plot('ss2','window n1=2 | dd type=complex | window | ' + fdmod.cgraph('wantscalebar=y symbol=o plotcol=1',par)) # ------------------------------------------------------------ ## Make wavelets/source functions # dipole source fdmod.wavelet('hor0',par['frq'],par) fdmod.wavelet('ver0',par['frq'],par) # Make wavelet for vertical dipole Flow('hor','hor0','add scale=0') Flow('ver','ver0','add scale=1') Flow('wave0','ver hor','cat axis=2 space=n \${SOURCES[1:2]}') Flow('wave','wave0', ''' transp plane=12 | transp plane=23 | transp plane=12 ''') Plot('ver','wave','window n2=1 f2=0 | window n1=500 |' + fdmod.waveplot('title="Elastic vertical source"',par)) Plot('hor','wave','window n2=1 f2=1 | window n1=500 |' + fdmod.waveplot('title="Elastic horizontal source"',par)) Plot('wave','hor ver','Movie',view=1) # Make wavelet for horizontal dipole Flow('horh','hor0','add scale=1') Flow('verh','ver0','add scale=0') Flow('wave0h','verh horh','cat axis=2 space=n \${SOURCES[1:2]}') Flow('waveh','wave0h', ''' transp plane=12 | transp plane=23 | transp plane=12 ''') Plot('verh','waveh','window n2=1 f2=0 | window n1=500 |' + fdmod.waveplot('title="Elastic vertical source"',par)) Plot('horh','waveh','window n2=1 f2=1 | window n1=500 |' + fdmod.waveplot('title="Elastic horizontal source"',par)) Plot('waveh','horh verh','Movie',view=1) # ------------------------------------------------------------ # plot final velocity/parameter models Plot('vp', fdmod.cgrey('wantscalebar=y allpos=y color=j bias=2.0 pclip=99.9 title="vp"',par)) Plot('vs', fdmod.cgrey('wantscalebar=y allpos=y color=j bias=0.0 pclip=96 title="vs"',par)) Plot('vph', fdmod.cgrey('wantscalebar=y allpos=y color=j bias=2.0 pclip=99.9 title="vph"',par)) Plot('vsh', fdmod.cgrey('wantscalebar=y allpos=y color=j bias=0.0 pclip=100 title="vsh"',par)) Plot('roh', fdmod.cgrey('wantscalebar=y allpos=y color=j bias=2000000 pclip=100 title="roh"',par)) Plot('ro', fdmod.cgrey('wantscalebar=y allpos=y color=j bias=1500000 pclip=100 title="density"',par)) Plot('epsilon',fdmod.cgrey('wantscalebar=y allpos=y pclip=100',par)) Plot('delta', fdmod.cgrey('wantscalebar=y allpos=y pclip=100',par)) Result('vp', ['vp', 'ss','rr'],'Overlay') Result('vs', ['vs', 'ss','rr'],'Overlay') Result('vph', ['vph', 'ss','rr'],'Overlay') Result('vsh', ['vsh', 'ss','rr'],'Overlay') Result('roh', ['roh', 'ss','rr'],'Overlay') Result('ro', ['ro', 'ss','rr'],'Overlay') Result('epsilon',['epsilon','ss','rr'],'Overlay') Result('delta', ['delta', 'ss','rr'],'Overlay') # Compute Cij's from velocities, densities & Thomsen parameters fdmod.anisotropic('cc','vph','vsh','roo','epso','delo',par) fdmod.anisotropic('cchomo','vpo','vso','roo','epso','delo',par) #-------------------------------------------------------------------------------------------- # vector-acoustic modeling - MAKING DATA #-------------------------------------------------------------------------------------------- ## # particle velocity field without free surface, par['ssou']='y' rules.emodel('de','we','wave','cc','roo','ss','rr','ssou=%(ssou)s opot=n' % par,par) Flow('def','de', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) for i in range(2): Flow('we'+str(i+1),'we', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('we'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('def'+str(i+1),'def', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Flow('weall','we1 we2','cat axis=1 space=n \${SOURCES[1]}') Plot('weall', ''' grey title="Vector components" wantaxis=y screenratio=%f screenht=8 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.25 g2num=0.25 ''' % (2*par['ratio']),view=1) Flow('wall','we1 we2','cat axis=1 space=n \${SOURCES[1]}') Plot('wall', ''' grey title="" wantaxis=y screenratio=%f screenht=10 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.1 g2num=0.1 ''' % (3*par['ratio']),view=1) # # # vertical point-force with free surface, no scattering in the subsurface. Used to compute "THE ANSWER". rules.emodel('deo','weo','wave','cchomo','roo','ss','rr','ssou=%(ssou)s opot=n' % par,par) Flow('deof','deo', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) # get scattered waves Flow('desf','def deof','add \${SOURCES[1]} scale=1,-1 out=stdout') Flow('wes','we weo','add \${SOURCES[1]} scale=1,-1') for i in range(2): Flow('weo'+str(i+1),'weo', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('weo'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('deof'+str(i+1),'deof', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Plot('desf'+str(i+1),'desf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) ## # pressure fields par['ssou']='y' rules.emodel('dep','wep','wave','cc','roo','ss','rr','ssou=%(ssou)s opot=y' % par,par) Flow('depf','dep', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) for i in range(2): Flow('wep'+str(i+1),'wep', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('wep'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('depf'+str(i+1),'depf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Flow('wepall','wep1 wep2','cat axis=1 space=n \${SOURCES[1]}') Plot('wepall', ''' grey title="Vector components" wantaxis=y screenratio=%f screenht=8 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.25 g2num=0.25 ''' % (2*par['ratio']),view=1) Flow('wpall','wep1 wep2','cat axis=1 space=n \${SOURCES[1]}') Plot('wpall', ''' grey title="" wantaxis=y screenratio=%f screenht=10 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.1 g2num=0.1 ''' % (3*par['ratio']),view=1) # # # vertical point-force with free surface, no scattering in the subsurface. Used to compute "THE ANSWER". rules.emodel('depo','wepo','wave','cchomo','roo','ss','rr','ssou=%(ssou)s opot=y' % par,par) Flow('depof','depo', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) # get scattered waves Flow('depsf','depf depof','add \${SOURCES[1]} scale=1,-1 out=stdout') Flow('weps','wep wepo','add \${SOURCES[1]} scale=1,-1') for i in range(2): Flow('wepo'+str(i+1),'wepo', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('wepo'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('depof'+str(i+1),'depof', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Plot('depsf'+str(i+1),'depsf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) ## SECOND RECEIVER # particle velocity field without free surface, par['ssou']='y' rules.emodel('der','wer','wave','cc','roo','ss2','rr','ssou=%(ssou)s opot=n' % par,par) Flow('derf','der', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) for i in range(2): Flow('wer'+str(i+1),'wer', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('wer'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('derf'+str(i+1),'derf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Flow('werall','wer1 wer2','cat axis=1 space=n \${SOURCES[1]}') Plot('werall', ''' grey title="Vector components" wantaxis=y screenratio=%f screenht=8 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.25 g2num=0.25 ''' % (2*par['ratio']),view=1) Flow('wrall','wer1 wer2','cat axis=1 space=n \${SOURCES[1]}') Plot('wrall', ''' grey title="" wantaxis=y screenratio=%f screenht=10 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.1 g2num=0.1 ''' % (3*par['ratio']),view=1) # # # vertical point-force with free surface, no scattering in the subsurface. Used to compute "THE ANSWER". rules.emodel('deor','weor','wave','cchomo','roo','ss2','rr','ssou=%(ssou)s opot=n' % par,par) Flow('deorf','deor', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) # get scattered waves Flow('desrf','derf deorf','add \${SOURCES[1]} scale=1,-1 out=stdout') Flow('wesr','wer weor','add \${SOURCES[1]} scale=1,-1') for i in range(2): Flow('weor'+str(i+1),'weor', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('weor'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('deorf'+str(i+1),'deorf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Plot('desrf'+str(i+1),'desrf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) ## # pressure fields par['ssou']='y' rules.emodel('depr','wepr','wave','cc','roo','ss2','rr','ssou=%(ssou)s opot=y' % par,par) Flow('deprf','depr', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) for i in range(2): Flow('wepr'+str(i+1),'wepr', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('wepr'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('deprf'+str(i+1),'deprf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Flow('weprall','wepr1 wepr2','cat axis=1 space=n \${SOURCES[1]}') Plot('weprall', ''' grey title="Vector components" wantaxis=y screenratio=%f screenht=8 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.25 g2num=0.25 ''' % (2*par['ratio']),view=1) Flow('wprall','wepr1 wepr2','cat axis=1 space=n \${SOURCES[1]}') Plot('wprall', ''' grey title="" wantaxis=y screenratio=%f screenht=10 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.1 g2num=0.1 ''' % (3*par['ratio']),view=1) # # # vertical point-force with free surface, no scattering in the subsurface. Used to compute "THE ANSWER". rules.emodel('depor','wepor','wave','cchomo','roo','ss2','rr','ssou=%(ssou)s opot=y' % par,par) Flow('deporf','depor', '''put n1=%(nx)d d1=%(dx)g o1=%(ox)d n2=6 d2=1 o2=0 n3=2 d3=1 o3=0 n4=%(nt)d d4=%(dt)g o4=%(ot)d out=stdout ''' % par ) # get scattered waves Flow('depsrf','deprf deporf','add \${SOURCES[1]} scale=1,-1 out=stdout') Flow('wepsr','wepr wepor','add \${SOURCES[1]} scale=1,-1') for i in range(2): Flow('wepor'+str(i+1),'wepor', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('wepor'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('deporf'+str(i+1),'deporf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Plot('depsrf'+str(i+1),'depsrf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) ## ## TRUE RESPONSE # pressure fields par['ssou']='y' rules.emodel('dept','wept','wave','cc','roo','ss2','ss','ssou=%(ssou)s opot=y' % par,par) Flow('deptf','dept', '''cat out=stdout ''') for i in range(2): Flow('wept'+str(i+1),'wept', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('wept'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) # Plot('deptf'+str(i+1),'deptf', # ''' # window n3=1 f3=%d squeeze=n | # transp plane=34 | transp plane=23 | transp plane=12 | # window f1=%d | put o1=%g | pad end1=%d | # ''' % (i,par['kt'],par['ot'],par['kt']) # + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Flow('weptall','wept1 wept2','cat axis=1 space=n \${SOURCES[1]}') Plot('weptall', ''' grey title="Vector components" wantaxis=y screenratio=%f screenht=8 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.25 g2num=0.25 ''' % (2*par['ratio']),view=1) Flow('wptall','wept1 wept2','cat axis=1 space=n \${SOURCES[1]}') Plot('wptall', ''' grey title="" wantaxis=y screenratio=%f screenht=10 gainpanel=a pclip=99 grid1=y grid2=y g1num=0.1 g2num=0.1 ''' % (3*par['ratio']),view=1) # # # vertical point-force with free surface, no scattering in the subsurface. Used to compute "THE ANSWER". rules.emodel('depot','wepot','wave','cchomo','roo','ss2','ss','ssou=%(ssou)s opot=y' % par,par) Flow('depotf','depot', '''cat out=stdout ''') # get scattered waves Flow('depstf','deptf depotf','add \${SOURCES[1]} scale=1,-1 out=stdout') Flow('wepst','wept wepot','add \${SOURCES[1]} scale=1,-1') for i in range(2): Flow('wepot'+str(i+1),'wepot', ''' window n3=1 f3=%d | window min1=%g max1=%g min2=%g max2=%g | scale axis=123 ''' % (i,par['zmin'],par['zmax'],par['xmin'],par['xmax'])) Plot('wepot'+str(i+1), fdmod.wgrey('title=u%s pclip=99' % str(i+1),par),view=1) Plot('depotf'+str(i+1),'depotf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) Plot('depstf'+str(i+1),'depstf', ''' window n3=1 f3=%d squeeze=n | transp plane=34 | transp plane=23 | transp plane=12 | window f1=%d | put o1=%g | pad end1=%d | ''' % (i,par['kt'],par['ot'],par['kt']) + fdmod.dgrey('title=u%s pclip=99 grid=y' %str(i+1),par),view=1) # Result('depstf', ''' window min2=0 max2=0 min3=0 max3=1.0 squeeze=n | transp plane=13 | window min1=0.098 max1=1.0 | pad end1=250 | put o1=0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="True" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) Plot('depstf', ''' window min2=0 max2=0 min3=0 max3=1.0 squeeze=n | transp plane=13 | window min1=0.098 max1=1.0 | pad end1=250 | put o1=0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="True" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) # Result('deptf', ''' window min2=0 max2=0 min3=0 max3=1.0 squeeze=n | transp plane=13 | window min1=0.098 max1=1.0 | pad end1=250 | put o1=0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="True" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) Plot('deptf', ''' window min2=0 max2=0 min3=0 max3=1.0 squeeze=n | transp plane=13 | window min1=0.098 max1=1.0 | pad end1=250 | put o1=0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="True" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) # Result('depotf', ''' window min2=0 max2=0 min3=0 max3=1.0 squeeze=n | transp plane=13 | window min1=0.098 max1=1.0 | pad end1=250 | put o1=0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="True" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) Plot('depotf', ''' window min2=0 max2=0 min3=0 max3=1.0 squeeze=n | transp plane=13 | window min1=0.098 max1=1.0 | pad end1=250 | put o1=0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="True" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) ##################################################################### #### Interferometry ##################################################################### ## Make frequency-domain data # QUESTION: Can you name these files (in comment lines) as G(A,B), G_S(A,B) and G_0(A,B)? for i in range(6): Flow('def_w_'+str(i+1),'def',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('depf_w_'+str(i+1),'depf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('derf_w_'+str(i+1),'derf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('deprf_w_'+str(i+1),'deprf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('depsrf_w_'+str(i+1),'depsrf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('deporf_w_'+str(i+1),'deporf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('desrf_w_'+str(i+1),'desrf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('deorf_w_'+str(i+1),'deorf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('depsf_w_'+str(i+1),'depsf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('depof_w_'+str(i+1),'depof',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('desf_w_'+str(i+1),'desf',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) Flow('deof_w_'+str(i+1),'deof',''' window min2=%d max2=%d min3=0 max3=0 squeeze=n | transp plane=41 | sfcausint | pad beg1=5001 | fft1 inv=n opt=y sym=n | transp plane=42 ''' % (i,i)) ## Compute frequency-domain integrands ## Analysis of TOP integral for FULL-FIELDS # Term: p_(r_A) v_^{*}(r_B) Flow('full_pAS_vB0_top',['deprf_w_1','def_w_1'],''' math y=\${SOURCES[0]} x=\${SOURCES[1]} output='y*(conj(x))' ''') Flow('full_pAS_vB0_topt','full_pAS_vB0_top',''' fft1 inv=y opt=y sym=y ''') # Term: v_(r_A) p_^{*}(r_B) Flow('full_vAS_pB0_top',['derf_w_1','depf_w_1'],''' math y=\${SOURCES[0]} x=\${SOURCES[1]} output='y*(conj(x))' ''') Flow('full_vAS_pB0_topt','full_vAS_pB0_top',''' fft1 inv=y opt=y sym=n ''') # #Result('vph_top', ['vph','ss','ss2','rr1'],'Overlay') #Plot('vph_top', ['vph','ss','ss2','rr1'],'Overlay') Result('full_pAS_vB0_topt', '''window j1=25 | grey title="top: p V*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_pAS_vB0_topt', '''window j1=25 | grey title="top: p V*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Result('full_vAS_pB0_topt', '''window j1=25 | grey title="top: V p*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_vAS_pB0_topt', '''window j1=25 | grey title="top: V p*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Result('full_terms_top','vph_top full_pAS_vB0_topt full_vAS_pB0_topt','SideBySideAniso') # total top integrand and integration: Flow('full_Integrand_top','full_pAS_vB0_topt full_vAS_pB0_topt','add \${SOURCES[1]} scale=1,1') Result('full_Integrand_top', '''window j1=10 | grey title="top: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_Integrand_top', ''' grey title="top: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) # taper for integration # QUESTION: why do we need it? par['ktxs']=0.25/par['dx'] + 1 par['ltxs']=1.75/par['dx'] + 1 Flow('taper',None, ''' spike nsp=1 mag=1 n1=10002 o1=-1.0002 d1=%(dt)g k1=1 l1=10002 n2=%(nx)d o2=%(ox)g d2=%(dx)g k2=%(ktxs)d l2=%(ltxs)d | smooth rect1=50 rect2=50 repeat=3 | put label1=z label2=x unit1=km unit2=km ''' % par) Result('taper','window j1=10 | ' + fdmod.cgrey('title="taper" wantscalebar=y',par)) # Apply taper to integrand Flow('full_Int_taper_top','full_Integrand_top taper','add \${SOURCES[1]} mode=p scale=1,-1') Result('full_Int_taper_top', '''window j1=25 | grey title="top: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_Int_taper_top', '''window j1=25 | grey title="top: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) # compute top contribution Flow('full_Contribution_top','full_Int_taper_top','stack axis=2') Result('full_Contribution_top', ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Top" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) Plot('full_Contribution_top', ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Top" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) Result('full_Cont_top','vph_top full_Int_taper_top full_Contribution_top','SideBySideAniso') ## Analysis of TOP integral for SCATTERING # Term: p_S(r_A) v_0^{*}(r_B) Flow('pAS_vB0_top',['depsrf_w_1','deof_w_1'],''' math y=\${SOURCES[0]} x=\${SOURCES[1]} output='y*(conj(x))' ''') Flow('pAS_vB0_topt','pAS_vB0_top',''' fft1 inv=y opt=y sym=y''') # Term: v_S(r_A) p_0^{*}(r_B) Flow('vAS_pB0_top',['desrf_w_1','depof_w_1'],''' math y=\${SOURCES[0]} x=\${SOURCES[1]} output='y*(conj(x))' ''') Flow('vAS_pB0_topt','vAS_pB0_top',''' fft1 inv=y opt=y sym=n ''') # Result('vph_top', ['vph','ss','ss2','rr1'],'Overlay') Plot('vph_top', ['vph','ss','ss2','rr1'],'Overlay') Result('pAS_vB0_topt', '''window j1=25 | grey title="top: p_S V_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('pAS_vB0_topt', '''window j1=25 | grey title="top: p_S V_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Result('vAS_pB0_topt', '''window j1=25 | grey title="top: v_S p_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('vAS_pB0_topt', '''window j1=25 | grey title="top: v_S p_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Result('terms_top','vph_top pAS_vB0_topt vAS_pB0_topt','SideBySideAniso') # total top integrand and integration: Flow('Integrand_top','pAS_vB0_topt vAS_pB0_topt','add \${SOURCES[1]} scale=1,1') Result('Integrand_top', '''window j1=10 | grey title="top: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('Integrand_top', ''' grey title="top: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) # taper Flow('Int_taper_top','Integrand_top taper','add \${SOURCES[1]} mode=p scale=1,1') Result('Int_taper_top', '''window j1=25 | grey title="top: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) Plot('Int_taper_top', '''window j1=25 | grey title="top: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (par['lt'],par['ut'],par['lx'],par['ux'])) # compute top contribution Flow('Contribution_top','Int_taper_top','stack axis=2 | causint') Result('Contribution_top', ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Top" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) Plot('Contribution_top', ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Top" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (par['lt'],par['ut'])) Result('Cont_top','vph_top Int_taper_top Contribution_top','SideBySideAniso') ### Analysis of bottom integral contributions # Plot bottom surfaces Flow('rrb','rr2 rr3 rr4 rr5 rr6','cat axis=2 space=n \${SOURCES[1]} \${SOURCES[2]} \${SOURCES[3]} \${SOURCES[4]}') Plot('rrb','window n1=2 | dd type=complex | window j2=10 | ' + fdmod.cgraph('wantscalebar=y symbol=x plotcol=1',par)) Result('vph_bottom', ['vph','ss','ss2','rrb'],'Overlay') Plot('vph_bottom', ['vph','ss','ss2','rrb'],'Overlay') ## SCATTERING terms... for i in range(5): # model and surface Result('vph_bot_'+str(i+2), ['vph','ss','ss2','rr'+str(i+2)],'Overlay') Plot('vph_bot_'+str(i+2), ['vph','ss','ss2','rr'+str(i+2)],'Overlay') # Term: p_S(r_A) v_0^{*}(r_B) Flow('pAS_vB0_bot_'+str(i+2),['depsrf_w_'+str(i+2),'deof_w_'+str(i+2)],''' math x=\${SOURCES[1]} output='input*(conj(x))' ''') Flow('pAS_vB0_bott_'+str(i+2),'pAS_vB0_bot_'+str(i+2),''' fft1 inv=y opt=y sym=y ''') # Term: v_S(r_A) p_0^{*}(r_B) Flow('vAS_pB0_bot_'+str(i+2),['desrf_w_'+str(i+2),'depof_w_'+str(i+2)],''' math x=\${SOURCES[1]} output='input*(conj(x))' ''') Flow('vAS_pB0_bott_'+str(i+2),'vAS_pB0_bot_'+str(i+2),''' fft1 inv=y opt=y sym=n ''') # Plot terms Result('pAS_vB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: p_S V_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('pAS_vB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: p_S V_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Result('vAS_pB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: v_S p_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('vAS_pB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: v_S p_0*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Result('terms_bot_'+str(i+2),['vph_bot_'+str(i+2),'pAS_vB0_bott_'+str(i+2),'vAS_pB0_bott_'+str(i+2)],'SideBySideAniso') # sum terms Flow('Integrand_bot_'+str(i+2),['pAS_vB0_bott_'+str(i+2),'vAS_pB0_bott_'+str(i+2)],'add \${SOURCES[1]} mode=a scale=1,-1') Result('Integrand_bot_'+str(i+2), '''window j1=10 | grey title="bottom #%d: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('Integrand_bot_'+str(i+2), ''' grey title="bottom #%d: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) # taper integrand Flow('Int_taper_bot_'+str(i+2),['Integrand_bot_'+str(i+2),'taper'],'add \${SOURCES[1]} mode=p scale=1,-1') Result('Int_taper_bot_'+str(i+2), '''window j1=25 | grey title="bottom #%d: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('Int_taper_bot_'+str(i+2), '''window j1=25 | grey title="bottom #%d: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) # get final contribution Flow('Contribution_bot_'+str(i+2),'Int_taper_bot_'+str(i+2),'stack axis=2') Result('Contribution_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Plot('Contribution_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Flow('Contribution_top_bot_'+str(i+2),['Contribution_top','Contribution_bot_'+str(i+2)],'cat \${SOURCES[1]} axis=2') Result('Contribution_top_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Plot('Contribution_top_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Result('Cont_bot_'+str(i+2),['vph_bot_'+str(i+2),'Int_taper_bot_'+str(i+2),'Contribution_top_bot_'+str(i+2)],'SideBySideAniso') # ## FULL-FIELD terms... for i in range(5): # Term: p(r_A) v^{*}(r_B) Flow('full_pAS_vB0_bot_'+str(i+2),['deprf_w_'+str(i+2),'def_w_'+str(i+2)],''' math x=\${SOURCES[1]} output='input*(conj(x))' ''') Flow('full_pAS_vB0_bott_'+str(i+2),'full_pAS_vB0_bot_'+str(i+2),''' fft1 inv=y opt=y sym=y ''') # Term: v_S(r_A) p_0^{*}(r_B) Flow('full_vAS_pB0_bot_'+str(i+2),['derf_w_'+str(i+2),'depf_w_'+str(i+2)],''' math x=\${SOURCES[1]} output='input*(conj(x))' ''') Flow('full_vAS_pB0_bott_'+str(i+2),'full_vAS_pB0_bot_'+str(i+2),''' fft1 inv=y opt=y sym=n ''') # Plot terms Result('full_pAS_vB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: p V*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_pAS_vB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: p V*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Result('full_vAS_pB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: V p*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_vAS_pB0_bott_'+str(i+2), '''window j1=25 | grey title="bottom #%d: V p*" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Result('full_terms_bot_'+str(i+2),['vph_bot_'+str(i+2),'full_pAS_vB0_bott_'+str(i+2),'full_vAS_pB0_bott_'+str(i+2)],'SideBySideAniso') # sum terms Flow('full_Integrand_bot_'+str(i+2),['full_pAS_vB0_bott_'+str(i+2),'full_vAS_pB0_bott_'+str(i+2)],'add \${SOURCES[1]} mode=a scale=1,1') Result('full_Integrand_bot_'+str(i+2), '''window j1=10 | grey title="bottom #%d: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_Integrand_bot_'+str(i+2), ''' grey title="bottom #%d: total integrand" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) # taper integrand Flow('full_Int_taper_bot_'+str(i+2),['full_Integrand_bot_'+str(i+2),'taper'],'add \${SOURCES[1]} mode=p scale=1,1') Result('full_Int_taper_bot_'+str(i+2), '''window j1=25 | grey title="bottom #%d: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) Plot('full_Int_taper_bot_'+str(i+2), '''window j1=25 | grey title="bottom #%d: total integrand (w/ taper)" pclip=99 grid=y labelrot=n wantaxis=y label1=%s unit1=%s label2=%s unit2=%s ''' % (i+2,par['lt'],par['ut'],par['lx'],par['ux'])) # get final contribution Flow('full_Contribution_bot_'+str(i+2),'full_Int_taper_bot_'+str(i+2),'stack axis=2') Result('full_Contribution_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Plot('full_Contribution_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Flow('full_Contribution_top_bot_'+str(i+2),['full_Contribution_top','full_Contribution_bot_'+str(i+2)],'cat \${SOURCES[1]} axis=2') Result('full_Contribution_top_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Plot('full_Contribution_top_bot_'+str(i+2), ''' window min1=0 max1=1.0 | wiggle pclip=100 poly=y screenht=10 screenratio=4 transp=y title="Bottom #%d" labelrot=n wantaxis=y label1=%s unit1=%s yreverse=y ''' % (i+2,par['lt'],par['ut'])) Result('full_Cont_bot_'+str(i+2),['vph_bot_'+str(i+2),'full_Int_taper_bot_'+str(i+2),'full_Contribution_top_bot_'+str(i+2)],'SideBySideAniso') # ------------------------------------------------------------ # ------------------------------------------------------------ #rules.test('vp','vs','ro','epsilon','delta','ss','rr',par) # ------------------------------------------------------------ # ------------------------------------------------------------ End()```