from rsf.proj import *

#adapt the command-line version to be SConstruct.

#create 2D Gaussian function
#sfmath output="exp(-(x1*x1+x2*x2)/(2*1.5*1.5))" n1=200 \
#d1=0.1 o1=-10. n2=200 d2=0.1 o2=-10. | sfput label1=z \
#unit1=km label2=x unit2=km > gg.rsf

Flow('gg', None,
        '''
        math output="exp(-(x1*x1+x2*x2)/(2*1.5*1.5))" n1=200 
        d1=0.1 o1=-10. n2=200 d2=0.1 o2=-10. | 
        put label1=z unit1=km label2=x unit2=km
        ''')

#plot the 2D Gaussian function 
#< gg.rsf sfgrey pclip=100 screenratio=1 scalebar=y | xtpen
Result('gg','grey pclip=100 screenratio=1 scalebar=y wanttitle=n')


#############################################################
# Note that single quote only applies to one-line command. 
# Triple-quotes means multi-line documenting. sfprog will 
# omit the prefix 'sf' in Flow(). Madagascar adds the prefix 
# 'sf' automatically!
#############################################################


#extract 1D subset from the 2D Gaussian function
#< gg.rsf sfwindow n2=1 f2=100 | sfgraph | xtpen
Plot('gg','window n2=1 f2=100 | sfgraph')

#create a velocity model
#< gg.rsf sfscale rscale=-1. | sfadd add=3 > vel.rsf
Flow('vel','gg','sfscale rscale=-1. | sfadd add=3')

#< vel.rsf sfgrey title="" pclip=100 screenratio=1 \
#bias=3 scalebar=y|xtpen
Result('vel',
        ''' 
        grey title= pclip=100 screenratio=1 bias=3 scalebar=y
        ''')

#compute traveltimes with an eikonal solver
#< vel.rsf sfeikonal zshot=-10 yshot=0 > fme.rsf
#< fme.rsf sfcontour title="" nc=200 screenratio=1 |xtpen
Flow('fme','vel','eikonal zshot=-10 yshot=0')
Result('fme','contour title= nc=200 screenratio=1')

#compute rays and wavefronts
#< vel.rsf sfhwt2d xsou=0 zsou=-10 nt=1000 ot=0 dt=0.01 \
#ng=1801 og=-90 dg=0.1 > hwt.rsf
#< hwt.rsf sfwindow j1=20 j2=20 | sfgraph title="" \
#yreverse=y screenratio=1 min1=-10 max1=+10 min2=-10 \
#max2=+10 |xtpen
Flow('hwt','vel',
        '''
        sfhwt2d xsou=0 zsou=-10 nt=1000 ot=0 dt=0.01 ng=1801 
        og=-90 dg=0.1
        ''')
Result('hwt',
        '''
        sfwindow j1=20 j2=20 | sfgraph title="" yreverse=y 
        screenratio=1 min1=-10 max1=+10 min2=-10 max2=+10
        ''')

#############################################################
# Plot() and Result() have similar function for visualization 
# of your result: Here, you may use sfgrey, sfgraph, sfgrey3, 
# sfcontour, sfwiggle etc. The outcome of Result() will be 
# used as a final figure appeared in your paper. Two figures
# generated by Plot() may be overlaid in Result().
#############################################################

Plot('vel','grey title= pclip=100 screenratio=1 bias=3')
Plot('hwt',
        '''
        sfwindow j1=20 j2=20 | sfgraph title="" yreverse=y 
        screenratio=1 min1=-10 max1=+10 min2=-10 max2=+10
        ''')
Result('combined','vel hwt','Overlay')

End()