up [pdf]
from rsf.proj import *

Fetch('marmvel.hh','marm')
Flow('vel','marmvel.hh',
    '''
    dd form=native |
    put label1=Z label2=X label=Velocity unit=m/s |
    window j1=2 j2=4 |
    smooth rect1=20 rect2=20 repeat=1 | put d1=25 d2=25
    ''')

Result('mvel','vel',
       '''
       grey color=j allpos=y bias=1500 scalebar=y title="Marmousi Model" 
       barlabel=Velocity barreverse=y barunit=m/s
       labelfat=6 titlefat=6 labelsz=9 titlesz=10
       unit1=m unit2=m
       ''')

Flow('grad1','vel','igrad square=n adj=n | math output="input/25" ')
Flow('grad2','vel','transp| igrad square=n adj=n | transp | math output="input/25" ')
Result('grad1','grey scalebar=y color=j title="Velocity gradient in Z direction" ')
Result('grad2','grey scalebar=y color=j title="Velocity gradient in X direction" ')

time=1.8

nt=121
dt=0.015

dtt=0.0005
factor=dt/dtt
ntt=(nt-1)*factor+1
ktt=(ntt-1)/10+1

Flow('source1',None,
     '''
     spike n1=%d d1=%g k1=%d |
     ricker1 frequency=15
     '''%(ntt,dtt,ktt))
Flow('real','source1','math "output=0"')
Flow('imag','source1','envelope hilb=y | halfint | halfint | math output="input/2" ')

Flow('csource1','real imag','cmplx ${SOURCES[1]}')
Flow('csource','csource1','window j1=%d'% factor)
Flow('source','source1','window j1=%d'% factor)
Result('source','graph  title="Source Wavelet" ')

Flow('fft','vel','rtoc | fft3 axis=1 pad=1 | fft3 axis=2 pad=1')

Flow('refl','vel','spike k1=150 k2=288 | smooth rect1=2 rect2=2 repeat=2')

### traveltime
#Flow('ray','vel','spray axis=3 n=1 d=5 o=0 | eikonal zshot=1220 yshot=1250 xshot=0')
#Plot('ray','contour nc=1 c0=0.78 wantaxis=n wanttitle=n dash=1 plotfat=5 plotcol=3')

### gradient
Flow('right left','vel grad1 grad2 fft',
     '''
     cisolr2grad seed=2010 dt=%g grad1=${SOURCES[1]} grad2=${SOURCES[2]} fft=${SOURCES[3]} left=${TARGETS[1]} npk=50 eps=1e-4
     ''' % dt)

Flow('wave','csource refl left right',
     '''
     cfftwave2omp ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y cmplx=n
     ''')

Plot('wave',
     '''
     window j3=25  |
     grey label2="Z" label1="X" title="Isotropic"
     yreverse=y gainpanel=all 
     ''')

Result('wavesnap','wave',
     '''
     window n3=1 min3=%g |
     grey title="With Gradient"
     yreverse=y gainpanel=all 
     label1=Z unit1=m label2=X unit2=m
     labelfat=6 titlefat=6 labelsz=9 titlesz=10 
     ''' % time)

#Result('comp','wavesnap ray','Overlay')

Plot('wavecross','wave',
     '''
     window n3=1 min3=%g n1=1 min1=8000 |
     graph label2="Amplitude" label1="X" title="Isotropic" min2=-0.0004 max2=0.0005
     ''' % time)

### without gradient
Flow('right1 left1','vel fft',
     '''
     cisolr2 seed=2010 dt=%g fft=${SOURCES[1]} left=${TARGETS[1]} npk=50 eps=1e-4 
     ''' % dt)

Flow('wave1','csource refl left1 right1',
     '''
     cfftwave2 ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y cmplx=n
     ''')

Plot('wave1',
     '''
     window j3=250  |
     grey label2="X" label1="Z" title="Isotropic"
     yreverse=y gainpanel=all 
     ''')

Result('wave1snap','wave1',
     '''
     window n3=1 min3=%g |
     grey title="Without Gradient"
     yreverse=y gainpanel=all 
     label1=Z unit1=m label2=X unit2=m
     labelfat=6 titlefat=6 labelsz=9 titlesz=10 
     ''' % time)

#Result('comp1','wave1snap ray','Overlay')

Plot('wavecross1','wave1',
     '''
     window n3=1 min3=%g n1=1 min1=8000 |
     graph label2="Amplitude" label1="X" title="Isotropic" min2=-0.0004 max2=0.0005
     ''' % time)

### reference
nt=1201
dt=0.0015

dtt=0.0005
factor=dt/dtt
ntt=(nt-1)*factor+1
ktt=(ntt-1)/10+1

