from rsfproj import *
import rwezo,zomig,fdmod

par = {
    'nz':1201, 'dz':0.025,  'oz':0,     'lz':'label1=z', 'uz':'km',
    'nx':2133, 'dx':0.0375, 'ox':10.925,'lx':'label2=x', 'ux':'km',
    'nt':30000,'dt':0.0005, 'ot':0,     'lt':'label1="\F10 t\F3 " unit1="s"', 'ut':'s',
    'ng':800,  'dg':1,      'og':0,     'lg':'label2="\F10 g\F3 " unit2="."',
    'kt':100,    # wavelet delay 
    'nT':1600, 'dT':0.004,  'oT':0,
    'ow':1,
    'nbz':250,'nbx':250
    }
par['ox']*=0.3048
par['dx']*=0.3048
par['oz']*=0.3048
par['dz']*=0.3048

# compute parameters
rwezo.param(par)
zomig.migpar(par)
fdmod.param(par)

# ------------------------------------------------------------
# velocity
velo = 'sigsbee2a_migvel.sgy'
Fetch(velo,'sigsbee')

Flow('zvelo tzvelo',velo,
     '''
     segyread
     tape=$SOURCE
     tfile=${TARGETS[1]}
     hfile=${TARGETS[2]}
     bfile=${TARGETS[3]}
     ''',stdin=0)

# ------------------------------------------------------------
par['nxcut']=1600
par['nxpad']=2000
par['nzpad']=799

Flow('vraw','zvelo',
     '''
     scale rscale=0.001 |
     scale rscale=0.3048 |
     put d1=%(dz)g o2=%(ox)g d2=%(dx)g label1=z label2=x |
     window n2=%(nxcut)d
     ''' % par)

par['ox']=par['ox']-par['nxpad']*par['dx']
par['nx']=par['nxcut']+par['nxpad']
par['nz']=par['nz']+par['nzpad']

par['xmin']=par['ox']
par['xmax']=par['ox'] + (par['nx']-1) * par['dx']
par['zmin']=par['oz']
par['zmax']=par['oz'] + (par['nz']-1) * par['dz']

Flow('padx','vraw',
     '''
     window n2=1 |
     spray axis=2 n=%(nxpad)d d=%(dx)g o=%(ox)g |
     cat ${SOURCES[0]} axis=2 space=n
     ''' % par)

Flow('padz','padx',
     '''
     window n1=1 f1=1200 |
     spray axis=1 n=%(nzpad)d d=%(dz)g o=0 |
     math "output=input+x1*0.25"
     ''' % par)

Flow('velo','padx padz',
     '''
     cat axis=1 space=n ${SOURCES[1]}
     ''')
Flow('dens','velo','math output=2.5')

Plot  ('velo',rwezo.cgrey('allpos=y pclip=95',par))
Result('velo',rwezo.cgrey('allpos=y pclip=95',par))

Flow(  'salt','velo','add add=-4.26 | clip clip=0.3048')
#Result('salt',rwezo.cgrey('allpos=y pclip=100',par))

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

Flow(  'edgez','salt','         igrad square=y')
Flow(  'edgex','salt','transp | igrad square=y | transp')
Flow(  'edge','edgez edgex',
       'math z=${SOURCES[0]} x=${SOURCES[1]} output="z+x"')
Result('edge','sfdd type=float |' + rwezo.cgrey('pclip=100',par))

Flow(  'mask','edge','mask min=0.0006')
Result('mask','sfdd type=float |' + rwezo.cgrey('pclip=100',par))

for x in ('x1','x2'):
    if(x=='x1'): o='zs'
    if(x=='x2'): o='xs'

    Flow(o,'edge mask',
         '''
         math output=%s |
         put n1=1 n2=%d |
         headerwindow mask=${SOURCES[1]} |
         window
         ''' % (x,par['nz']*par['nx']))

Flow('rs','edge mask',
     '''
     scale axis=123 |
     put n1=1 n2=%d |
     headerwindow mask=${SOURCES[1]} |
     window
     ''' % (par['nz']*par['nx']))

Flow('ss',['xs','zs','rs'],
     '''
     cat axis=2 space=n
     ${SOURCES[0]} ${SOURCES[1]} ${SOURCES[2]} | transp |
     window j2=2
     ''', stdin=0)

