from rsf.proj import *
import rsf.recipes.beg as beg

pl_Op=' labelsz=16 labelfat=6 titlesz=16 titlefat=6 crowd=0.7'

## Create RSF header for my binary data files 
color=1

tmin=3.0; #3.0
tmax=5.5; #5.5
dtau=0.003
color=1
colorPCK=5
linefat=15
nCMP="2"
path="../DATA/cmp" + nCMP + "/"

for prof in ('vN','eta'):

    bin = path + prof + '.bin'
    bin = prof + '.bin'
    Fetch(bin,'lorenzo',beg.server)

    Flow(prof,bin,
      '''       
      echo n1=2334 d1=%g o1=0 n2=1 d2=1 o2=1 in=${SOURCES[0]} label1="time" unit1="s" data_format=native_float |   
      window min1=%g | scale dscale=%g
      ''' % (dtau,tmin,(0.001,1)[prof=='eta']))
### EFFECTIVE PARAMETER FROM CONVENTIONAL ANALYSIS
Plot('vN',
     '''
     window max2=1 max1=%g|
     graph transp=y yreverse=y pad=n min2=0.995 max2=2.005 plotcol=%d wanttitle=n wantaxis=n plotfat=%d
     '''%(tmax,color,linefat)+pl_Op)
Plot('eta',
     '''
     window max2=1 max1=%g|
     graph transp=y yreverse=y pad=n min2=-0.5 max2=0.5 plotcol=%d wanttitle=n wantaxis=n plotfat=%d
     '''%(tmax,color,linefat)+pl_Op)

Flow('vH','vN eta','math e=${SOURCES[1]} output="input*sqrt(1+2*e)"')
Plot('vH',
     '''
     window max2=1 max1=%g|
     graph transp=y yreverse=y pad=n min2=0.995 max2=2.505 plotcol=%d wanttitle=n wantaxis=n plotfat=%d
     '''%(tmax,color,linefat)+pl_Op)

Result('effSEMB','vN vH eta','SideBySideAniso')
### INTERVAL PARAMETER FROM CONVENTIONAL ANALYSIS

for prof in ('vNi','etai'):

    bin = path + prof + '.bin'
    bin = prof + '.bin'
    Fetch(bin,'lorenzo',beg.server)

    Flow(prof,bin,
      '''       
      echo n1=2334 d1=%g o1=0 n2=1 d2=1 o2=1 in=${SOURCES[0]} label1="time" unit1="s" data_format=native_float |   
      window min1=%g | scale dscale=%g
      ''' % (dtau,tmin,(0.001,1)[prof=='etai']))

Plot('vNi',
     '''
     window max2=1 max1=%g|
     graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d wanttitle=n wantaxis=n plotfat=%d 
     '''%(tmax,color,linefat)+pl_Op)

Plot('etai',
     '''
     window max2=1 max1=%g|
     graph transp=y yreverse=y pad=n min2=-0.5 max2=0.5 plotcol=%d wanttitle=n wantaxis=n plotfat=%d 
     '''%(tmax,color,linefat)+pl_Op)

Flow('vHi','vNi etai','math e=${SOURCES[1]} output="input*sqrt(1+2*e)"')
Plot('vHi',
     '''
     window max2=1 max1=%g|
     graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d wanttitle=n wantaxis=n plotfat=%d 
     '''%(tmax,color,linefat)+pl_Op)

Result('intSEMB','vNi vHi etai','SideBySideAniso')

#### READ THE DATA

for data in ('cmp','cmpFX'):
    bin = path + data+'.bin'
    bin = data + '.bin'
    Fetch(bin,'lorenzo',beg.server)

    nx = (465,59)[data=='cmp']
    dx = (0.0125,0.1)[data=='cmp']

    nt = (1334,1267)[nCMP=='1']
    ot = (3.0,3.201)[nCMP=='1']

    Flow(data,bin,
         '''    
         echo n1=%d d1=%g o1=%g n2=%d o2=0.475 d2=%g in=${SOURCES[0]} label1="time" unit1="s" label2="offset" unit2="km" data_format=native_float |
                 pad beg1=101
         '''% (nt,dtau,ot,nx,dx))

Plot('cmp',
     '''
     grey title="data Original (a)" pclip=80
     ''' + pl_Op)

