# Spatial energy spectrum for 2-D fluctuations (data and synthetic fGn)
#
# July 2008
#
# Thomas Jules Browaeys 
# Bureau of Economic Geology
# University of Texas at Austin
# mailto:jules.browaeys@beg.utexas.edu


from rsf.proj import *
from rsf.recipes.beg import server as private
from math import pi

# ----------------------------------------------------------------------
# RESERVOIR MODEL VS DATA 2-D CROSSLINE SECTION
# ----------------------------------------------------------------------

# Image bias
pcb = {}

# Grid
pgrid = {'nz':609, 'oz':4., 'dz':2., 'nx':521, 'ox':272., 'dx':17.}

# VS model
dfile = 'coronation_xl590_vs.HH'
Fetch(dfile,'apache',private)

# Plot crossline VS model  
Flow('dfile',dfile,'dd form=native | math output="input" | put label1=z label2=x unit1=m unit2=m')
Result('dfile','grey color=a title="VS velocity variations" parallel2=n')


# ----------------------------------------------------------------------
# SCALE AND WINDOW VARIATIONS
# ----------------------------------------------------------------------

# Window
pgrid.update({'ozw':200., 'nzw':351, 'oxw':272., 'nxw':521})

# Bias for scaled image
pcb['bvv'] = 0.0

# Window 
Flow('wdfile','dfile','window min1=%(ozw)g n1=%(nzw)d min2=%(oxw)g n2=%(nxw)d' % pgrid)

# Shift for zero mean scaled signal
Flow('mscale','wdfile','stack axis=2 norm=n | stack axis=1 norm=n | math output="input*%g"' % (1./(pgrid['nzw']*pgrid['nxw'])))
Flow('mdfile','mscale',
     '''
     window | spray axis=1 n=%(nzw)d o=%(ozw)g d=%(dz)g |
     spray axis=2 n=%(nxw)d o=%(oxw)g d=%(dx)g | put label1=z
     ''' % pgrid)
Flow('signal',['wdfile','mdfile'],'math r=${SOURCES[0]} m=${SOURCES[1]} output="r/m-1."')

# Plot
Result('signal',
       '''
       grey color=b bias=%(bvv)g 
       title="\F5 V\_\s70 \F2 S\s100 \^ velocity variations" font=2
       parallel2=n n2tic=8 o2num=200 d2num=100
       label1="\F2 z" label2="\F2 x"
       unit1="\F2 m" unit2="\F2 m"
       ''' % pcb)


# ----------------------------------------------------------------------
# SPATIAL FREQUENCY SPECTRUM
# ----------------------------------------------------------------------

# Parameters for clipping minimum frequency fz (1/m)
pcb.update({'fzc':0.0025})

# Fourier transform
fft2 = 'fft1 sym=y | fft3 sym=y'
Flow('fdfile','signal',fft2)

# Spatial frequency energy spectrum
Flow('f2dfile','fdfile','math output="input*conj(input)" | real | window min1=%(fzc)g' % pcb)

# Plot energy spectrum
# Result('f2dfile',
#       '''
#       put label1=fz label2=fx unit1=1/m unit2=1/m |
#       grey allpos=y title="Energy density spectrum" parallel2=n n2tic=5 o2num=0.05 d2num=0.05
#       ''')

Result('f2dfile',
       '''
       put label1=fz label2=fx unit1=1/m unit2=1/m |
       grey color=H allpos=y
       parallel2=n n2tic=5 o2num=0.05 d2num=0.05 format2=%3.2f
       title="\F2 Energy density spectrum" font=2
       label1="\F2 fz" label2="\F2 fx"
       unit1="\F2 m\^\s50 -1\s100 \_" unit2="\F2 m\^\s50 -1\s100 \_"
       ''')

# ----------------------------------------------------------------------
# SPATIAL FREQUENCY 2-D GRID AND RANDOM MEDIUM
# ----------------------------------------------------------------------


