```from rsf.proj import * nt=351 dt=0.004 Flow('source',None, ''' spike n1=%d d1=%g k1=101 | ricker1 frequency=15 '''%(nt,dt)) #Parameter specification Flow('velx',None, ''' spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 | math output="1500+30*(x2-1.5)*(x2-1.5)+30*(x3-1)*(x3-1)+40*(x1)*(x1)" | scale dscale=0.001 ''') # Result('velx', # ''' # byte clip=820 bias=860 bar=bar1.rsf barreverse=y | # grey3 transp=n pclip=100 screenration=1 color=j # label1="X" label2="Y" label3="Z" scalebar=y poly=y bias=1000 # frame1=200 frame2=550 frame3=600 # unit1= unit2= title="V1" # ''') Flow('vely',None, ''' spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 | math output="1500+40*(x2-1.5)*(x2-1.5)+40*(x3-1)*(x3-1)+60*(x1)*(x1)" | scale dscale=0.001 ''') # Result('vely', # ''' # byte clip=820 bias=860 bar=bar2.rsf barreverse=y | # grey3 transp=n pclip=100 screenration=1 color=j # frame1=200 frame2=550 frame3=600 # unit1= unit2= title="V2" # ''') Flow('velz',None, ''' spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 | math output="1500+35*(x2-1.5)*(x2-1.5)+40*(x3-1)*(x3-1)+50*(x1)*(x1)" | scale dscale=0.001 ''') #Result('velz', # ''' # byte clip=820 bias=860 bar=bar.rsf barreverse=y | # grey3 transp=n pclip=100 screenration=1 color=j # frame1=200 frame2=550 frame3=600 # label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" # unit1=km unit2=km unit3=km flat=y screenratio=0.8 # scalebar=y barlabel="V\_p" barunit="km/s" # wanttitle=n # ''') Flow('eta1',None, ''' spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 | math output=0.3 ''') Flow('eta2',None, ''' spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 | math output=0.1 ''') Flow('gamma',None, ''' spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 | math output=1.03 ''') # Conversion to cij Flow('c33','velz', ''' math output='input*input' ''') Flow('c44','vely', ''' math output='0.5*input*input' ''') Flow('c55','velx', ''' math output='0.5*input*input' ''') Flow('c66','velx', ''' math output='0.25*input*input' ''') Flow('c13','velx c33 c55', ''' math c33=\${SOURCES[1]} c55=\${SOURCES[2]} output='-c55+sqrt(c55*c55+(input*input*(c33-c55)-c33*c55))' ''') Flow('c23','vely c33 c44', ''' math c33=\${SOURCES[1]} c44=\${SOURCES[2]} output='-c44+sqrt(c44*c44+(input*input*(c33-c44)-c33*c44))' ''') Flow('c11','eta1 c33 c13 c55', ''' math c33=\${SOURCES[1]} c13=\${SOURCES[2]} c55=\${SOURCES[3]} output='(input+0.5)*(2*c13*(c13+2*c55)+2*c33*c55)/(c33-c55)' ''') Flow('c22','eta2 c33 c23 c44', ''' math c33=\${SOURCES[1]} c23=\${SOURCES[2]} c44=\${SOURCES[3]} output='(input+0.5)*(2*c23*(c23+2*c44)+2*c33*c44)/(c33-c44)' ''') Flow('c12','gamma c11 c66', ''' math c11=\${SOURCES[1]} c66=\${SOURCES[2]} output='sqrt((input*input-1)*(c11*(c11-c66))+(c11-c66)*(c11-c66))-c66' ''') # Calculation of q Flow('q12','c11 c13 c33 c55',''' math c13=\${SOURCES[1]} c33=\${SOURCES[2]} c55=\${SOURCES[3]} output='(c55*(input-c55)+(c55+c13)*(c55+c13))/(c33*(input-c55))' ''') Flow('q32','c11 c13 c33 c55',''' math c13=\${SOURCES[1]} c33=\${SOURCES[2]} c55=\${SOURCES[3]} output='(c55*(c33-c55)+(c55+c13)*(c55+c13))/(input*(c33-c55))' ''') Flow('q21','c22 c23 c33 c44',''' math c13=\${SOURCES[1]} c33=\${SOURCES[2]} c55=\${SOURCES[3]} output='(c55*(input-c55)+(c55+c13)*(c55+c13))/(c33*(input-c55))' ''') Flow('q31','c22 c23 c33 c44',''' math c13=\${SOURCES[1]} c33=\${SOURCES[2]} c55=\${SOURCES[3]} output='(c55*(c33-c55)+(c55+c13)*(c55+c13))/(input*(c33-c55))' ''') Flow('q23','c22 c12 c11 c66',''' math c13=\${SOURCES[1]} c33=\${SOURCES[2]} c55=\${SOURCES[3]} output='(c55*(input-c55)+(c55+c13)*(c55+c13))/(c33*(input-c55))' ''') Flow('q13','c22 c12 c11 c66',''' math c13=\${SOURCES[1]} c33=\${SOURCES[2]} c55=\${SOURCES[3]} output='(c55*(c33-c55)+(c55+c13)*(c55+c13))/(input*(c33-c55))' ''') Flow('refl',None, ''' spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=128 k2=128 k3=128 | smooth rect1=2 rect2=2 rect3=3 repeat=3 ''') Flow('seta1','c11','math output=45') Flow('seta2','c11','math output=45') #Flow('seta','seta1 seta2','cat \${SOURCES[1]} axis=4') Flow('vp','c33','''math output='sqrt(input)' ''') # Lowrank decomposition ##################################################### Flow('fftc','vp','fft1 pad=1 | fft3 axis=2 pad=1 | fft3 axis=3 pad=1') # Exact ###################################################################### Flow('right0 left0','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2', ''' ortholr3 seed=2015 npk=10 eps=0.00001 dt=%g approx=0 tilt=y mode=0 fft=\${SOURCES[1]} c12=\${SOURCES[2]} c13=\${SOURCES[3]} c22=\${SOURCES[4]} c23=\${SOURCES[5]} c33=\${SOURCES[6]} c44=\${SOURCES[7]} c55=\${SOURCES[8]} c66=\${SOURCES[9]} seta1=\${SOURCES[10]} seta2=\${SOURCES[11]} left=\${TARGETS[1]} right=\${TARGETS[0]} ''' % dt) Flow('wave0 snaps0','source refl left0 right0', ''' fftwave3p ref=\${SOURCES[1]} left=\${SOURCES[2]} right=\${SOURCES[3]} verb=y snap=50 snaps=\${TARGETS[1]} ''') Result('wave0', ''' window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all | grey3 frame1=92 frame2=92 frame3=92 screenratio=1 labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" ''') Plot('snaps0', ''' byte gainpanel=all | grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7 wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n ''') # Zone approximation ###################################################################### Flow('right1 left1','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2', ''' ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=1 relation=3,3,3 tilt=y mode=0 fft=\${SOURCES[1]} c12=\${SOURCES[2]} c13=\${SOURCES[3]} c22=\${SOURCES[4]} c23=\${SOURCES[5]} c33=\${SOURCES[6]} c44=\${SOURCES[7]} c55=\${SOURCES[8]} c66=\${SOURCES[9]} seta1=\${SOURCES[10]} seta2=\${SOURCES[11]} left=\${TARGETS[1]} right=\${TARGETS[0]} ''' % dt) Flow('wave1 snaps1','source refl left1 right1', ''' fftwave3p ref=\${SOURCES[1]} left=\${SOURCES[2]} right=\${SOURCES[3]} verb=y snap=50 snaps=\${TARGETS[1]} ''') Result('wave1', ''' window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all | grey3 frame1=92 frame2=92 frame3=92 screenratio=1 labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" ''') Plot('snaps1', ''' byte gainpanel=all | grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7 wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n ''') # Acoustic approximation ###################################################################### Flow('right2 left2','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2', ''' ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=2 tilt=y mode=0 fft=\${SOURCES[1]} c12=\${SOURCES[2]} c13=\${SOURCES[3]} c22=\${SOURCES[4]} c23=\${SOURCES[5]} c33=\${SOURCES[6]} c44=\${SOURCES[7]} c55=\${SOURCES[8]} c66=\${SOURCES[9]} seta1=\${SOURCES[10]} seta2=\${SOURCES[11]} left=\${TARGETS[1]} right=\${TARGETS[0]} ''' % dt) Flow('wave2 snaps2','source refl left2 right2', ''' fftwave3p ref=\${SOURCES[1]} left=\${SOURCES[2]} right=\${SOURCES[3]} verb=y snap=50 snaps=\${TARGETS[1]} ''') Result('wave2', ''' window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all | grey3 frame1=92 frame2=92 frame3=92 screenratio=1 labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" ''') Plot('snaps2', ''' byte gainpanel=all | grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7 wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n ''') # Weak Anisotropy approximation ###################################################################### Flow('right3 left3','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2', ''' ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=3 tilt=y mode=0 fft=\${SOURCES[1]} c12=\${SOURCES[2]} c13=\${SOURCES[3]} c22=\${SOURCES[4]} c23=\${SOURCES[5]} c33=\${SOURCES[6]} c44=\${SOURCES[7]} c55=\${SOURCES[8]} c66=\${SOURCES[9]} seta1=\${SOURCES[10]} seta2=\${SOURCES[11]} left=\${TARGETS[1]} right=\${TARGETS[0]} ''' % dt) Flow('wave3 snaps3','source refl left3 right3', ''' fftwave3p ref=\${SOURCES[1]} left=\${SOURCES[2]} right=\${SOURCES[3]} verb=y snap=50 snaps=\${TARGETS[1]} ''') Result('wave3', ''' window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 | byte gainpanel=all | grey3 frame1=92 frame2=92 frame3=92 screenratio=1 labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" ''') Plot('snaps3', ''' byte gainpanel=all | grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7 wanttitle=n flat=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" ''') # Zone approximation (nine) ###################################################################### Flow('right4 left4','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2', ''' ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=1 nine=y tilt=y mode=0 fft=\${SOURCES[1]} c12=\${SOURCES[2]} c13=\${SOURCES[3]} c22=\${SOURCES[4]} c23=\${SOURCES[5]} c33=\${SOURCES[6]} c44=\${SOURCES[7]} c55=\${SOURCES[8]} c66=\${SOURCES[9]} seta1=\${SOURCES[10]} seta2=\${SOURCES[11]} left=\${TARGETS[1]} right=\${TARGETS[0]} ''' % dt) Flow('wave4 snaps4','source refl left4 right4', ''' fftwave3p ref=\${SOURCES[1]} left=\${SOURCES[2]} right=\${SOURCES[3]} verb=y snap=50 snaps=\${TARGETS[1]} ''') Result('wave4', ''' window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all | grey3 frame1=92 frame2=92 frame3=92 screenratio=1 labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" ''') Plot('snaps4', ''' byte gainpanel=all | grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7 wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n ''') clip2=0.0003 clip1=0.00015 for case in range(4): clip=clip1 if case+1 == 1: method = 'proposed approximation' elif case+1 == 2: method = 'acoustic approximation' elif case+1 == 3: method = 'Weak-anisotropy approximation' clip=clip2 else: method = 'proposed approximation (nine)' Flow('error%d' % (case+1),'wave0 wave%d' % (case+1),'add scale=-1,1 \${SOURCES[1]}| window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3') Result('error%d' % (case+1), ''' byte gainpanel=all bar=bar.rsf clip=%g allpos=y maxval=%g minval=0.0| grey3 label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" unit1=km unit2=km unit3=km screenratio=0.8 scalebar=y barlabel="Absolute Error" title="Amplitude error of the %s" frame1=92 frame2=92 frame3=92 flat=y wanttitle=n labelfat=6 titlefat=6 labelsz=9 titlesz=9 ''' % (clip,clip,method) ) End()```