Plot('cmpFX',
       '''
       grey title="data FX (a)" pclip=80
       ''' + pl_Op)

##########################################################################################################################
# RADON TRANFORMATION OF THE DATA
#################################'

kp=4
dp= 0.00225/kp
np=151*kp
op=0.05
pmax=0.32
pmin=0.09

tstart=7.0
#slope0=5.425
slope0=5.0

Flow('cmp_R','cmp',
     '''        
     put unit1="s" label2="p" unit2="s/km" | pad beg1=101 | 
     radon np=%d dp=%g p0=%g verb=y inv=n nx=59 ox=0.436 dx=0.1  |
         mutter t0=%g slope0=-%g x0=%g inner=y 
     '''%(np,dp,op,tstart,slope0,op))
Plot('cmp_R','grey title="data Radon (b)" label1="\F10 t\F3" unit1="s" label2="p" unit2="s/km" pclip=80 ' + pl_Op)

Flow('cmpFX_R','cmpFX',
     '''        
         pad beg1=101 | 
     radon np=%d dp=%g p0=%g verb=y inv=n nx=465 ox=0.436 dx=0.0125 | 
     mutter t0=%g slope0=-%g x0=%g inner=y |
     put  unit1="s" label2="p" unit2="s/km"
     '''%(np,dp,op,tstart,slope0,op))
Plot('cmpFX_R','grey title="data FX Radon (b)" pclip=80 label1="\F10 t\F3" ' + pl_Op )

Result('cmp','cmp cmp_R','SideBySideAniso')
Result('cmpFX','cmpFX cmpFX_R','SideBySideAniso')

## COMPUTATION OF AN INITIAL SLOPES FIELD
vN0=1.6
vH0=1.8


Flow('vN0','vN',
     '''
     math  output="%g*1+0.005*(x1-3.2)" label1="\F10 t\F3" unit1="s" 
     ''' %(vN0))
Flow('vH0','vH',
     '''
     math  output="%g*1+0.005*(x1-3.2)" label1="\F10 t\F3" unit1="s" 
     ''' %(vH0))

Flow('eta0','vN0 vH0','math vh=${SOURCES[1]} output=".5*(vh^2/input^2-1)"')

Flow('vH0mat','vH0',
    '''
     spray n=%g d=%g o=%g axis=2  |             
     put label2="p" unit2="s/km" 
     ''' % (np,dp,op))
Flow('vN0mat','vN0',
     '''
     spray n=%g d=%g o=%g axis=2 |
     put label2="p" unit2="s/km" 
     ''' % (np,dp,op))
#Plot('vH0mat','grey color=j scalebar=y bias=2.2')
#Plot('vN0mat','grey color=j scalebar=y bias=2.0')
#Result('mat','vH0mat vN0mat','SideBySideAniso')

a = dp/dtau
Flow('dip0','cmpFX_R vN0mat vH0mat',
     '''
         window min1=%g | 
     math vN=${SOURCES[1]} vH=${SOURCES[2]} output="%g*(-x1*vN^2*x2)/((1-x2^2*vH^2)^(1/2))/(1-x2^2*(vH^2-vN^2))^(3/2)"
     ''' %(tmin, a) )
Flow('dip0w','dip0','window min1=%g max1=%g'%(tmin,tmax))

Flow('dip00','dip0 vN0 eta0',
     '''
     itaupmo interval=n velocity=${SOURCES[1]} eta=${SOURCES[2]} | 
     window min1=%g max1=%g    
     '''%(tmin,tmax))
Plot('dip0','grey color=j scalebar=y label1="\F10 t\F3" title=dip0')
Plot('dip00','grey color=j scalebar=y label1="\F10 t\F3" title=dip00')
Result('dip0','dip0 dip00','SideBySideAniso')

## DIP ESTIMATION 
#L1=150
#L2=60
L1=200
L2=100
# slope estimation (dimensionless)
Flow('dips','cmpFX_R dip0w',
     '''
     window  min1=%g max1=%g |
     dip both=y idip=${SOURCES[1]} rect1=%d rect2=%d order=5 pmax=0 transp=y 
     '''%(tmin,tmax,L1,L2))
Flow('dip','dips','stack axis=3 scale=0.5,-0.5 norm=n')
Flow('dip1','dips','window n3=1')



