## 
 # Angle-gathers after wave-equation migration (ziggy model)
 ##

from rsf.proj import *
import pplot

# ------------------------------------------------------------
# MIGRATION parameters
# ------------------------------------------------------------
par = {
    'ox':10925,'dx':150, 'nx':500,
    'oz':6100, 'dz':25,
    'oh':-3000,'dh':150
    }

# ------------------------------------------------------------
# from ft to km
ft2km = 0.3024/1000

par['ox']=par['ox']*ft2km
par['dx']=par['dx']*ft2km

par['oz']=par['oz']*ft2km
par['dz']=par['dz']*ft2km

par['oh']=par['oh']*ft2km
par['dh']=par['dh']*ft2km

# ------------------------------------------------------------

par['zmin']=2
par['zmax']=9

par['xmin']=4
par['xmax']=21

IRATIO = ' labelsz=6 screenratio=0.5 screenht=7 min2=%g max2=%g ' % (par['xmin'],par['xmax'])
CRATIO = ' labelsz=5 screenratio=3.5 screenht=10.5 '
LABEL = " wantaxis=n "
# ------------------------------------------------------------

def igrey(custom,par):
    return '''
    grey labelrot=n wantaxis=y wanttitle=n grid=y gridcol=6
    title="" pclip=97 label1="z" unit1=km label2="x" unit2=km
    min1=%g max1=%g
    %s
    ''' % (par['zmin'],par['zmax'],custom)

def p1x3(plot,p0,p1,p2,ys,xs,x0,x1,x2):
    j0 = '_' + plot+ p0
    j1 = '_' + plot+ p1
    j2 = '_' + plot+ p2

    Plot(j0,p0,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x0))
    Plot(j1,p1,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x1))
    Plot(j2,p2,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x2))

    Plot  (plot,[j0,j1,j2],'Overlay')
    Result(plot,[j0,j1,j2],'Overlay')

def p1x4(plot,p0,p1,p2,p3,ys,xs,x0,x1,x2,x3):
    j0 = '_' + plot+ p0
    j1 = '_' + plot+ p1
    j2 = '_' + plot+ p2
    j3 = '_' + plot+ p3

    Plot(j0,p0,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x0))
    Plot(j1,p1,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x1))
    Plot(j2,p2,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x2))
    Plot(j3,p3,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,x3))

    Plot  (plot,[j0,j1,j2,j3],'Overlay')
    Result(plot,[j0,j1,j2,j3],'Overlay')

def p1x6(plot,p0,p1,p2,p3,p4,p5,ys,xs,xc):
    j0 = '_' + plot+ p0
    j1 = '_' + plot+ p1
    j2 = '_' + plot+ p2
    j3 = '_' + plot+ p3
    j4 = '_' + plot+ p4
    j5 = '_' + plot+ p5

    Plot(j0,p0,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,   0-1))
    Plot(j1,p1,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,  xc-1))
    Plot(j2,p2,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,2*xc-1))
    Plot(j3,p3,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,3*xc-1))
    Plot(j4,p4,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,4*xc-1))
    Plot(j5,p5,'Overlay',vppen='yscale=%f xscale=%f ycenter=0 xcenter=%f'% (ys,xs,5*xc-1))

    Plot  (plot,[j0,j1,j2,j3,j4,j5],'Overlay')
    Result(plot,[j0,j1,j2,j3,j4,j5],'Overlay')

# ------------------------------------------------------------

# ------------------------------------------------------------
# from SGY to RSF
# ------------------------------------------------------------
velo = 'sigsbee2a_migvel.sgy'
vstr = 'sigsbee2a_stratigraphy.sgy'

Fetch(velo,'sigsbee')
Fetch(vstr,'sigsbee')

Flow('zvelo tzvelo ./vhead ./bvhead',velo,
     '''
     segyread
     tape=$SOURCE
     tfile=${TARGETS[1]}
     hfile=${TARGETS[2]}
     bfile=${TARGETS[3]}
     ''',stdin=0)
