# Comparison of 1D Fourier spectra from sonic log data and synthetic fBm
#
# 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 math import pi
from rsf.recipes.beg import server as private

# ----------------------------------------------------------------------
# WELL LOG DATA 
# ----------------------------------------------------------------------

# WELLS C
# DATA (04): z(col=0), rho(col=1), VP(col=2), VS(col=3)
dfile = '04-02-37-10.txt'

Fetch(dfile,'apache',private)

# Grid 1-D 
pgrid = {'nx':6921, 'ny':4, 'ox':142.5, 'dx':0.125}

# Frequency restriction and plot parameters
# Well C - VS
optrfs = 'title="Sonic log \F5 V\_\s70 \F2 S\s100 \^ spectrum \k100 \k100 Well 3" n2tic=6 o2num=0.0 d2num=0.5 format2=%3.1f'
soniclog = 'C2'
pgrid['col'] = 3

# Files
sonic    = 'sonic'    + soniclog
mscale   = 'mscale'   + soniclog
msonic   = 'msonic'   + soniclog
signal   = 'signal'   + soniclog
fsignal  = 'fsignal'  + soniclog
rfsignal = 'rfsignal' + soniclog

# Sonic log signal
Flow(sonic,dfile,
     '''
     echo n1=%(ny)d n2=%(nx)d o2=%(ox)g d2=%(dx)g in=$SOURCE data_format=ascii_float |
     dd form=native | transp |
     window n1=%(nx)d n2=1 f2=%(col)d |
     put label1=Depth unit1=m label2=Velocity unit2=m/s
     ''' % pgrid)

# Plot
Result(sonic,'graph title="Sonic well log"')

# Shift and scale by mean signal
Flow(mscale,sonic,'stack axis=1 norm=n | math output="input*%g"' % (1./pgrid['nx']))
Flow(msonic,mscale,'spray axis=1 n=%d o=%g' % (pgrid['nx'],pgrid['ox']))
Flow(signal,[sonic,msonic],'math r=${SOURCES[0]} m=${SOURCES[1]} output="r/m-1."')

# Plot scaled signal
Result(signal,
       '''
       graph min1=%g max1=%g min2=%g max2=%g label2= unit2= %s parallel2=n
       ''' % (pgrid['ox'],pgrid['ox']+(pgrid['nx']-1.)*pgrid['dx'],-1.4,1.4,'title="Sonic log V\_S\^ - Well 3" format2=%3.1f  n2tic=5 o2num=-1.0 d2num=0.5' ))

# Discrete Fourier Transform (fft) (k=0,N-1)
Flow(fsignal,signal,'fft1 sym=y')
Flow(rfsignal,fsignal,'add abs=y | real | put label1=f unit1=1/m label2=Amplitude unit2=m')

# Plot data spectrum 
Result(rfsignal,
       '''
       graph font=2 unit1="m\^\s50 -1\s100 \_"
       min1=0. min2=0. max1=1. max2=2.5 parallel2=n %s 
       ''' % (optrfs))


# ----------------------------------------------
# SYNTHETIC fBm spectrum
# ----------------------------------------------

# Spatial grid (m)
pgridf = {'nx':4056, 'ox':1., 'dx':0.25}

# Create spatial grid
Flow('spacegrid',None,'spike nsp=1 mag=1 n1=%(nx)d d1=%(dx)g o1=%(ox)g | put label1=z unit1=m' % pgridf)

# Signal parameters (Gaussian white noise mean and range)
psyn = {'wmu':0.0, 'wrn':1.0}

# Windowing maximum depth (m)
wdcp = {'xmax':1000.}

# Fractional Brownian motion fBm (nu > 0)
hnux = '025'
psyn.update({
    'b':5.,        # correlation length (m)
    'nu':0.25,     # H exponent
    'ca':0.5991,   # H spectrum constant
    'wdd':0.3,     # Gaussian white noise standard deviation
    'wsd':317      # random generator seed
    })

# Title
titlefit = 'title="Synthetic fBm spectrum \F5 b\F2 \k70 =\k70 5\k10 m\g \G \F5 H\F2 \k60 =\k60 0.5" n2tic=6 o2num=0.0 d2num=0.5 format2=%3.1f'

# Variance of Gaussian white noise
psyn['wvr'] = psyn['wdd']*psyn['wdd']

wgauss     =  'wgauss'   + hnux
fwgauss    = 'fwgauss'   + hnux
vkfilt     = 'vkfilt'    + hnux
fcgauss    = 'fcgauss'   + hnux
rfcgauss   = 'rfcgauss'  + hnux
irfcgauss  = 'irfcgauss' + hnux
fitfiltb   = 'fitfiltb'  + hnux

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

# Fourier transform of Gaussian noise
Flow(fwgauss,wgauss,'fft1 sym=y')

# Von Karman 1D filter in spectral domain
Flow(vkfilt,fwgauss,'math output="sqrt((%g))*(1.+((%g)*x1)^2)^(%g)"' % (2.*psyn['b']*psyn['ca'],2.*pi*psyn['b'],-0.25-0.5*psyn['nu']))

# Stochastic process von Karman 1D filtering and integration
Flow(fcgauss,[vkfilt,fwgauss],'math r=${SOURCES[0]} p=${SOURCES[1]} type=complex output="-I*r*p"')
Flow(rfcgauss,fcgauss,'add abs=y | real | put label1=f unit1=1/m')

# Separable least squares for exponent and amplitude
# Newton algorithm on nonlinear parameter b*b    
Flow(irfcgauss,rfcgauss,'karman verb=y niter=100 x0=1.')

# Plot synthetic spectrum and fit curve
Result(fitfiltb,[rfcgauss, irfcgauss],
           '''
           cat ${SOURCES[0:2]} axis=2 |
           put label1=f unit1=1/m label2=Amplitude unit2=m |
           graph min1=0. min2=0. max1=1. max2=2.5
           unit1="\F2 m\^\s50 -1\s100 \_" font=2
           parallel2=n plotfat=1,5 %s
           ''' % (titlefit),stdin=0)

End()