# Low frequency scalar wave scattering attenuation and penetration in 3-D isotropic von Karman medium
#
# December 2007
#
# Thomas Jules Browaeys 
# Bureau of Economic Geology
# University of Texas at Austin
# mailto:jules.browaeys@beg.utexas.edu


from rsf.proj import *
from math import pi


# ----------------------------------------------
# FREQUENCY 1-D GRID 
# ----------------------------------------------
# nx = points number
# dx = frequency steps

pgrid = {'nx':300, 'ox':1., 'dx':1.}

# Create frequency grid
Flow('fgrd',None,'spike nsp=1 mag=1 n1=%(nx)d d1=%(dx)g o1=%(ox)g | put label1=f unit1=Hz' % pgrid)

# Logarithmic scale
Flow('lfgrd','fgrd','math output="log(x1)/log(10.)"')


# ----------------------------------------------
# STOCHASTIC MEDIUM
# ----------------------------------------------
# hu = H exponent > -0.5 
# cu = H spectrum constant

# Exponents and constants
pgrids = {'M025':{'hu':-0.25, 'cu':1.311,  'shu':'\F15 8\F2 0.25'},
           '025':{'hu':0.25,  'cu':0.5991, 'shu':'0.25'},
            '05':{'hu':0.5,   'cu':1.0,    'shu':'0.5'},
           '075':{'hu':0.75,  'cu':1.3317, 'shu':'0.75'},
             '1':{'hu':1.0,   'cu':1.57,   'shu':'1.0'}
          }

# Correlation lengths (m)
bvalues = {'b25':2.5, 'b05':5.0, 'b10':10.0}
bvalstr = {'b25':'2.5', 'b05':'5', 'b10':'10'}

# VP VS variations and velocity
parq = {
    'sp':0.3,    # P wave standard deviation
    'ss':0.3,    # S wave standard deviation
    'cp':2700.,  # P wave velocity
    'cs':1230.   # S wave velocity
    }

# Variances
parq['vrs2'] = parq['ss']*parq['ss']
parq['vrp2'] = parq['sp']*parq['sp']


