#!/bin/ksh

# Given a strain x ($1) and a guess at the equilibrium value of the
# internal parameter u ($2) for that value of x,
# calculates a new guess for the equilibrium
# value of x 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
static=/disk/dumpy/mehl/skfits/many/static.105/static

# Get x into skin file, using awk to provide a
#  standard format.  Note that we've cut down on
#  the size of x to fit the label into the proper field
#  along with u.

xlabel=`echo $1 | awk '{printf "%7.4f\n", $1}' | sed "s/ //g"`
# echo $xlabel

# the output file:

skfile=skeng.$xlabel
# echo $skfile

# make sure it exists:

touch $skfile

sed "s/xxxx/$xlabel/" skin.1 > skin.proto

# Central value of u.

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

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

u1=`echo $u2 | awk '{printf "%13.9f\n", $1-.001}'`
u3=`echo $u2 | awk '{printf "%13.9f\n", $1+.001}'`

# For the internal positions we need ui+1/3, and 2/3-ui

u1p=`echo $u1 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u2p=`echo $u2 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u3p=`echo $u3 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u1m=`echo $u1 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`
u2m=`echo $u2 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`
u3m=`echo $u3 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`

# 'Step size'

hv=`echo $u2 $u1 | awk '{printf "%13.9f\n", $1-$2}'`

# echo $u1 $u2 $u3

sed "s/lpar1/$u1/" skin.proto | sed "s/lpar2/$u2/" | sed "s/lpar3/$u3/" > skin.2
sed "s/upar1/$u1p/" skin.2 | sed "s/upar2/$u2p/" | sed "s/upar3/$u3p/" > skin.3
sed "s/mpar1/$u1m/" skin.3 | sed "s/mpar2/$u2m/" | sed "s/mpar3/$u3m/" > SKIN

$static > /dev/null

# Print the energies to the screen:

cat SKENG

# and to the file:

cat SKENG >> $skfile

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

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

# 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:

unew=`echo $u2 $ep $e2p | awk '{printf "%13.9f\n", $1-$2/$3}'`

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

echo $xlabel $unew $enew

# Start again:
# Central value of u.

u2=`echo $unew | awk '{printf "%13.9f\n", $1}'`

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

u1=`echo $u2 | awk '{printf "%13.9f\n", $1-.001}'`
u3=`echo $u2 | awk '{printf "%13.9f\n", $1+.001}'`

# For the internal positions we need ui+1/3, and 2/3-ui

u1p=`echo $u1 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u2p=`echo $u2 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u3p=`echo $u3 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u1m=`echo $u1 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`
u2m=`echo $u2 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`
u3m=`echo $u3 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`

# 'Step size'

hv=`echo $u2 $u1 | awk '{printf "%13.9f\n", $1-$2}'`

# echo $u1 $u2 $u3

sed "s/lpar1/$u1/" skin.proto | sed "s/lpar2/$u2/" | sed "s/lpar3/$u3/" > skin.2
sed "s/upar1/$u1p/" skin.2 | sed "s/upar2/$u2p/" | sed "s/upar3/$u3p/" > skin.3
sed "s/mpar1/$u1m/" skin.3 | sed "s/mpar2/$u2m/" | sed "s/mpar3/$u3m/" > SKIN

$static > /dev/null

# Print the energies to the screen:

cat SKENG

# and to the file:

cat SKENG >> $skfile

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

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

# 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:

unew=`echo $u2 $ep $e2p | awk '{printf "%13.9f\n", $1-$2/$3}'`

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

echo $xlabel $unew $enew

# Start again:
# Central value of u.

u2=`echo $unew | awk '{printf "%13.9f\n", $1}'`

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

u1=`echo $u2 | awk '{printf "%13.9f\n", $1-.001}'`
u3=`echo $u2 | awk '{printf "%13.9f\n", $1+.001}'`

# For the internal positions we need ui+1/3, and 2/3-ui

u1p=`echo $u1 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u2p=`echo $u2 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u3p=`echo $u3 | awk '{printf "%13.9f\n", 0.33333333333333333+$1}'`
u1m=`echo $u1 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`
u2m=`echo $u2 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`
u3m=`echo $u3 | awk '{printf "%13.9f\n", 0.66666666666666667-$1}'`

# 'Step size'

hv=`echo $u2 $u1 | awk '{printf "%13.9f\n", $1-$2}'`

# echo $u1 $u2 $u3

sed "s/lpar1/$u1/" skin.proto | sed "s/lpar2/$u2/" | sed "s/lpar3/$u3/" > skin.2
sed "s/upar1/$u1p/" skin.2 | sed "s/upar2/$u2p/" | sed "s/upar3/$u3p/" > skin.3
sed "s/mpar1/$u1m/" skin.3 | sed "s/mpar2/$u2m/" | sed "s/mpar3/$u3m/" > SKIN

$static > /dev/null

# Print the energies to the screen:

cat SKENG

# and to the file:

cat SKENG >> $skfile

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

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

# 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:

unew=`echo $u2 $ep $e2p | awk '{printf "%13.9f\n", $1-$2/$3}'`

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

echo $xlabel $unew $enew

# Clean up after ourselves:

rm SKIN skin.proto skin.2 skin.3 SKENG SKOUT