#Flow('dip','data','dip rect1=30 rect2=5 order=3 pmax=0 transp=y')
Plot('dip','grey color=j title="slope (a)" scalebar=y bias=-0.5 label1="\F10 t\F3" polarity=y pclip=80' + pl_Op)
# scaled derivatives(dimensionless)

smoothdert = 100;
smoothderp = 10;

Flow('dipx','dip','transp | deriv | smooth rect1=%d rect2=%d | transp'%(smoothderp,smoothdert))

Flow('dipt','dip',
     '''
     deriv |
     sfsmooth rect1=%d rect2=%d
     '''%(smoothdert,smoothderp))

Flow('curv','dip dipx dipt',
     '''
     math dx=${SOURCES[1]} dt=${SOURCES[2]} output="input*dt+dx" |
     sfsmooth rect1=%d rect2=%d   
     '''%(smoothdert,smoothderp))

Plot('curv','grey color=j title="curvature (b)" scalebar=y bias=-.002 maxval=0 label1="\F10 t\F3" polarity=y pclip=90' + pl_Op)

Result('data','cmpFX_R dip curv','SideBySideAniso')

# windowing of the data... take just the portion when curvature estimate is reliable enough
#sfwindow < in.rsf > out.rsf verb=n squeeze=y j#=(1,...) d#=(d1,d2,...) min#=(o1,o2,,...) max#=(o1+(n1-1)*d1,o2+(n1-1)*d2,,...)

Flow('datawin','cmpFX_R','window min1=%f max1=%f min2=%f max2=%f'%(tmin,tmax,pmin,pmax))
Flow('dipwin','dip','window min2=%f max2=%f'%(pmin,pmax))
Flow('dipwin1','dip1','window min2=%f max2=%f'%(pmin,pmax))
Flow('curvwin','curv','window min2=%f max2=%f'%(pmin,pmax))
Flow('diptwin','dipt','window min2=%f max2=%f'%(pmin,pmax))
Plot('datawin','grey title="data (a)"  label1="\F10 t\F3" ' + pl_Op)
Plot('dipwin','grey title="slope (b)"  color=j scalebar=y bias=-0.5 label1="\F10 t\F3" polarity=y pclip=80 ' + pl_Op)
Plot('curvwin','grey title="curvature (c)" color=j scalebar=y bias=-.002 maxval=0 label1="\F10 t\F3" polarity=y pclip=90 ' + pl_Op)
Result('datawin','datawin dipwin curvwin','SideBySideAniso')

tau0str="\F9 t\_\F3 0\^(\F9 t\F3 ,p) \F3 "

# slope independent pNMO (isotropic)
Flow('nmo vel2','datawin dipwin diptwin','ptaupmo dip=${SOURCES[1]} dipt=${SOURCES[2]} vel2=${TARGETS[1]}')
Plot('nmo','grey title="ISO NMO (b)" grid2=y gridfat=3 gridcol=5 g2num=.5  label1="\F10 t\F3" ' + pl_Op ) ### slides

# slope independent pNMO (VTI)
Flow('nmoVTI tau0','datawin dipwin curvwin  ','ptaupmoVTI dip=${SOURCES[1]} ddip=${SOURCES[2]} tau0=${TARGETS[1]}')
Plot('nmoVTI','grey title="VTI NMO (c)" grid2=y gridfat=3 gridcol=5 g2num=.5 label1="\F10 t\F3" ' + pl_Op) ### slides
Result('dataNMO','datawin nmo nmoVTI','SideBySideAniso')

Plot('tau0','tau0','sfgrey color=j scalebar=y pclip=90 bias=4 minval=2 maxval=6  title="'+tau0str+' (b)" label1="\F10 t\F3" '+ pl_Op)
Plot('tau01','tau0','contour nc=50 c0=%g dc=%g plotcol=5 plotfat=3 wanttitle=n wantaxis=n' % (tmin,(tmax-tmin)/30) + pl_Op)
Plot('tau02','datawin tau01','Overlay')
Result('dataNMO1','tau02 tau0 nmoVTI','SideBySideAniso')

Flow('curvtwin','curvwin','deriv | smooth rect1=%d rect2=%d' % (smoothdert,smoothderp))

## VELOCITY ESTIMATION: EFFECTIVE VELOCITIES