Flow('zvstr tzvstr ./shead ./bshead',vstr,
     '''
     segyread
     tape=$SOURCE
     tfile=${TARGETS[1]}
     hfile=${TARGETS[2]}
     bfile=${TARGETS[3]}
     ''',stdin=0)

# ------------------------------------------------------------
# VELOCITY
# ------------------------------------------------------------
# stratigraphic slowness
Flow('vstr','zvstr',
     '''
     scale rscale=%g |
     put d1=%g o2=%g d2=%g label1=z label2=x
     '''% (ft2km,25*ft2km, 10025*ft2km, 25*ft2km))
Flow('vstrpad','vstr',
     '''
     window n1=1 f1=1200 |
     spray axis=1 n=100 o=0 d=%g
     ''' % (25*ft2km))
Flow('s0','vstr vstrpad',
     '''
     cat axis=1 space=n ${SOURCES[1]} |
     window |
     math "output=1/input" |
     transp |
     spray axis=2 n=1 o=0 d=1 |
     put label2=y |
     window squeeze=n f3=244
     ''' % par)

# migration slowness
Flow('velo','zvelo',
     '''
     scale rscale=%g |
     put d1=%g o2=%g d2=%g label1=z label2=x
     '''% (ft2km,25*ft2km, 10025*ft2km, 37.5*ft2km))
Flow('velopad','velo',
     '''
     window n1=1 f1=1200 |
     spray axis=1 n=100 o=0 d=%g
     ''' % (25*ft2km))
Flow('s1','velo velopad',
     '''
     cat axis=1 space=n ${SOURCES[1]} |
     window |
     math "output=1/input" |
     transp |
     spray axis=2 n=1 o=0 d=1 |
     put label2=y |
     window squeeze=n f3=244
     ''')

# wrong migration velocity
Flow('sm','s1',
     '''
     window squeeze=n f3=350 |
     math output=0.1 |
     pad beg3=350    |
     math output=input+1 |
     smooth rect3=31
     ''')
Flow('s2','s1 sm',
     'math s=${SOURCES[0]} m=${SOURCES[1]} output=s*m')

# ------------------------------------------------------------
# MIGRATION
# ------------------------------------------------------------

SLOlist = '0','2'
CIGlist = '7','9','11','13','15','17'
IMGlist = 'x','t'

cigpar = {'npt':500,'opt':4000,'dpt':80,
           'na':400, 'oa':-2.0, 'da':0.01, 'mina':0, 'maxa':1.3 }
cigpar['opt']=cigpar['opt']*ft2km
cigpar['dpt']=cigpar['dpt']*ft2km

# ------------------------------------------------------------

# CIG to slant-stack
def tcig2ssk():
    return '''
    slant adj=y p0=%g np=%d dp=%g
    ''' % (cigpar['opt'],cigpar['npt'],cigpar['dpt'])
def hcig2ssk():
    return '''
    slant adj=y p0=%g np=%d dp=%g
    ''' % (cigpar['oa'],cigpar['na'],cigpar['da'])

# slant-stack to angle
def tssk2ang():
    return '''
    tshift cos=y a0=0 na=150 da=0.45
    velocity=${SOURCES[1]} dip=${SOURCES[2]}
    '''
def hssk2ang():
    return '''
    tan2ang a0=0 na=150 da=0.45
    '''
# ------------------------------------------------------------

