from rsf.proj import *

raw=['marmvel.hh']

# fetch Marmousi
Fetch('marmvel.hh',"marm")

Flow('marmvel','marmvel.hh',
     '''
     dd form=native | scale rscale=.001 | 
     put label1=Depth label2=Position unit1=km unit2=km d1=0.004 d2=0.004
     ''')

# smooth model
Flow('marm1','marmvel','window n1=7')
Flow('marm2','marmvel','window f1=7 n1=744 | smooth rect1=5 rect2=5 repeat=2')
Flow('marm','marm1 marm2','cat axis=1 ${SOURCES[1]}')

Plot('marm',
     '''
     window $SOURCE | grey color=j allpos=y scalebar=y 
     title="Smoothed Marmousi" barlabel=Velocity barunit=km\/s
     labelsz=10 titlesz=12 barreverse=y titlefat=6 labelfat=6
     ''')

# interpolation
###############

zpos = 500
xpos = 825
slice = 26

# shot0
Flow('ys0',None,'math n1=2 o1=3 d1=0.2 output=x1')
Flow('ss0','ys0','math output=0.')
Flow('s0','ss0 ys0','cat axis=2 ${SOURCES[1]} ${SOURCES[0]} | transp')

Flow('st0 stdl0 stds0','marm s0',
     '''
     put n3=1 d3=0.004 o3=0 label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} |
     put o4=3 d4=0.2
     ''')

# interpolation
for what in ('h','l','p'):
    Flow(what+'interp','st0 stds0',
         'tinterp deriv=${SOURCES[1]} type=%s ns=51 ds=0.004 os=3' % what)

    if what=='h':
        color=4
        dash=2
    elif what=='l':
        color=2
        dash=3
    else:
        color=7
        dash=6

    Plot(what+'c',what+'interp',
         '''
         window n1=1 f1=%d n2=1 f2=%d |
         graph plotfat=5 dash=%d plotcol=%d labelsz=7 titlesz=9
         wanttitle=n wantaxis=n
         ''' % (zpos,xpos,dash,color))

Plot('etime','st1',
     '''
     window n4=1 f4=%d | 
     grey color=j scalebar=y allpos=y barreverse=y minval=0 pclip=100
     labelsz=10 titlesz=12 title=First-arrivals barlabel=Traveltime barunit=s
     titlefat=6 labelfat=6
     ''' % slice)

Plot('hslice','hinterp st1',
     '''
     add scale=1,-1 ${SOURCES[1]} | window n4=1 f4=%d |
     grey color=j scalebar=y labelsz=10 titlesz=12 title="Error Hermite"
     minval=-0.025 maxval=0.025 clip=0.025 barlabel=Traveltime barunit=s
     titlefat=6 labelfat=6
     ''' % slice)
Plot('lslice','linterp st1',
     '''
     add scale=1,-1 ${SOURCES[1]} | window n4=1 f4=%d |
     grey color=j scalebar=y labelsz=10 titlesz=12 title="Error linear"
     minval=-0.025 maxval=0.025 clip=0.025 barlabel=Traveltime barunit=s
     titlefat=6 labelfat=6
     ''' % slice)
Plot('pslice','pinterp st1',
     '''
     add scale=1,-1 ${SOURCES[1]} | window n4=1 f4=%d |
     grey color=j scalebar=y labelsz=10 titlesz=12 title="Error shift"
     minval=-0.025 maxval=0.025 clip=0.025 barlabel=Traveltime barunit=s
     titlefat=6 labelfat=6
     ''' % slice)

Result('vel','marm hslice','OverUnderAniso')
Result('slice','lslice pslice','OverUnderAniso')

# compare
Flow('yc1',None,'math n1=51 o1=3 d1=0.004 output=x1')
Flow('sc1','yc1','math output=0.')
Flow('sc','sc1 yc1','cat axis=2 ${SOURCES[1]} ${SOURCES[0]} | transp')

Flow('st1','marm sc',
     '''
     put n3=1 d3=0.004 o3=0. label3= unit3= |
     eikonal shotfile=${SOURCES[1]} |
     put o4=3 d4=0.004
     ''')

Plot('c0','st1',
     '''
     window n1=1 f1=%d n2=1 f2=%d |
     graph plotfat=5 dash=0 plotcol=6 labelsz=6 titlesz=8
     label1="Source Position" unit1=km label2=Traveltime unit2=s
     title="Traveltime Interpolation" titlefat=4 labelfat=4
     ''' % (zpos,xpos))

Result('curve','c0 hc lc pc','Overlay')

# semi-recursive migration
##########################

dz=0.004
dx=0.004
xmin=3.
xmax=8.

# smooth model
Flow('marm1','marmvel','window n1=7')
Flow('marm2','marmvel','window f1=7 n1=744 | smooth rect1=5 rect2=5 repeat=2')
Flow('marm','marm1 marm2','cat axis=1 ${SOURCES[1]}')