Flow('vNeff vHeff etaeff','tau0 datawin dipwin curvwin',
     '''
     pveltranVTI method=e velH=${TARGETS[1]} eta=${TARGETS[2]}
     v0=1.0 dv=0.005 nv=201 nvh=301  nw=4 map=n
     cmp=${SOURCES[1]} dip=${SOURCES[2]} curv=${SOURCES[3]} 
     ''')
Flow('vNmapE vHmapE etamapE','tau0 dipwin curvwin',
     '''
     pveltranVTI method=e velH=${TARGETS[1]} eta=${TARGETS[2]}
     v0=1.0 dv=0.005 nv=201 nvh=301 nw=4 map=y
     dip=${SOURCES[1]} curv=${SOURCES[2]} 
     ''')

Plot('vNmapE','grey color=j scalebar=y clip=.5 minval=0 maxval=3 bias=1.5 title="effective \F5 V\_N \^ \F3 (a)" ' + pl_Op)
Plot('vHmapE','grey color=j scalebar=y clip=.5 minval=0 maxval=3 bias=1.5 title="effective \F5 V\_H \^ \F3 (b)" ' + pl_Op)
Plot('etamapE','grey color=j scalebar=y  clip=0.25 minval=-0.5 maxval=0.5 bias=0 title="effective \F9 h \F3 (c)" ' + pl_Op)
Result('mapE','vNmapE vHmapE etamapE','SideBySideAniso')

Plot('vNeff','vNeff',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s 
     title="effective \F5 V\_N \^ \F3 (a)" label1="\F10 t\F3"  
    ''' + pl_Op)

Plot('vNeff2','vNeff vN','Overlay')

Plot('vHeff','vHeff',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s 
     title="effective   \F5 V\_H \^ \F3 (b)" label1="\F10 t\F3"  
     ''' + pl_Op)

Plot('vHeff2','vHeff vH','Overlay')

Plot('etaeff','etaeff',
    '''
     envelope |
     grey allpos=y label2='' unit2=''
     title="effective  \F9 h \F3  (c)" label1="\F10 t\F3" 
    ''' + pl_Op)

Plot('etaeff2','etaeff eta','Overlay')

Result('eff','vNeff2 vHeff2 etaeff2','SideBySideAniso')

## VELOCITY ESTIMATION: INTERVAL VELOCITIES 
## STRIPPING PROCESSING

Flow('vNint vHint etaint','tau0 datawin diptwin curvtwin ',
     '''
     pveltranVTI method=s velH=${TARGETS[1]} eta=${TARGETS[2]} 
     v0=1.0 dv=0.010 nv=301 interval=y nw=4 map=n
     cmp=${SOURCES[1]} dipt=${SOURCES[2]} curvt=${SOURCES[3]} 
     ''')
Flow('vNmapS vHmapS etamapS','tau0 diptwin curvtwin ',
     '''
     pveltranVTI method=s velH=${TARGETS[1]} eta=${TARGETS[2]} 
     v0=1.0 dv=0.010 nv=301 interval=y nw=4 map=y
     dipt=${SOURCES[1]} curvt=${SOURCES[2]} 
     ''')

Plot('vNmapS','grey color=j scalebar=y clip=1.5 minval=0 maxval=4 bias=2 title="interval \F5 V\_N \^ \F3 (a)" ' + pl_Op)
Plot('vHmapS','grey color=j scalebar=y clip=1.5 minval=0 maxval=4 bias=2 title="interval \F5 V\_H \^ \F3 (b)" ' + pl_Op)
Plot('etamapS','grey color=j scalebar=y  clip=0.4 minval=-0.5 maxval=0.5 bias=0 title="interval \F9 h \F3 (c)" ' + pl_Op)
Result('mapS','vNmapS vHmapS etamapS','SideBySideAniso')

Plot('vNint','vNint',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_N \^ \F3 (a)" label1="\F10 t\F3" 
    ''' + pl_Op)

Flow('vNintpk','vNint','envelope | scale axis=2 | pick rect1=10 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('vNintpk','graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op)


Plot('vHint','vHint',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_H \^ \F3 (b)" label1="\F10 t\F3" 
    ''' + pl_Op)

#Plot('vHint2','vHint vHintE','Overlay')

Plot('etaint','etaint',
    '''
     envelope |
     grey allpos=y label2='' unit2=''
     title="interval \F9 h \F3 (c)" label1="\F10 t\F3" 
    ''' + pl_Op)