for i in SLOlist:
    for k in IMGlist:
        SRcig = 'SRC' + i + k
        Fetch(     SRcig+'.hh','sigsbee')

        if(k=='t'):
            Flow(SRcig,SRcig+'.hh',
                 '''
                 window squeeze=n |
                 put o1=%(ox)g d1=%(dx)g
                     o3=%(oz)g d3=%(dz)g
                 ''' % par)
        if(k=='x'):
            Flow(SRcig,SRcig+'.hh',
                 '''
                 window squeeze=n |
                 put o1=%(ox)g d1=%(dx)g
                     o3=%(oz)g d3=%(dz)g
                     o4=%(oh)g d4=%(dh)g
                 ''' % par)

        if(0):
            cig = 'cig' + i + k
            Flow(cig,SRcig,'transp plane=15 memsize=500 | window')
            env = 'env' + i + k
            Flow(env,cig,'envelope')
            pik = 'pik' + i + k
            Flow(pik,env,'pick rect1=10 rect2=100')
            slc = 'slc' + i + k
            Flow(slc,[cig,pik],'slice pick=${SOURCES[1]} | window')
            Result(slc,igrey('pclip=95'+IRATIO+LABEL,par))

# ------------------------------------------------------------
# RESULTS
# ------------------------------------------------------------
for i in SLOlist:
    slo = 's' + i
    Plot  (slo,'window      | transp |'
           +igrey('pclip=99 allpos=y bias=0.22'+IRATIO+LABEL,par))
    Result(slo,'window      | transp |'
           +igrey('pclip=99 allpos=y bias=0.22'+IRATIO+LABEL,par))

    vel = 'v' + i
    Flow(vel,slo,'window | transp | smooth rect1=51 | math output=1/input')

    for k in IMGlist:

        # image / CIGs
        SRimg = 'SRI' + i + k
        SRcig = 'SRC' + i + k
        Flow(SRimg,SRcig,'window n4=1 min4=0 | transp')
        Plot  (SRimg,igrey('pclip=95'+IRATIO+LABEL,par))
        Result(SRimg,igrey('pclip=95'+IRATIO+LABEL,par))

        IMGSLO = 'IMGSLO' + i + k
        pplot.p2x1(IMGSLO,SRimg,slo,0.75,0.75,-6.5)

        # dip field
        SRdip = 'DIP' + i + k
        Flow(SRdip,SRimg,'dip rect1=3 rect2=3 order=3 liter=100')
        Result(SRdip,igrey('pclip=99 color=j'+IRATIO+LABEL,par))

        SRmap = 'MAP' + i + k
        Flow(SRmap,[SRdip,slo],'transp | remap1 pattern=${SOURCES[1]} | transp')

        # image after plane-wave distruction
        SRpwd = 'PWD' + i + k
        Flow(SRpwd,[SRimg,SRdip],'pwd dip=${SOURCES[1]}')
        Result(SRpwd,igrey('pclip=99'+IRATIO+LABEL,par))

        # dip-corrected velocity
        vco = 'v' + i + k
        Flow(vco,[SRmap,vel],
             '''
             math v=${SOURCES[1]} output="v/cos(atan(input*%(dz)g/%(dx)g))"
             ''' % par)

        for l in CIGlist:
            ikl = i + k + l

            # image
            z = 'z' + ikl
            lmin = float(l) - 20 * float('%(dx)g' % par)
            lmax = float(l) + 20 * float('%(dx)g' % par)
            Flow(z,SRcig,
                 '''
                 window n4=1 min4=0 min1=%g max1=%g |
                 transp |
                 bandpass fhi=10
                 ''' % (lmin,lmax) )
            Plot(z,z,igrey('labelsz=5 screenratio=6 screenht=10.5'+LABEL,par))

            # velocity
            v = 'v' + ikl
            Flow(v,vco,'window n2=1 min2=%g' % float(l) )

            # structural dip
            d = 'd' + ikl
            Flow(d,SRdip,'window n2=1 min2=%g' % float(l) )

            # offset gathers
            h = 'h' + ikl
            Flow(h,SRcig,
                 '''
                 window n1=1 min1=%g |
                 bandpass fhi=10
                 ''' % float(l))

            # angle gathers
            a    = 'a'    + ikl    # slant-stack
            c    = 'c'    + ikl    # angle-gather

            cdip = 'cdip' + ikl    # 
            ccln = 'ccln' + ikl    # dip noise separation

            if(k=='x'):
                Flow(a,h,'bandpass fhi=10 | '+hcig2ssk())
                Flow(c,a,                     hssk2ang())

                Plot(h,igrey('label2="h" unit2=km'                   +CRATIO+LABEL,par))
                Plot(a,igrey('label2="tan" unit2="\F10 q\F3 " min2=0'+CRATIO+LABEL,par))

            if(k=='t'):
                Flow(a,h,'bandpass fhi=10 | '+tcig2ssk())
                Flow(c,[a,v,d],               tssk2ang())

                Plot(h,igrey('label2="\F10 t\F3 " unit2=s '           +CRATIO+LABEL,par))
                Plot(a,igrey('label2="\F10 n\F3 " unit2="km/s" max2=5'+CRATIO+LABEL,par))

            # clean-up noise
            Flow(cdip,c,'twodip2 eps=10 lam=1 p0=-1 q0=0 both=n')
            Flow(ccln,[c,cdip],'pwdsigk dips=${SOURCES[1]} verb=y eps=1 | window n3=1 f3=1')

            Plot  (c   ,igrey('label2="\F10 q\F3 " unit2="\^o\_"'      +CRATIO+LABEL,par))
            Plot  (ccln,igrey('label2="\F10 q\F3 " unit2="\^o\_"'      +CRATIO+LABEL,par))