# Spatial grid (m)
pgrids = {'nz':701,  'oz':200., 'dz':1.,'nx':8841, 'ox':272., 'dx':1.}

# Create spatial grid
Flow('spacegrid',None,'spike mag=1 n1=%(nz)d n2=%(nx)d d1=%(dz)g d2=%(dx)g o1=%(oz)g o2=%(ox)g label1="z" label2="x" unit1="m" unit2="m"' % pgrids)

# Medium heterogeneities 
pmd = {
    'bx':2500.,       # correlation length in x (m)
    'bz':10.,         # correlation length in z (m)
    'nu':-0.5,        # exponent H
    'wst':0.1         # standard deviation
    }

# Gaussian white noise
pwhite = {
    'wmu':0.,                     # mean
    'wvr':pmd['wst']*pmd['wst'],  # variance
    'wrn':30.,                    # noise range
    'wsd':7                       # seed for random generator
    }

# Generate white noise
Flow('wgauss','spacegrid','noise mean=%(wmu)g range=%(wrn)g seed=%(wsd)g var=%(wvr)g rep=y type=y' % pwhite)

# Frequency window clipping
fcb = {
    'fzc':0.0025,       # min frequency fz (1/m)
    'fzm':0.25,         # max frequency fz (1/m)
    'fxc':-0.0294,      # min frequency fx (1/m)
    'fxm':0.0294,       # max frequency fx (1/m)
    'bcg':0.0           # bias for final plot
    }

# Plot white noise 
Result('wgauss','window j2=4 | grey color=j parallel2=n bias=%g title="White noise"' % (fcb['bcg']))

# Fourier transform of white noise
Flow('fwgauss','wgauss',fft2)

# Zero frequency amplitude
specm = pmd['nu']*4.*pi*pmd['bx']*pmd['bz']

# Unitary von Karman spectrum 
Flow('specsde','fwgauss','math output="(1.+(x1*2.*(%g))^2+(x2*2.*(%g))^2)^(%g)"' % (pi*pmd['bz'],pi*pmd['bx'],-pmd['nu']-1.))

# Fractional Gaussian noise (fGn)
Flow('fcgaussb',['specsde','fwgauss'],'math r=${SOURCES[0]} p=${SOURCES[1]} type=complex output="sqrt(abs(%g)*r)*p"'% (specm))

# Inverse Fourier transform for fGn (von Karman medium)
ifft2 = 'sffft3 sym=y inv=y | sffft1 sym=y inv=y'
Flow('cgaussb','fcgaussb',ifft2)

# Plot correlated spatial random medium spectrum

Flow('fcgaussbp','fcgaussb',
     '''
     math output="input*conj(input)" | real | window min1=%(fzc)g max1=%(fzm)g min2=%(fxc)g max2=%(fxm)g |
     put label1=fz label2=fx unit1=1/m unit2=1/m
     ''' % fcb)

Result('fcgaussbp',
       '''
       grey color=H allpos=y
       title="\F2 Energy density spectrum" font=2
       parallel2=n n2tic=5 o2num=0.05 d2num=0.05
       label1="\F2 fz" label2="\F2 fx" format2=%3.2f
       unit1="\F2 m\^\s50 -1\s100 \_" unit2="\F2 m\^\s50 -1\s100 \_"
       ''')

# Plot correlated spatial random medium
Result('cgaussb',
       '''
       window j2=4 |
       grey color=b bias=%g 
       parallel2=n n2tic=8 o2num=200 d2num=100
       title="\F5 b\_x\^\F2 \k70 =\k70 2500\k10 m\g \G \F5 b\_z\^\F2 \k60 =\k60 10\k10 m\g \G \F5 H\F2 \k60 =\k60 \F15 8\F2 0.5" font=2
       label1="\F2 z" label2="\F2 x"
       unit1="\F2 m" unit2="\F2 m"
       ''' % (fcb['bcg']))


End()