# fetch data
Fetch('marmrefl.hh','marm')

Plot('zoom','marmvel',
     '''
     window n1=251 f1=500 n2=751 f2=1250 |
     grey title=Velocity screenratio=0.7
     labelsz=10 titlesz=12 mean=y
     titlefat=6 labelfat=6
     ''')

Flow('data','marmrefl.hh',
     '''
     dd form=native | 
     put label1=Time unit1=s 
     d2=0.025 o2=-2.575 label2=Offset unit2=km 
     d3=0.025 o3=3      label3=Shot   unit3=km
     ''')

# source
ns=31
ds=0.2
nr=43
dr=0.2

Flow('ys',None,'math n1=%d o1=3 d1=%g output=x1' % (ns,ds))
Flow('zs','ys','math output=0.008')
Flow('xs','ys','math output=0')
Flow('sht','zs ys xs','cat axis=2 ${SOURCES[1]} ${SOURCES[2]} | transp')

# receiver
Flow('yr',None,'math n1=%d o1=0.425 d1=%g output=x1' % (nr,dr))
Flow('zr','yr','math output=0.012')
Flow('xr','yr','math output=0')
Flow('rcv','zr yr xr','cat axis=2 ${SOURCES[1]} ${SOURCES[2]} | transp')

# mute refraction
Flow('mute','data','mutter half=n v0=1.6 tp=0.1')

##########
# standard
##########

# traveltime
Flow('time0s tdl0s tds0s','marm sht',
     '''
     put d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} b1=2 b2=2 |
     put o4=3 d4=%g | window
     ''' % (dx,ds))

Flow('time0r tdl0r tds0r','marm rcv',
     '''
     put d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} b1=2 b2=2 |
     put o4=0.425 d4=%g | window
     ''' % (dx,dr))

# migration
Flow('dmig0','mute time0s tds0s time0r tds0r',
     '''
     kirmigsr aperture=5 antialias=0.25 tau=0.0565 
     stable=${SOURCES[1]} sderiv=${SOURCES[2]} 
     rtable=${SOURCES[3]} rderiv=${SOURCES[4]}
     ''')

# plot
Result('dmig0',
       '''
       window f2=%d n2=%d |
       grey title=Standard screenht=7 screenwd=13
       labelsz=6 titlesz=8 label2=Position pclip=95
       titlefat=4 labelfat=4
       ''' % (xmin/dx,(xmax-xmin)/dx+1))

Plot('zoom0','dmig0',
     '''
     window n1=251 f1=500 n2=751 f2=1250 |
     grey title=Standard screenratio=0.7
     labelsz=10 titlesz=12 label2=Position pclip=90
     titlefat=6 labelfat=6
     ''')

# dense eikonal
Flow('ysd',None,'math n1=241 o1=3 d1=0.025 output=x1')
Flow('zsd','ysd','math output=0.008')
Flow('xsd','ysd','math output=0')
Flow('shtd','zsd ysd xsd','cat axis=2 ${SOURCES[1]} ${SOURCES[2]} | transp')

Flow('yrd',None,'math n1=337 o1=0.425 d1=0.025 output=x1')
Flow('zrd','yrd','math output=0.012')
Flow('xrd','yrd','math output=0')
Flow('rcvd','zrd yrd xrd','cat axis=2 ${SOURCES[1]} ${SOURCES[2]} | transp')

Flow('time0sd tdl0sd tds0sd','marm shtd',
     '''
     put d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} b1=2 b2=2 |
     put o4=3 d4=0.025 | window
     ''' % dx)

Flow('time0rd tdl0rd tds0rd','marm rcvd',
     '''
     put d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} b1=2 b2=2 |
     put o4=0.425 d4=0.025 | window
     ''' % dx)

Flow('dmig0d','mute time0sd tds0sd time0rd tds0rd',
     '''
     kirmigsr aperture=5 antialias=0.25 tau=0.0565 
     stable=${SOURCES[1]} sderiv=${SOURCES[2]} 
     rtable=${SOURCES[3]} rderiv=${SOURCES[4]}
     ''')

Result('dmig0d',
       '''
       window f2=%d n2=%d |
       grey title="Standard (without interpolation)" screenht=7 screenwd=13
       labelsz=6 titlesz=8 label2=Position pclip=95
       titlefat=4 labelfat=4
       ''' % (xmin/dx,(xmax-xmin)/dx+1))

#############
# recursive 1
#############

datum1=1.5
chunk1=2.

# chunk1
########

# traveltime
Flow('time1s_1','time0s','window n1=%d' % (chunk1/dz+1))
Flow('tds1s_1','tds0s','window n1=%d' % (chunk1/dz+1))