## VELOCITY ESTIMATION: INTERVAL VELOCITIES
## DIX PROCESSING

Flow('vNint_dix vHint_dix etaint_dix','tau0 datawin dipwin curvwin diptwin curvtwin ',
     '''
     pveltranVTI method=d velH=${TARGETS[1]} eta=${TARGETS[2]} 
     v0=1.0 dv=0.010 nv=301 interval=y nw=4 map=n
     cmp=${SOURCES[1]} dip=${SOURCES[2]} curv=${SOURCES[3]} dipt=${SOURCES[4]} curvt=${SOURCES[5]} 
     ''')


Plot('vNint_dix','vNint_dix',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_N \^ \F3 (a)" label1="\F10 t\F3" 
    ''' + pl_Op)

Flow('vNintpk_dix','vNint_dix','envelope | scale axis=2 | pick rect1=10 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('vNintpk_dix','graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op)

Plot('vHint_dix','vHint_dix',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_H \^ \F3 (b)" label1="\F10 t\F3" 
    ''' + pl_Op)

#Plot('vHint_dix2','vHint_dix vHintE','Overlay')

Plot('etaint_dix','etaint_dix',
    '''
     envelope |
     grey allpos=y label2='' unit2=''
     title="interval \F9 h \F3 (c)" label1="\F10 t\F3" 
    '''+ pl_Op)

## CLAREBOUT-FOWLER PROCESSING  #######################

Flow('tau0t','tau0',
     '''
     deriv | smooth rect1=%d rect2=%d
     '''
     % (smoothdert,smoothderp))

Flow('vNint_fowler vHint_fowler etaint_fowler','tau0 datawin tau0t diptwin',
     '''
     pveltranVTI method=f velH=${TARGETS[1]} eta=${TARGETS[2]}
     cmp=${SOURCES[1]} tau0t=${SOURCES[2]} dipt=${SOURCES[3]} 
      v0=1.0 dv=0.010 nv=301 interval=y nw=4 map=n
     ''')

Flow('vNintpk_fowler','vNint_fowler','envelope | scale axis=2 | pick rect1=10 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('vNintpk_fowler','graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op)

Plot('vNint_fowler','vNint_fowler',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_N \^ \F3  (a)" label1="\F10 t\F3" 
    ''' + pl_Op)

#Plot('vNint_fowler2','vNint_fowler vNintE','Overlay')

Plot('vHint_fowler','vHint_fowler',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_H \^ \F3 (b)" label1="\F10 t\F3" 
    ''' + pl_Op)

#Plot('vHint_fowler2','vHint_fowler vHintE','Overlay')

Plot('etaint_fowler','etaint_fowler',
    '''
     envelope |
     grey allpos=y label2='' unit2=''
     title="interval \F9 h \F3 (c)" label1="\F10 t\F3" 
    ''' + pl_Op)


#Plot('etaint_fowler2','etaint_fowler etaintE','Overlay')
#Result('temp','tempeff tempint','SideBySideAniso')

#RESULTS!!!

Plot('vNint2','vNint vNintpk_dix','Overlay')
Plot('vNint2_dix','vNint_dix vNintpk_dix','Overlay')
Plot('vNint2_fowler','vNint_fowler vNintpk_dix','Overlay')

Result('int','vNint2 vHint etaint','SideBySideAniso')
Result('int_dix','vNint2_dix vHint_dix etaint_dix','SideBySideAniso')
Result('int_fowler','vNint2_fowler vHint_fowler etaint_fowler','SideBySideAniso')
Result('int_VNMO','vNint2_dix vNint2_fowler vNint2','SideBySideAniso')

# slope estimation for painting (dimensionless)
Flow('dipPAINT','cmpFX_R dip1',
     '''
     window  min1=%g max1=%g |
     dip idip=${SOURCES[1]} rect1=%d rect2=%d order=4 pmax=0 transp=y 
     '''%(tmin,tmax,50,20))

Flow('dipPAINTwin','dipPAINT','window min2=%f max2=%f'%(pmin,pmax))
Plot('dipPAINTwin','dipPAINTwin','grey title="slope (b)" color=j scalebar=y label1="\F10 t\F3" color=j polarity=y bias=-0.5 pclip=80 maxval=0 ' + pl_Op)

