from rsf.proj import *

Fetch('marmvel.hh','marm')

Flow('raw','marmvel.hh',
     '''
     dd form=native | window j1=2 j2=2 f2=1000 n2=501 |
     scale dscale=0.001 |
     put o2=4 d2=0.008 d1=0.008
     label1=Depth unit1=km label2=Distance unit2=km 
     ''')
Flow('smooth','raw','smooth rect1=5 rect2=5 repeat=2')

for case in ('raw','smooth'):
    Plot(case,
         '''
         grey title="Marmousi Model: %s Velocity" allpos=y bias=1
         ''' % ('True','Smoothed')[case=='smooth'])
    Plot('b-'+case,case,'grey wanttitle=n allpos=y bias=1')

    # Ray tracing
    rays = 'rays-'+case
    Flow(rays,case,
         '''
         rays2 dt=0.01 nt=130 nr=201 a0=130 amax=230 zshot=0 yshot=6 sym=n |
         reverse which=2
         ''')
    Plot(rays,
         '''
         window f2=5 |
         graph min2=4 max2=8 min1=0 max1=3
         transp=y yreverse=y plotcol=0
         title="Ray Tracing" wantaxis=n wheretitle=b wherexlabel=t
         ''')

    # Hygens wavefront tracing
    huyg = 'huyg-'+case
    Flow(huyg,case,'hwt2d dt=0.01 nt=131 ng=201 og=-50 dg=0.5 xsou=6 zsou=0')
    Plot(huyg,
         '''
         window f1=5 |
         transp |
         graph min1=4 max1=8 min2=0 max2=3 yreverse=y plotcol=0
         title="Wavefront Tracing" wantaxis=n wheretitle=b wherexlabel=t
         ''')
    for meth in (rays,huyg):
        Plot('m-'+meth,['b-'+case,meth],'Overlay')

    Result('m-velrw-'+case,['m-rays-'+case,'m-huyg-'+case],'SideBySideIso')

Result('m-velocity','raw smooth','SideBySideIso')

Flow('diff','huyg-smooth rays-smooth',
     '''
     dd type=float | put n1=2 n2=201 n3=131 |
     reverse which=1 | put n1=402 o1=-50 d1=0.5 n2=131 n3=1 |
     dd type=complex | transp |
     add scale=1,-1 ${SOURCES[1]} |
     window n1=81 f2=5 |
     math output="abs(input)" | real
     ''')
Result('m-diff-smooth','diff',
       '''
       grey allpos=y scalebar=y 
       label1=Time unit1=s 
       label2="Take-off Angle" unit2="\^o\_"
       barlabel=Distance barunit=km
       wanttitle=n
       ''')

End()