from rsf.proj import *

# Make a reflector model 
Flow('refl',None,
     '''
     math n1=10001 o1=-4 d1=0.001 output="sqrt(1+x1*x1)" 
     ''')
Flow('model','refl',
     '''
     window min1=-3 max1=3 j1=5  |
     unif2 n1=401 d1=0.005 v00=1,2 |
     put label1=Depth unit1=km label2=Lateral unit2=km |
     smooth rect1=3
     ''')

Plot('model',
     '''
     window min2=-1 max2=2 |
     grey allpos=y bias=1 title=Model 
     ''')
Plot('refl',
     '''
     graph wanttitle=n wantaxis=n yreverse=y
     min1=-1 max1=2 min2=0 max2=2 
     plotfat=3 plotcol=4 
     ''')
Result('model','model refl','Overlay')

# Kirchoff zero-offset modeling
Flow('dip','refl','math output="x1/input" ')

Flow('data','refl dip',
     '''
     kirmod nt=1501 dt=0.004 freq=10
     ns=1201 s0=-3 ds=0.005 nh=1 dh=0.01 h0=0
     twod=y vel=1 dip=${SOURCES[1]} verb=y |
     window | costaper nw2=200 
     ''')
Result('data',
     '''
     window min2=-2 max2=2 min1=1 max1=4 |
     grey title="Zero-Offset Data" grid=y 
     label1=Time unit1=s label2=Distance unit2=km
     ''')

# Reverse-time migration
proj = Project()
prog = proj.Program('rtm.c')

Flow('image wave','model %s data' % prog[0],
     '''
     pad beg1=200 end1=200 | math output=0.5 |
     ./${SOURCES[1]} data=${SOURCES[2]} 
     wave=${TARGETS[1]} jt=20 n0=200
     ''')

# Display wave propagation movie
Plot('wave',
     '''
     window min2=-1 max2=2 min1=0 max1=2 f3=33 |
     grey wanttitle=n gainpanel=all
     ''',view=1)

# Display image
Plot('image',
     '''
     window min2=-1 max2=2 min1=0 max1=2 |
     grey title=Image 
     ''')
Result('image','image refl','Overlay')

End()