Plot('ss','window n1=2 | dd type=complex | window | '
     + fdmod.cgraph('symbol=. plotcol=0 symbolsz=8',par))
Result('ss',['velo','ss'],'Overlay')

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

# ------------------------------------------------------------
# receiver (flat)

Flow('ro',None,'math n1=%(nx)d d1=%(dx)g o1=%(ox)g output=0' % par)

Flow('zo','ro','math output="%g"' % par['oz'])
Flow('xo','ro','math output="x1"')
Flow('oo',['xo','zo'],
     '''
     cat axis=2 space=n
     ${SOURCES[0]} ${SOURCES[1]} | transp
     ''', stdin=0)

Plot('oo','window n1=2 | dd type=complex | window j2=10 | '
     + fdmod.cgraph('symbol=* plotcol=5',par))
#Result('oo',['velo','oo'],'Overlay')

Result('model','velo ss oo','Overlay')

Flow('surf',None,
     '''
     spike nsp=1 mag=1 n1=%(nz)d d1=%(dz)g o1=%(oz)g k1=2 |
     spray axis=2       n=%(nx)d  d=%(dx)g  o=%(ox)g
     ''' % par)

# ------------------------------------------------------------
# wavelet
Flow('wav',None,
     '''
     spike nsp=1 mag=1 n1=%(nt)d d1=%(dt)g o1=%(ot)g k1=%(kt)d |
     ricker1 frequency=20 | 
     put label1=t label2=x label3=y 
     ''' % par)
Result('wav','wav','window n1=1000 | graph title=" "')

# ------------------------------------------------------------
# modeling
fdmod.amodel('do','wo','wav','velo','dens','ss','oo','free=n',par)

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

Result('wo',
       'window j3=20 min1=%(oz)g n1=%(nz)d min2=%(ox)g n2=%(nx)d |' % par
       + rwezo.cgrey('pclip=99 gainpanel=a',par))
Result('do','window j2=4 | transp |' + rwezo.dgrey('',par))

for j in range(20,300,20):
    par['j'] = j

    Plot('wo'+str(j/20),'wo',
         '''
         window n3=1 f3=%(j)d
         min1=%(oz)g n1=%(nz)d min2=%(ox)g n2=%(nx)d |
         ''' % par
         + rwezo.cgrey('pclip=99',par))
    Result('wo'+str(j/20),['wo'+str(j/20),'cos1','ss'],'Overlay')

Flow(  'dd','do','window f2=%(kt)d | put o2=0 d2=0.001 | window j2=4' % par)
Result('dd','transp |' + rwezo.dgrey('',par))

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

# ------------------------------------------------------------
Flow('vel1','velo','smooth rect1=351 rect2=351 repeat=3')

Flow('s_',None,'spike nsp=1 n1=%(ng)d d1=%(dg)g o1=%(og)g' % par)

# coordinate system: rays(g,t)
Flow('mzs1','s_','math output="01.0+0.009*x1"' % par)
Flow('mxs1','s_','math output="17.0+0.008*x1"' % par)
Flow('sou1','mxs1 mzs1',
     'cat axis=2 space=n ${SOURCES[1]} | transp | dd type=complex | window')

# wavefronts (by HWT)
Flow('hwt1','vel1 sou1',
     '''
     hwtex verb=n sou=${SOURCES[1]}
     nt=%(nt)d ot=%(ot)g dt=%(dt)g
     ''' % par)
Flow('cos1','hwt1',
     '''
     reverse which=2 |
     put o2=0 d2=%(dt)g |
     window j1=2 j2=8
     ''' % par)

# ------------------------------------------------------------
Flow('vel2','velo','math output=2.75')

# coordinate system: rays(g,t)
Flow('mzs2','s_','math output="-06.5+0.025*x1"' % par)
Flow('mxs2','s_','math output="-17.0-0.005*x1"' % par)
Flow('sou2','mxs2 mzs2',
     'cat axis=2 space=n ${SOURCES[1]} | transp | dd type=complex | window')

# wavefronts (by HWT)
Flow('hwt2','vel2 sou2',
     '''
     hwtex verb=n sou=${SOURCES[1]}
     nt=%(nt)d ot=%(ot)g dt=%(dt)g
     ''' % par)
Flow('cos2','hwt2','window j2=4')