Flow('source1-b',None,
     '''
     spike n1=%d d1=%g k1=%d |
     ricker1 frequency=15
     '''%(ntt,dtt,ktt))
Flow('real-b','source1-b','math "output=0"')
Flow('imag-b','source1-b','envelope hilb=y | halfint | halfint | math output="input/2" ')

Flow('csource1-b','real-b imag-b','cmplx ${SOURCES[1]}')
Flow('csource-b','csource1-b','window j1=%d'% factor)
Flow('source-b','source1-b','window j1=%d'% factor)
Result('source-b','graph  title="Source Wavelet" ')

Flow('right1-b left1-b','vel fft',
     '''
     cisolr2 seed=2010 dt=%g fft=${SOURCES[1]} left=${TARGETS[1]} npk=50 eps=1e-4 
     ''' % dt)

Flow('wave1-b','csource-b refl left1-b right1-b',
     '''
     cfftwave2 ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y cmplx=n
     ''')

Plot('wave1-b',
     '''
     window j3=25  |
     grey label2="X" label1="Z" title="Isotropic"
     yreverse=y gainpanel=all 
     ''')

Result('wave1snap-b','wave1-b',
     '''
     window n3=1 min3=%g |
     grey title="Reference"
     yreverse=y gainpanel=all 
     label1=Z unit1=m label2=X unit2=m
     labelfat=6 titlefat=6 labelsz=9 titlesz=10 
     ''' % time)

#Result('comp1-b','wave1snap-b ray','Overlay')

Plot('wavecross1-b','wave1-b',
     '''
     window n3=1 min3=%g n1=1 min1=8000 |
     graph label2="Amplitude" label1="X" title="Isotropic" min2=-0.0004 max2=0.0005
     ''' % time)

Flow('last','wave','window n3=1 min3=%g | scale axis=3 ' % time)
Flow('last1','wave1','window n3=1 min3=%g | scale axis=3 '% time)
Flow('last2','wave1-b','window n3=1 min3=%g | scale axis=3 '% time)
Flow('diff','last last2','math truth=last2.rsf output="abs(input-truth)" ')
Flow('diff1','last1 last2','math truth=last2.rsf output="abs(input-truth)" ')
Result('mdiff','diff','grey maxval=0.15 minval=0 clip=0.15 scalebar=y label1=Z unit1=m label2=X unit2=m title="Error (with gradient)" labelfat=6 titlefat=6 labelsz=9 titlesz=10 barlabel=Normalized\ Amplitude')
Result('mdiff1','diff1','grey maxval=0.15 minval=0 clip=0.15 scalebar=y label1=Z unit1=m label2=X unit2=m title="Error (without gradient)" labelfat=6 titlefat=6 labelsz=9 titlesz=10 barlabel=Normalized\ Amplitude ')

pos=8800

Flow('trace','wave','window n1=1 min1=%d min3=%g n3=1 min2=4000 max2=10000 ' % (pos,time))
Flow('trace1','wave1','window n1=1 min1=%d min3=%g n3=1 min2=4000 max2=10000 ' % (pos,time))
Flow('trace2','wave1-b','window n1=1 min1=%d min3=%g n3=1 min2=4000 max2=10000 ' % (pos,time))

Flow('comp','trace trace1 trace2','cat axis=2 ${SOURCES[1:3]} | scale axis=1')
Result('mcompp','trace trace2','cat axis=2 ${SOURCES[1]} | scale axis=1 | graph dash=0,1 plotcol=6,7 plotfat=3 label1=X unit1=m label2="Normalized Amplitude" unit2= wanttitle=y title="With gradient" screenwd=14.5 screenht=5.5 labelfat=4 labelsz=5.5')
Result('mcompp1','trace1 trace2','cat axis=2 ${SOURCES[1]} | scale axis=1 | graph dash=4,1 plotcol=5,7 plotfat=3 label1=X unit1=m label2="Normalized Amplitude" unit2= wanttitle=y title="Without gradient" screenwd=14.5 screenht=5.5 labelfat=4 labelsz=5.5')

#Result('compp1','comp','graph dash=0,4,1 plotcol=6,5,7 plotfat=2 label1=Distance unit1=m label2="Normalized Amplitude" unit2= wanttitle=n labelfat=4 labelsz=9')


End()

sfdd
sfput
sfwindow
sfsmooth
sfgrey
sfigrad
sfmath
sftransp
sfspike
sfricker1
sfenvelope
sfhalfint
sfcmplx
sfgraph
sfrtoc
sffft3
sfcisolr2grad
sfcfftwave2omp
sfcisolr2
sfcfftwave2
sfscale
sfcat

data/marm/marmvel.hh