#
#  this script shows how to manipulate arrays from feff.dat
#  files and how to reconstruct chi(k) from this data.  
#  
#  it also compares the simple XAFS formula to that used by
#  ff2chi(), in which the full complex k-dependence of the
#  cumulants is included.
#

macro path_calcs 1 x1
  'explicitly generate full chi(k) for a path using simple XAFS formula'
  get_path($1, prefix= p1 ,do_arrays)

  # the arrays generated by this are:
  #  p1.k, p1.amp, p1.phase, p1.caps, p1.rep, p1.lambda
  #  are on the 'feff grid', and reflect the data from feff.dat
  #    amp  = [col 3] * [col 5]    : Scattering Amplitude
  #    caps = [col 2]              : Central Atom Phase Shift
  #    phase= [col 2] + [col 4]    : Total Phase Shift

  #
  # we'll generate chi(k) from these data
  # first, interpolate these onto the regular k grid

  set $2.k     = range(0,20,0.05)
  set $2.amp   = linterp(p1.k, p1.amp,    $2.k)
  set $2.phase = linterp(p1.k, p1.phase,  $2.k)
  set $2.caps  = linterp(p1.k, p1.caps,   $2.k)
  set $2.rep   = linterp(p1.k, p1.rep,    $2.k)
  set $2.lambda= linterp(p1.k, p1.lambda, $2.k)

  #
  # now construct XAFS according to XAFS formula:
  #   note that we carefully avoid 'divide-by-zero' at k=0  
  set $2.kr2   = (($2.k + 1.e-17 ) * p1_reff * p1_reff)
  set $2.ss2   = exp(-2*p1_sigma2*$2.k*$2.k)
  set $2.mag   = $2.amp * $2.ss2 * exp(-2*p1_reff /$2.lambda) / $2.kr2
  set $2.pha   = $2.phase + 2 * (p1_reff) * $2.k

  # using this correction would help some:
  # set $2.pha   = $2.phase + 2 * (p1_reff - 2 *p1_sigma2/p1_reff) * $2.k

  set $2.chi   = p1_degen * p1_S02 * $2.mag * sin($2.pha)
end macro

macro compare ss2  scale
   'compare ff2chi and simple XAFS Formula for a value of sigma2'

   pause == compare FF2CHI and simple XAFS Formula with sigma2=$1 ==

   set ss2 = $1
   path_calcs 1 x1
   ff2chi(1, prefix=f1)
   set  x1.diff = (x1.chi-f1.chi) * $2

   newplot( x1.k, (x1.chi-0.25),  key = 'XAFS formula', 
            title= 'ff2chi v. simple XAFS formula with \gs\u2\d = $1' )
   plot f1.k, (f1.chi-0.25),  key = 'FF2CHI'

   plot x1.k, x1.diff,    key = 'diff * $2', ymax=0.5,ymin=-1
end macro


#------------------------------------------------------
#
# step 1: define path
set ss2  = 0.0
path(1, ../feff/feff0001.dat, label ='Path 1', sigma2=ss2)

# now compare ff2chi v. simple XAFS formula as sigma2 increases:

compare 0.000    1.e14
compare 0.00001  1000
compare 0.001    10
compare 0.003    10
compare 0.020    10


#
# We could now fit the reconstructed data to the feff path...
# an exercise for the reader

# guess ss2 = 0.002
# 
# set kmin  = 1, kmax=18, dk=1, kweight=2, rmin=1, rmax=3
# 
# feffit(chi=x1.chi, 1, group=fit)
# 
# show @variables