# ------------------------------------------------------------
# time data CC: (t,x)
Flow  ('datCC','dd','transp')
Result('datCC','grey pclip=99 wanttitle=n label1=t unit1=s label2=x unit2=km')

# ------------------------------------------------------------
# WEM image CC: (z,x)
zomig.image('imgCC','sloCC','frqCC',par)
Plot  ('imgCC','window | transp |' + rwezo.cgrey('pclip=99',par))
Result('imgCC',['imgCC','cos1'],'Overlay')

# loop over coordinate systems
# 1 - Riemannian coordinate system
# 2 - Cartesian  coordinate system
for i in ('1','2'):
    # coordinate system plot
    rwezo.cos(   'cos'+i,20,100,par)
    Result('cos'+i,['velo','cos'+i],'Overlay')

    # slowness
    rwezo.slow(  'sloCC','sloRC'+i,'velo','cos'+i,par)
    Result('sloRC'+i,'transp |'
           + rwezo.rgrey('pclip=99 allpos=y color=j bias=0.061',par))

    # surface and salt
    Result('edgeRC'+i,['edge','surf','cos'+i],
           '''
           add ${SOURCES[1]} | transp |
           c2r rays=${SOURCES[2]} adj=n | transp | bandpass fhi=10 |
           scale axis=123 |
           '''+ rwezo.rgrey('pclip=99 scalebar=y',par))

    # ------------------------------------------------------------
    # coordinate system coefficients (A,B,M=mask) and references (A,B)
    rwezo.abm(    'abmRC'+i,'abrRC'+i,'sloRC'+i,'cos'+i,par)
    rwezo.abmplot('abmRC'+i,par)

    # ------------------------------------------------------------
    # frequency data RC: (g,w)
    rwezo.frq('frqRC'+i,'frqCC','datCC','cos'+i,par)

    # ------------------------------------------------------------
    # RWE image
    rwezo.mig('migCC'+i,'migRC'+i,'frqRC'+i,'abmRC'+i,'abrRC'+i,'cos'+i,par)

# ------------------------------------------------------------
# PLOTS
#plots(par)

# ------------------------------------------------------------
zoompar = {
    'nz':7,    'dz':1,      'oz':2,     'lz':'label1=z unit1=km',
    'nx':10,   'dx':1,      'ox':12,    'lx':'label2=x unit2=km',
    'nt':30000,'dt':0.0005, 'ot':0,     'lt':'label1="\F10 t\F3 " unit1="s"',
    'ng':400,  'dg':1,      'og':0,     'lg':'label2="\F10 g\F3 " unit2="."',
    'nT':1600, 'dT':0.004,  'oT':0
    }
rwezo.param(zoompar)

Plot('zoomvel','velo',rwezo.cgrey('',zoompar))
Flow('vscaleCC','sloCC',
     'window | transp | math output=1/input | scale axis=123 | add add=-0.6')

for i in ('1','2'):
    Flow('vscaleRC'+i,'sloRC'+i,
         'transp | math output=1/input | scale axis=123 | add add=-0.6')

    Flow(     'zoomcos'+i,'cos'+i,'window')
    rwezo.cos('zoomcos'+i,5,25,zoompar)
    Result(   'zoomcos'+i,  ['zoomvel','zoomcos'+i],'Overlay')
    Plot('zoomvel-ovl'+i,['zoomvel','zoomcos'+i],'Overlay')

    for j in (['SSF','FFD','PSC',
               'F15','F45','F60']):
        sfx = i + '-' + j

        Flow('iscaleCC'+sfx,'migCC'+sfx,'transp | scale axis=123')
        Flow('compCC'+sfx,['vscaleCC','iscaleCC'+sfx],
             'math a=${SOURCES[0]} b=${SOURCES[1]} output="0.5*a+b"')
        Plot(  'compCC'+sfx,rwezo.cgrey('pclip=100',zoompar))
        Result('compCC'+sfx,['compCC'+sfx,'zoomcos'+i],'Overlay')

        Flow('iscaleRC'+sfx,'migRC'+sfx,'transp | scale axis=123')
        Flow('compRC'+sfx,['vscaleRC'+i,'iscaleRC'+sfx],
             'math a=${SOURCES[0]} b=${SOURCES[1]} output="0.5*a+b"')
        Result('compRC'+sfx,rwezo.rgrey('pclip=100 min1=10',zoompar))

End()