#!/bin/ksh

# Given a volume ($1) and a guess at the equilibrium c/a ratio
# for that volume ($2), calculates a new guess for the equilibrium
# c/a ratio using a quasi-Newton method.

# Usual disclaimers apply.

# Enter your location for the static executable here:

static=/usr/local/mehl/skfits/p6/static.105/static

# Get volume into skin file, using awk to provide a
#  standard format.

vol=`echo $1 | awk '{printf "%10.5f\n", $1}'`
# echo $vol

# Make sure the file skeng.$vol exists:

touch skeng.$1

sed "s/vol/$1/" skin.1 > skin.proto

# Center coa ratio:

coa2=`echo $2 | awk '{printf "%12.8f\n", $1}'`

# Use awk as a calculator.  "Eval" is better, if you have it:

coa1=`echo $coa2 | awk '{printf "%12.8f\n", 0.99*$1}'`
coa3=`echo $coa2 | awk '{printf "%12.8f\n", 1.01*$1}'`

# 'Step size'

hv=`echo $coa2 | awk '{printf "%12.8f\n", 0.01*$1}'`

# echo $coa1 $coa2 $coa3

sed "s/coa1/$coa1/" skin.proto | sed "s/coa2/$coa2/" | sed "s/coa3/$coa3/" > SKIN

$static > /dev/null

# Print the energies to the screen:

cat SKENG

# and to the file:

cat SKENG >> skeng.$1

# Now pick off the energies corresponding to coa[1-3]:

e1=`head -1 SKENG | awk '{print $4}'`
e2=`head -2 SKENG | tail -1 | awk '{print $4}'`
e3=`tail -1 SKENG | awk '{print $4}'`

# echo $hv $e1 $e2 $e3

# Find the first derivative.  I'm using awk because
# it's portable, not because it is pretty.  Any perl
# programmer who wants to rewrite all this is welcome.

ep=`echo $e1 $e3 $hv | awk '{printf "%12.8f\n", .5*($2-$1)/$3}'`

e2p=`echo $e1 $e2 $e3 $hv | awk '{printf "%12.8f\n", ($1+$3-2*$2)/($4*$4)}'`

# Now use these numbers to estimate the position of the equilibrium
#  and to get an estimate of the new equilibrium energy:

cnew=`echo $coa2 $ep $e2p | awk '{printf "%12.8f\n", $1-$2/$3}'`

enew=`echo $e2 $ep $e2p | awk '{printf "%15.9f\n", $1-.5*$2*$2/$3}'`

echo $vol $cnew $enew

coa2=$cnew

# Use awk as a calculator.  "Eval" is better, if you have it:

coa1=`echo $coa2 | awk '{printf "%12.8f\n", 0.99*$1}'`
coa3=`echo $coa2 | awk '{printf "%12.8f\n", 1.01*$1}'`

# 'Step size'

hv=`echo $coa2 | awk '{printf "%12.8f\n", 0.01*$1}'`

# echo $coa1 $coa2 $coa3

sed "s/coa1/$coa1/" skin.proto | sed "s/coa2/$coa2/" | sed "s/coa3/$coa3/" > SKIN

$static > /dev/null

# Print the energies to the screen:

cat SKENG

# and to the file:

cat SKENG >> skeng.$1

# Now pick off the energies corresponding to coa[1-3]:

e1=`head -1 SKENG | awk '{print $4}'`
e2=`head -2 SKENG | tail -1 | awk '{print $4}'`
e3=`tail -1 SKENG | awk '{print $4}'`

# echo $e1 $e2 $e3

# Find the first derivative.  I'm using awk because
# it's portable, not because it is pretty.  Any perl
# programmer who wants to rewrite all this is welcome.

ep=`echo $e1 $e3 $hv | awk '{printf "%12.8f\n", .5*($2-$1)/$3}'`

e2p=`echo $e1 $e2 $e3 $hv | awk '{printf "%12.8f\n", ($1+$3-2*$2)/($4*$4)}'`

# Now use these numbers to estimate the position of the equilibrium
#  and to get an estimate of the new equilibrium energy:

cnew=`echo $coa2 $ep $e2p | awk '{printf "%12.8f\n", $1-$2/$3}'`

enew=`echo $e2 $ep $e2p | awk '{printf "%15.9f\n", $1-.5*$2*$2/$3}'`

echo $vol $cnew $enew

coa2=$cnew

# Use awk as a calculator.  "Eval" is better, if you have it:

coa1=`echo $coa2 | awk '{printf "%12.8f\n", 0.99*$1}'`
coa3=`echo $coa2 | awk '{printf "%12.8f\n", 1.01*$1}'`

# 'Step size'

hv=`echo $coa2 | awk '{printf "%12.8f\n", 0.01*$1}'`

# echo $coa1 $coa2 $coa3

sed "s/coa1/$coa1/" skin.proto | sed "s/coa2/$coa2/" | sed "s/coa3/$coa3/" > SKIN

$static > /dev/null

# Print the energies to the screen:

cat SKENG

# and to the file:

cat SKENG >> skeng.$1

# Now pick off the energies corresponding to coa[1-3]:

e1=`head -1 SKENG | awk '{print $4}'`
e2=`head -2 SKENG | tail -1 | awk '{print $4}'`
e3=`tail -1 SKENG | awk '{print $4}'`

# echo $e1 $e2 $e3

# Find the first derivative.  I'm using awk because
# it's portable, not because it is pretty.  Any perl
# programmer who wants to rewrite all this is welcome.

ep=`echo $e1 $e3 $hv | awk '{printf "%12.8f\n", .5*($2-$1)/$3}'`

e2p=`echo $e1 $e2 $e3 $hv | awk '{printf "%12.8f\n", ($1+$3-2*$2)/($4*$4)}'`

# Now use these numbers to estimate the position of the equilibrium
#  and to get an estimate of the new equilibrium energy:

cnew=`echo $coa2 $ep $e2p | awk '{printf "%12.8f\n", $1-$2/$3}'`

enew=`echo $e2 $ep $e2p | awk '{printf "%15.9f\n", $1-.5*$2*$2/$3}'`

echo $vol $cnew $enew



# Clean up after ourselves:

rm SKIN skin.proto SKENG SKOUT