```from rsf.proj import * dx = 0.5 for case in range(3): letter = 'npc'[case] velocity = 'g%s-velocity' % letter Flow(velocity,None, ''' math o1=0.0 n1=751 d1=0.004 label1=z unit1=km o2=4.0 n2=501 d2=0.008 label2=x unit2=km label=Velocity unit=km/s output="%g+(%g)*exp(-%g*((x1-1)^2+(x2-6)^2))" ''' % ((3,2,2)[case],(-2,2,0)[case],1.0/(dx*dx))) Result(velocity, ''' grey allpos=y bias=1 wanttitle=n wantscalebar=y bartype=h ''') Plot(velocity, 'grey allpos=y bias=1 wanttitle=n') rays = 'rays%d' % case Flow(rays,velocity, ''' rays2 dt=0.01 nt=1000 nr=101 a0=130 amax=230 zshot=0 yshot=6 sym=n | reverse which=2 ''') Plot(rays, ''' window n1=131 n2=90 f2=5 | graph min2=4 max2=8 min1=0 max1=3 transp=y yreverse=y plotcol=%d title="Ray Tracing" wantaxis=n wheretitle=b wherexlabel=t ''' % (7,0,0)[case]) Plot('o'+rays,[velocity,rays],'Overlay') hwt = 'hwt%d' % case Flow(hwt,velocity,'hwt2d dt=0.01 nt=1001 ng=101 og=-50 dg=1 xsou=6 zsou=0') Plot(hwt, ''' window n2=131 n1=90 f1=5 | transp | graph min1=4 max1=8 min2=0 max2=3 yreverse=y plotcol=%d title="Wavefront Tracing" wantaxis=n wheretitle=b wherexlabel=t ''' % (7,0,0)[case]) Plot('o'+hwt,[velocity,hwt],'Overlay') comp = 'g%s-velrw' % letter Result(comp,['o'+rays,'o'+hwt],'SideBySideIso') dif = 'g%s-diff' % letter Flow(dif,[hwt,rays], ''' dd type=float | put n1=2 n2=101 n3=1001 | reverse which=1 | put n1=202 o1=-50 d1=1 n2=1001 n3=1 | dd type=complex | transp | add scale=1,-1 \${SOURCES[1]} | window n1=81 n2=90 f2=5 | math output="abs(input)" | real ''') Result(dif, ''' grey allpos=y scalebar=y label1=Time unit1=s label2="Take-off Angle" unit2="\^o\_" barlabel=Distance barunit=km wanttitle=n ''') End()```