Flow('seed','tau0','window n2=1')
Flow('t0paint','dipPAINTwin seed','pwpaint seed=${SOURCES[1]} verb=y')

Plot('t0paint1','t0paint','contour nc=50 c0=%g dc=%g plotcol=5 plotfat=3 wanttitle=n wantaxis=n' % (tmin,(tmax-tmin)/30) + pl_Op)
Result('paint','datawin t0paint1','Overlay')

# FLATTENING

Flow('flat','datawin t0paint','iwarp warp=${SOURCES[1]}')
Plot('nmoflat','flat','grey title="flattened (d)" grid2=y gridfat=3 gridcol=5  g2num=.5 label1="\F10 t\F3" ' + pl_Op)

Plot('t0paint','grey color=j scalebar=y clip=1.4 bias=4.5 minval=3.0  label1="\F10 t\F3" maxval=6 title="'+tau0str+' (c)" ' + pl_Op)
Result('flat','t0paint nmoflat','SideBySideAniso')
Plot('t0paint_con','t0paint','contour nc=38 c0=%g dc=%g plotcol=7 plotfat=6 scalebar=y wanttitle=n wantaxis=n' % (tmin,(tmax-tmin)/30) + pl_Op)

Plot('t0paint_new','t0paint','sfgrey color=j scalebar=y clip=1.4 bias=4.5 minval=3.0 title="painted \F9 t\_\F3 0 \^(c) " ' + pl_Op)
Plot('t0paintw','t0paint_new t0paint_con','Overlay')


Flow('dipPAINTt','dipPAINTwin','deriv | smooth rect1=%d rect2=%d '%(smoothdert,smoothderp))
Flow('t0paintt','t0paint','deriv | smooth rect1=%d rect2=%d' % (smoothdert,smoothderp))



Result('dataP','datawin dipPAINTwin t0paint nmoflat','SideBySideAniso')
Result('dataPw','datawin dipPAINTwin t0paintw nmoflat','SideBySideAniso')

# VELOCITY
#Flow('vNint_painting vHint_painting etaint_painting  vNmapP vHmapP etamapP','datawin t0paint t0paintt dipPAINTt',
#     '''
#     pveltranVTI_FOWLER tau0=${SOURCES[1]} tau0t=${SOURCES[2]} dipt=${SOURCES[3]} velH=${TARGETS[1]} eta=${TARGETS[2]}
#      vNmap=${TARGETS[3]} vHmap=${TARGETS[4]} etamap=${TARGETS[5]} v0=1.0 dv=0.010 nv=301 nw=4
#     ''')
Flow('vNint_painting vHint_painting etaint_painting','t0paint datawin t0paintt dipPAINTt',
     '''
     pveltranVTI method=f velH=${TARGETS[1]} eta=${TARGETS[2]}
     cmp=${SOURCES[1]} tau0t=${SOURCES[2]} dipt=${SOURCES[3]} 
     v0=1.0 dv=0.010 nv=301 interval=y nw=4 map=n
     ''')
Flow('vNmapP vHmapP etamapP','t0paint t0paintt dipPAINTt',
     '''
     pveltranVTI method=f velH=${TARGETS[1]} eta=${TARGETS[2]}
     tau0t=${SOURCES[1]} dipt=${SOURCES[2]} 
     v0=1.0 dv=0.010 nv=301 nw=4 map=y
     ''')

Plot('vNmapP','grey color=j scalebar=y clip=1.5 minval=0 maxval=4 bias=2 title="interval \F5 V\_N \^ \F3 (a)" ' + pl_Op)
Plot('vHmapP','grey color=j scalebar=y clip=2 minval=0 maxval=4 bias=2 title="interval \F5 V\_H \^ \F3 (b)" ' + pl_Op)
Plot('etamapP','grey color=j scalebar=y  clip=0.5 minval=-0.5 maxval=0.5 bias=0 title="interval \F9 h \F3 (c)" ' + pl_Op)
Result('mapP','vNmapP vHmapP etamapP','SideBySideAniso')

Flow('vNintpk_painting','vNint_painting','envelope | scale axis=2 | pick rect1=10 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('vNintpk_painting','graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op )