Flow('time1r_1','time0r','window n1=%d' % (chunk1/dz+1))
Flow('tds1r_1','tds0r','window n1=%d' % (chunk1/dz+1))

# migration
Flow('dmig1_1','mute time1s_1 tds1s_1 time1r_1 tds1r_1',
     '''
     kirmigsr aperture=5 antialias=0.25 tau=0.0565 
     stable=${SOURCES[1]} sderiv=${SOURCES[2]} 
     rtable=${SOURCES[3]} rderiv=${SOURCES[4]}
     ''')

# redatum
#########

Flow('sgreen1','time1s_1 tds1s_1',
     '''
     transp plane=34 memsize=4096 |
     tinterp deriv=${SOURCES[1]} os=3 ds=0.025 ns=240 |
     window n1=1 f1=%d |
     spline o1=3 d1=0.025 n1=240
     ''' % (datum1/dz))

Flow('rgreen1','time1r_1 tds1r_1',
     '''
     transp plane=34 memsize=4096 |
     tinterp deriv=${SOURCES[1]} os=0.425 ds=0.025 ns=335 |
     window n1=1 f1=%d |
     spline o1=0.425 d1=0.025 n1=335
     ''' % (datum1/dz))

Flow('redat1','mute sgreen1 rgreen1',
     '''
     reverse which=1 opt=i |
     kirdatsr sdatum=%g rdatum=%g aperture=300 taper=0 length=0.05
     sgreen=${SOURCES[1]} rgreen=${SOURCES[2]} |
     reverse which=1 opt=i
     ''' % (datum1-0.008,datum1-0.012))

# chunk2
########

Flow('chunk2','marm','window n1=%d f1=%d | put o1=0.'
     % (datum1/dz+1,(3.-datum1)/dz))

# source
Flow('ys1',None,'math n1=%d o1=3 d1=%g output=x1' % (ns,ds))
Flow('zs1','ys1','math output=0')
Flow('sht1','zs1 ys1','cat axis=2 ${SOURCES[1]} ${SOURCES[0]} | transp')

# receiver
Flow('yr1',None,'math n1=%d o1=0.425 d1=%g output=x1' % (nr,dr))
Flow('zr1','yr1','math output=0')
Flow('rcv1','zr1 yr1','cat axis=2 ${SOURCES[1]} ${SOURCES[0]} | transp')

# traveltime
Flow('time1s_2 tdl1s_2 tds1s_2','chunk2 sht1',
     '''
     put d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} |
     put o4=3 d4=%g | window
     ''' % (dx,ds))

Flow('time1r_2 tdl1r_2 tds1r_2','chunk2 rcv1',
     '''
     put d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} |
     put o4=0.425 d4=%g | window
     ''' % (dx,dr))

# migration
Flow('dmig1_2','redat1 time1s_2 tds1s_2 time1r_2 tds1r_2',
     '''
     kirmigsr aperture=2 antialias=0.25 tau=0.0565 
     stable=${SOURCES[1]} sderiv=${SOURCES[2]}
     rtable=${SOURCES[3]} rderiv=${SOURCES[4]}
     ''')

# plot
######

scale1=2.25

Flow('cut1','dmig1_2','window n1=%d f1=%d | scale dscale=%g'
     % ((3.-chunk1)/dz+1,(chunk1-datum1)/dz,scale1))
Flow('dmig1','dmig1_1 cut1','cat axis=1 ${SOURCES[1]}')

Result('dmig1',
       '''
       window f2=%d n2=%d |
       grey title="Recursive 1 step" screenht=7 screenwd=13
       labelsz=6 titlesz=8 label2=Position pclip=96
       titlefat=4 labelfat=4
       ''' % (xmin/dx,(xmax-xmin)/dx+1))

Plot('zoom1','dmig1',
     '''
     window n1=251 f1=500 n2=751 f2=1250 |
     grey title="Redatum 1 step" screenratio=0.7 
     labelsz=10 titlesz=12 label2=Position pclip=90
     titlefat=6 labelfat=6
     ''')

#############
# recursive 2
#############

datum2=0.5

# first
Flow('sgreen2_1','time1s_1 tds1s_1',
     '''
     transp plane=34 memsize=4096 |
     tinterp deriv=${SOURCES[1]} os=3 ds=0.025 ns=240 |
     window n1=1 f1=%d |
     spline o1=3 d1=0.025 n1=240
     ''' % (datum2/dz))

Flow('rgreen2_1','time1r_1 tds1r_1',
     '''
     transp plane=34 memsize=4096 |
     tinterp deriv=${SOURCES[1]} os=0.425 ds=0.025 ns=335 |
     window n1=1 f1=%d |
     spline o1=0.425 d1=0.025 n1=335
     ''' % (datum2/dz))