for lb in bvalues.keys():

    parq['b'] = bvalues[lb]
    bstring = bvalstr[lb]

    ksb  = 'ksb'   + lb
    ksb2 = 'ksb2'  + lb
    kpb  = 'kpb'   + lb
    kpb2 = 'kpb2'  + lb

    # Wave vector k*b for VS
    Flow(ksb,'fgrd','math output="(%g)*x1/(%g)"' % (2.*pi*parq['b'],parq['cs']))
    # Squared (k*b) for VS
    Flow(ksb2,ksb,'math output="input*input"')

    # Wave vector k*b for VP
    Flow(kpb,'fgrd','math output="(%g)*x1/(%g)"' % (2.*pi*parq['b'],parq['cp']))
    # Squared (k*b) for VP
    Flow(kpb2,kpb,'math output="input*input"')

    # -----------------------
    # Exponential correlation
    # -----------------------

    # S waves
    dvse    = 'dvse'     + lb
    qse     = 'qse'      + lb
    depthse = 'depthse'  + lb

    # Dispersion velocity
    Flow(dvse,ksb2,'math k2=${SOURCES[0]} output="1./(1.+(%g)*(0.5+4.*k2)/(1.+4.*k2))"' % (parq['vrs2']))
    # Attenuation 1/Q
    Flow(qse,[ksb2, ksb],'math k2=${SOURCES[0]} k=${SOURCES[1]} output="(%g)*(k2*k)/(1.+4.*k2)"' % (8.*parq['vrs2']))
    # Depth penetration
    Flow(depthse,ksb2,'math k2=${SOURCES[0]} output="(%g)/(%g)*(1.+4.*k2)/(k2*k2)"' % (parq['b'],4.*parq['vrs2']))

    # P waves
    dvpe    = 'dvpe'     + lb
    qpe     = 'qpe'      + lb
    depthpe = 'depthpe'  + lb

    # Dispersion velocity
    Flow(dvpe,kpb2,'math k2=${SOURCES[0]} output="1./(1.+(%g)*(0.5+4.*k2)/(1.+4.*k2))"' % (parq['vrp2']))
    # Attenuation 1/Q
    Flow(qpe,[kpb2, kpb],'math k2=${SOURCES[0]} k=${SOURCES[1]} output="(%g)*(k2*k)/(1.+4.*k2)"' % (8.*parq['vrp2']))
    # Depth penetration
    Flow(depthpe,kpb2,'math k2=${SOURCES[0]} output="(%g)/(%g)*(1.+4.*k2)/(k2*k2)"' % (parq['b'],4.*parq['vrp2']))

    lfdvse = 'lfdvse'   + lb
    lfdvpe = 'lfdvpe'   + lb

    lqse   = 'lqse'  + lb
    lqpe   = 'lqpe'  + lb
    lfqse  = 'lfqse' + lb
    lfqpe  = 'lfqpe' + lb

    lldve   = 'lldve'   + lb
    llqe    = 'llqe'    + lb
    depthle = 'depthle' + lb

    # Plots
    Flow(lfdvse,['lfgrd', dvse],'cmplx ${SOURCES[0:2]}',stdin=0)
    Flow(lfdvpe,['lfgrd', dvpe],'cmplx ${SOURCES[0:2]}',stdin=0)

    Flow(lqse,qse,'math output="log(input)/log(10.)"')
    Flow(lqpe,qpe,'math output="log(input)/log(10.)"')

    Flow(lfqse,['lfgrd', lqse],'cmplx ${SOURCES[0:2]}',stdin=0)
    Flow(lfqpe,['lfgrd', lqpe],'cmplx ${SOURCES[0:2]}',stdin=0)

    # Plot velocity dispersion
    Result(lldve,[lfdvpe, lfdvse],
           '''
           cat ${SOURCES[0:2]} axis=2 |
           put label1='log[f]' unit1='Hz' label2='c/c\_0\^' unit2= |
           graph dash=0,3 title="Dispersion velocity - H=0.5" min1=0. max1=2.5 min2=0.9 max2=1.
           ''',stdin=0)

    # Plot attenuation
    Result(llqe,[lfqpe, lfqse],
           '''
           cat ${SOURCES[0:2]} axis=2 |
           put label1='log[f]' unit1='Hz' label2='log[1/Q]' unit2= |
           graph dash=0,3 title="Attenuation - H=0.5" min1=0. max1=2.5 max2=0. min2=-6.
           ''',stdin=0)

    # Plot depth penetration
    Result(depthle,[depthpe, depthse],
           '''
           cat ${SOURCES[0:2]} axis=2 |
           put label2='d' unit2='m' |
           graph dash=0,3 title="Penetration depth - H=0.5" min1=0. max1=100. min2=0. max2=2500.
           ''',stdin=0)

    # -------------------
    # Fractal correlation
    # -------------------

    for fexp in ('M025','025','075'):

        parq.update(pgrids[fexp])

        # S waves
        hksb4   = 'hksb4'   + lb + fexp
        qsf     = 'qsf'     + lb + fexp
        depthsf = 'depthsf' + lb + fexp

        Flow(hksb4,[ksb2, ksb],'math k2=${SOURCES[0]} k=${SOURCES[1]} output="k-k/(1.+4.*k2)^(%g)"' % (0.5+parq['hu']))
        # Attenuation 1/Q
        Flow(qsf,hksb4,'math hk4=${SOURCES[0]} output="2.*(%g)*hk4"' % (parq['vrs2']*parq['cu']))
        # Depth penetration
        Flow(depthsf,[hksb4, ksb],'math hk4=${SOURCES[0]} k=${SOURCES[1]} output="(%g)/(%g)/(k*hk4)"' % (parq['b'],parq['vrs2']*parq['cu']))


        # P waves
        hkpb4   = 'hkpb4'   + lb + fexp
        qpf     = 'qpf'     + lb + fexp
        depthpf = 'depthpf' + lb + fexp

        Flow(hkpb4,[kpb2, kpb],'math k2=${SOURCES[0]} k=${SOURCES[1]} output="k-k/(1.+4.*k2)^(%g)"' % (0.5+parq['hu']))
        # Attenuation 1/Q
        Flow(qpf,hkpb4,'math hk4=${SOURCES[0]} output="2.*(%g)*hk4"' % (parq['vrp2']*parq['cu']))
        # Depth penetration
        Flow(depthpf,[hkpb4, kpb],'math hk4=${SOURCES[0]} k=${SOURCES[1]} output="(%g)/(%g)/(k*hk4)"' % (parq['b'],parq['vrp2']*parq['cu']))

        lqsf   = 'lqsf' + lb + fexp
        lqpf   = 'lqpf' + lb + fexp

        lfqsf  = 'lfqsf' + lb + fexp
        lfqpf  = 'lfqpf' + lb + fexp

        llqf    = 'llqf'    + lb + fexp
        depthlf = 'depthlf' + lb + fexp

        # Plots
        Flow(lqsf,qsf,'math output="log(input)/log(10.)"')
        Flow(lqpf,qpf,'math output="log(input)/log(10.)"')

        Flow(lfqsf,['lfgrd', lqsf],'cmplx ${SOURCES[0:2]}',stdin=0)
        Flow(lfqpf,['lfgrd', lqpf],'cmplx ${SOURCES[0:2]}',stdin=0)

        # Plot attenuation
        titleatt = 'title="Attenuation \k100 \k100 \F5 b\F2 \k70 =\k70 ' + bstring + '\k10 m\g \G \F5 H\F2 \k60 =\k60 ' + parq['shu'] + '"'
        Result(llqf,[lfqpf, lfqsf],
               '''
               cat ${SOURCES[0:2]} axis=2 |
               put label1='log[f]' unit1='Hz' label2='log[1/Q]' unit2= |
               graph dash=0,3 min1=0. max1=2.5 max2=0. min2=-6. font=2
               parallel2=n plotfat=4,4 n2tic=7 o2num=-6.0 d2num=1.0 %s
               ''' % (titleatt),stdin=0)

        # Plot depth penetration
        titledep = 'title="Penetration depth \k100 \k100 \F5 b\F2 \k70 =\k70 ' + bstring + '\k10 m\g \G \F5 H\F2 \k60 =\k60 ' + parq['shu'] + '"'
        Result(depthlf,[depthpf, depthsf],
               '''
               cat ${SOURCES[0:2]} axis=2 |
               put label2='d' unit2='m' |
               graph dash=0,3 min1=0. max1=100. min2=0. max2=2500. font=2
               parallel2=n plotfat=4,4 n2tic=6 o2num=0.0 d2num=500.0 %s
               ''' % (titledep),stdin=0)



End()