# < clean0t11.rsf sfwindow f3=1 > clean.rsf
# < clean.rsf sfenvelope | sfthreshold pclip=100 | sfmask min=0.1 | sfdd type=float | sfadd mode=p clean.rsf | sfstack | sfgraph | tube

# ------------------------------------------------------------

for i in SLOlist:
    for l in CIGlist:
        SRoff = 'SRoff' + i + '-' + l
        SRang = 'SRang' + i + '-' + l
        SRxic = 'SRxic' + i + '-' + l
        SRtic = 'SRtic' + i + '-' + l

        zx = 'z'    + i + 'x' + l
        zt = 'z'    + i + 't' + l

        hx = 'h'    + i + 'x' + l
        ht = 'h'    + i + 't' + l

        ax = 'a'    + i + 'x' + l
        at = 'a'    + i + 't' + l

        cx = 'ccln' + i + 'x' + l
        ct = 'ccln' + i + 't' + l

        p1x3(SRxic,hx,ax,cx,1,1,-1,-4,-7)
        p1x3(SRtic,ht,at,ct,1,1,-1,-4,-7)

        p1x3(SRoff,zx,hx,ht,1,1,-1,-3,-6)
        p1x3(SRang,zx,cx,ct,1,1,-1,-3,-6)

        SRt = 'SRt' +  i + '-' + l
        p1x4(SRt,zt,ht,at,ct,1,1,-1,-3,-6,-9)

        SRx = 'SRx' +  i + '-' + l
        p1x4(SRx,zx,hx,ax,cx,1,1,-1,-3,-6,-9)

#        SRic = 'SRic' + i + '-' + l
#        pplot.p2x1(SRic,SRxic,SRtic,0.5,0.5,-9)

#        SRcig = 'SRcig' + i + '-' + l
#        pplot.p2x1(SRcig,SRoff,SRang,0.5,0.5,-9)


p1x6('allxoff',   'h0x7',   'h0x9',   'h0x11',   'h0x13',   'h0x15',   'h0x17',0.75,0.75,-3)
p1x6('allxang','ccln0x7','ccln0x9','ccln0x11','ccln0x13','ccln0x15','ccln0x17',0.75,0.75,-3)

p1x6('alltoff',   'h0t7',   'h0t9',   'h0t11',   'h0t13',   'h0t15',   'h0t17',0.75,0.75,-3)
p1x6('alltang','ccln0t7','ccln0t9','ccln0t11','ccln0t13','ccln0t15','ccln0t17',0.75,0.75,-3)

pplot.p2x1('alloff','allxoff','alltoff',0.65,0.65,-7)
pplot.p2x1('allang','allxang','alltang',0.65,0.65,-7)

# ------------------------------------------------------------


End()