Flow('redat2_1','mute sgreen2_1 rgreen2_1',
     '''
     reverse which=1 opt=i |
     kirdatsr sdatum=%g rdatum=%g aperture=300 taper=0 length=0.05
     sgreen=${SOURCES[1]} rgreen=${SOURCES[2]} |
     reverse which=1 opt=i
     ''' % (datum2-0.008,datum2-0.012))

# second
Flow('time2s_2 tdl2s_2 tds2s_2','marm sht1',
     '''
     window n1=%d f1=%d | put o1=0. d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} |
     put o4=3 d4=%g
     ''' % (datum2/dz+1,datum2/dz,dx,ds))

Flow('time2r_2 tdl2r_2 tds2r_2','marm rcv1',
     '''
     window n1=%d f1=%d | put o1=0. d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} |
     put o4=0.425 d4=%g
     ''' % (datum2/dz+1,datum2/dz,dx,dr))

Flow('sgreen2_2','time2s_2 tds2s_2',
     '''
     tinterp deriv=${SOURCES[1]} os=3 ds=0.025 ns=240 |
     window n1=1 f1=%d |
     spline o1=3 d1=0.025 n1=240
     ''' % (datum2/dz))

Flow('rgreen2_2','time2r_2 tds2r_2',
     '''
     tinterp deriv=${SOURCES[1]} os=0.425 ds=0.025 ns=335 |
     window n1=1 f1=%d |
     spline o1=0.425 d1=0.025 n1=335
     ''' % (datum2/dz))

Flow('redat2_2','redat2_1 sgreen2_2 rgreen2_2',
     '''
     reverse which=1 opt=i |
     kirdatsr sdatum=%g rdatum=%g aperture=300 taper=0 length=0.05
     sgreen=${SOURCES[1]} rgreen=${SOURCES[2]} |
     reverse which=1 opt=i
     ''' % (datum2,datum2))

# third
Flow('time2s_3 tdl2s_3 tds2s_3','marm sht1',
     '''
     window n1=%d f1=%d | put o1=0. d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} |
     put o4=3 d4=%g
     ''' % (datum2/dz+1,2*datum2/dz,dx,ds))

Flow('time2r_3 tdl2r_3 tds2r_3','marm rcv1',
     '''
     window n1=%d f1=%d | put o1=0. d3=%g o3=0. label3= unit3= |
     eikods shotfile=${SOURCES[1]} tdl1=${TARGETS[1]} tds1=${TARGETS[2]} |
     put o4=0.425 d4=%g
     ''' % (datum2/dz+1,2*datum2/dz,dx,dr))

Flow('sgreen2_3','time2s_3 tds2s_3',
     '''
     tinterp deriv=${SOURCES[1]} os=3 ds=0.025 ns=240 |
     window n1=1 f1=%d |
     spline o1=3 d1=0.025 n1=240
     ''' % (datum2/dz))

Flow('rgreen2_3','time2r_3 tds2r_3',
     '''
     tinterp deriv=${SOURCES[1]} os=0.425 ds=0.025 ns=335 |
     window n1=1 f1=%d |
     spline o1=0.425 d1=0.025 n1=335
     ''' % (datum2/dz))

Flow('redat2_3','redat2_2 sgreen2_3 rgreen2_3',
     '''
     reverse which=1 opt=i |
     kirdatsr sdatum=%g rdatum=%g aperture=300 taper=0 length=0.05
     sgreen=${SOURCES[1]} rgreen=${SOURCES[2]} |
     reverse which=1 opt=i
     ''' % (datum2,datum2))

# migration
Flow('dmig2_2','redat2_3 time1s_2 tds1s_2 time1r_2 tds1r_2',
     '''
     kirmigsr aperture=2 antialias=0.25 tau=0.0565 
     stable=${SOURCES[1]} sderiv=${SOURCES[2]}
     rtable=${SOURCES[3]} rderiv=${SOURCES[4]}
     ''')

# plot
######

scale2=2.

Flow('cut2','dmig2_2','window n1=%d f1=%d | scale dscale=%g'
     % ((3.-chunk1)/dz+1,(chunk1-datum1)/dz,scale2))
Flow('dmig2','dmig1_1 cut2','cat axis=1 ${SOURCES[1]}')

Result('dmig2',
       '''
       window f2=%d n2=%d |
       grey title="Recursive 3 steps" screenht=7 screenwd=13
       labelsz=6 titlesz=8 label2=Position pclip=96
       titlefat=4 labelfat=4
       ''' % (xmin/dx,(xmax-xmin)/dx+1))

Plot('zoom2','dmig2',
     '''
     window n1=251 f1=500 n2=751 f2=1250 |
     grey title="Redatum 3 steps" screenratio=0.7
     labelsz=10 titlesz=12 label2=Position pclip=90
     titlefat=6 labelfat=6
     ''')

# zoom
Result('zoom','zoom zoom0 zoom1 zoom2','TwoRows')

End()