Plot('vNint_painting','vNint_painting',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_N \^ \F3 (a)" label1="\F10 t\F3" 
    ''' + pl_Op)

Plot('vNint2_painting','vNint_painting vNintpk_painting vNi','Overlay')

###

Plot('vHint_painting','vHint_painting',
    '''
     envelope |
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_H \^ \F3 (b)" label1="\F10 t\F3" 
    ''' + pl_Op)

Flow('vHintpk_painting','vHint_painting','envelope | scale axis=2 | pick rect1=10 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('vHintpk_painting','graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op )

Plot('vHint2_painting','vHint_painting vHintpk_painting vHi','Overlay')
##
Plot('etaint_painting','etaint_painting',
    '''
     envelope |
     grey allpos=y label2='' unit2=''
     title="interval \F9 h \F3 (c)" label1="\F10 t\F3" 
    ''' + pl_Op)

Flow('etaintpk_painting','etaint_painting','envelope | scale axis=2 | pick rect1=10 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('etaintpk_painting','graph transp=y yreverse=y pad=n min2=-0.5 max2=0.5 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op )

Plot('etaint2_painting','etaint_painting etaintpk_painting etai','Overlay')

Result('int-painting','vNint2_painting vHint2_painting etaint2_painting ','SideBySideAniso')
Result('interval','vNint2_painting vHint2_painting','SideBySideAniso')


## masking poor SNR areas
t0_mapP=5.5
op_mapP=0.2
slope0_mapP=7

Flow('vNmapP_mask','vNmapP','mutter half=n t0=%g slope0=-%g x0=%g abs=n nan=y inner=y'%(t0_mapP,slope0_mapP,op_mapP))
Flow('vHmapP_mask','vHmapP',
        '''
        mutter half=n t0=%g slope0=-%g x0=%g inner=y nan=y abs=n |
        mutter half=n t0=%g slope0=-1000 x0=%g inner=n nan=y abs=n 
        '''%(t0_mapP,slope0_mapP,op_mapP,t0_mapP,op_mapP+0.005))
Flow('etamapP_mask','etamapP',
        '''
        mutter half=n t0=%g slope0=-%g x0=%g inner=y nan=y abs=n |
        mutter half=n t0=%g slope0=-45 x0=%g inner=n nan=y abs=n
        '''%(t0_mapP,slope0_mapP,op_mapP,t0_mapP,op_mapP))


Plot('vNmapP_mask','grey color=j scalebar=y clip=1.5 minval=0 maxval=4 bias=2 title="interval \F5 V\_N \^ \F3 (a)" ' + pl_Op)
Plot('vHmapP_mask','grey color=j scalebar=y clip=2 minval=0 maxval=4 bias=2 title="interval \F5 V\_H \^ \F3 (b)" ' + pl_Op)
Plot('etamapP_mask','grey color=j scalebar=y  clip=0.5 minval=-0.5 maxval=0.5 bias=0 title="interval \F9 h \F3 (c)" ' + pl_Op)
Result('mapP_mask','vNmapP_mask vHmapP_mask etamapP_mask','SideBySideAniso')

Flow('vNint_painting_mask','flat vNmapP_mask','map2coh  map=${SOURCES[1]} v0=1.0 dv=0.010 nv=301 min2=0.12')
Flow('vHint_painting_mask','flat vHmapP_mask','map2coh  map=${SOURCES[1]} v0=1.0 dv=0.010 nv=301 ')
Flow('etaint_painting_mask','flat etamapP_mask','map2coh  map=${SOURCES[1]} v0=-0.5 dv=0.01 nv=101  ')

Plot('vNint_painting_mask','vNint_painting_mask',
    '''
     envelope | math output="input^1" | 
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_N \^ \F3 (a)" label1="\F10 t\F3" 
    ''' + pl_Op)

Flow('vNintpk_painting_mask','vNint_painting_mask','envelope | math output="input" | scale axis=2 | pick rect1=20 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('vNintpk_painting_mask','graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op )
Plot('vNint2_painting_mask','vNint_painting_mask vNintpk_painting_mask','Overlay')
Plot('vNint3_painting_mask','vNint_painting_mask vNintpk_painting_mask vNi','Overlay')

Plot('vHint_painting_mask','vHint_painting_mask',
    '''
     envelope | math output="input^1" | 
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F5 V\_H \^ \F3 (b)" label1="\F10 t\F3" 
    ''' + pl_Op)

Flow('vHintpk_painting_mask','vHint_painting_mask','envelope | math output="input" | scale axis=2 | pick rect1=20 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('vHintpk_painting_mask','graph transp=y yreverse=y pad=n min2=0.995 max2=4.005 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op )
Plot('vHint2_painting_mask','vHint_painting_mask vHintpk_painting_mask','Overlay')
Plot('vHint3_painting_mask','vHint_painting_mask vHintpk_painting_mask vHi','Overlay')

Plot('etaint_painting_mask','etaint_painting_mask',
    '''
     envelope | math output="input^1" | 
     grey allpos=y label2=Velocity unit2=km/s
     title="interval \F9 h \F3 (c)" label1="\F10 t\F3" 
    ''' + pl_Op)

Flow('etaintpk_painting_mask','etaint_painting_mask','envelope | math output="input" | scale axis=2 | pick rect1=20 | window min1=%g| remap1 pattern=$SOURCE order=2'%3.5)
Plot('etaintpk_painting_mask','graph transp=y yreverse=y pad=n min2=-0.5 max2=0.5 plotcol=%d plotfat=10 wanttitle=n wantaxis=n '%colorPCK + pl_Op )
Plot('etaint2_painting_mask','etaint_painting_mask etaintpk_painting_mask','Overlay')
Plot('etaint3_painting_mask','etaint_painting_mask etaintpk_painting_mask etai','Overlay')

Flow('eta_divn','vHintpk_painting_mask vNintpk_painting_mask','divn rect1=50 den=${SOURCES[1]} | math output=".5*(input^2-1)" ')
Plot('eta_divn','graph transp=y yreverse=y pad=n min2=-0.5 max2=0.5 plotcol=3 plotfat=10 wanttitle=n wantaxis=y ' + pl_Op )
Result('etaVS','eta_divn etai etaintpk_painting_mask','Overlay')

Result('int-painting-mask','vNint3_painting_mask vHint3_painting_mask etaint3_painting_mask ','SideBySideAniso')
Result('interval-mask','vNint3_painting_mask vHint3_painting_mask ','SideBySideAniso')

Result('dataPw+vel','datawin dipPAINTwin vNint2_painting_mask vHint2_painting_mask','SideBySideAniso')

## Reconstructed model

Flow('cmpFX_R_rec','cmpFX_R vNintpk_painting_mask eta_divn',
     '''
     window f2=0 n2=1 min1=%g max1=%g |
     spray n=%g d=%g o=%g axis=2 |
     itaupmo interval=y velocity=${SOURCES[1]} eta=${SOURCES[2]} |
     window min2=%g max2=%g |
     put label2="p" unit2="s/km" |
         mutter half=n t0=%g slope0=-%g x0=%g inner=y  | window squeeze=y
     ''' %(tmin,tmax,np,dp,op,pmin,pmax,tstart,slope0,op))

Plot('cmpFX_R_rec','grey title="reconstructed gather" pclip=80 ' + pl_Op )
Result('cmp_rec','datawin cmpFX_R_rec','SideBySideAniso')

## Physically flattened gather
Flow('cmpFX_R_pmo','datawin vNint_painting_mask vHint_painting_mask',
'''
taupmo velocity=${SOURCES[0]} interval=y
''')
Plot('cmpFX_R_pmo','grey title="physically flattened" grid2=y gridfat=3 gridcol=5  g2num=.5 ' + pl_Op )
Result('flatVS','nmoflat cmpFX_R_pmo','SideBySideAniso')

## Mesaure X/D in the dataset

Flow('D','vN','window max1=%f | math output="x1/2*input" | spray n=%g d=%g o=%g axis=2 | window min2=%f max2=%f'%(tmax,np,dp,dp,pmin,pmax))

Flow('X','dipwin t0paint','iwarp warp=${SOURCES[1]} | math output="(-1)*input*%f/%f"' %(dtau,dp))
Flow('XD','X D','math D=${SOURCES[1]} output="input/(D+1e-6)"')
Flow('theta','XD','math output="atan(input)/3.14*180" ')

Result('XD','XD','grey bias=0.75 clip=0.55 scalebar=y color=j minval=0.0 maxval=1.55')
Result('theta','theta',' grey bias=22.5 clip=20 scalebar=y color=j minval=0.0 maxval=45.0')




End()