+#!/bin/bash
+# -*- coding: utf-8 -*-
+#
+# Test the Kramers integral k(F) implementation on the cusp potential.
+#
+# The test is escape from a cusp potential from Kramers' 1940 paper.
+# He solves the cusp potential in the text immediately after Eq. 17 & Fig. 2
+# We can translate Kramers' notation to Hänggi's more modern notation with
+# Kramers Hängii Meaning Why I think so
+# 2πω ω Curvature at minimum Hänggi's footnote to Eq. 1.5
+# η γ Viscosity/damping Compare H. Eq.1.5 with K. Eq. 25
+# Q E_b Energy barrier height K. Fig. 1 (p 291)
+# T k_B·T≡1/β Absolute temperature K: “k_B = 1” (text after Eq. 5)
+# r k Reaction rate seems more reasonable than k ;)
+# With these translations, Kramers' Eq. is
+# k ≅ ω²/2πγ √(πβE_b) e^(-βE_b)
+# for E(x) = ½ω²(|x|-Δx)².
+# The barrier is at x_b=0 with minimums at x_ac=±Δx.
+# The barrier height E_b = E(0) = ½ω²Δx²
+#
+# As a force is applied to tip the potential,
+# E(F,x) = ½ω²(|x|-Δx)²-Fx.
+# Focusing only on the area left of the barrier (x<0), we have
+# E_L(F,x) = ½ω²(-x-Δx)²-Fx = ½ω²(x+Δx)²-Fx
+# E_L'(F,x) = ω²(x+Δx)-F
+# E_L"(F,x) = ω²,
+# just like the F=0 case. Kramers' solution still holds then, but with a
+# different barrier height. The minimum moves to
+# x_a = F/ω² - Δx.
+# so
+# Eb(F) = E(F,0) - E(F,x_a) = E(0,0) - E(F,x_a)
+# = ½ω²Δx² - ½ω²[(F/ω² - Δx)+Δx]²-F(F/ω² - Δx)
+# = ½ω²Δx² - ½F² - F²/ω² - FΔx
+# = ½ω²Δx² - FΔx - (½+1/ω²)F²
+# ----+----- ^- reduction due to decreasing barrier distance
+# ^- Bell model
+#
+# Usage: kramers_rate.sh
+
+# Since we're on the cluster, we'll go crazy and get really good statistics ;)
+FMAX=1
+K_UTILS=../../bin/k_model_utils
+DATA=kramers_rate.d
+HIST=kramers_rate.hist
+GPFILE=kramers_rate.gp
+PNGFILE=kramers_rate.png
+
+# set up the physics
+# E(x) = ½ω²(|x|-Δx)²
+OMEGA=3
+GAMMA=100
+KB=1.3806503e-23 # m²·kg/s²·K
+T=`python -c "print (1 / ($KB))"` # pick units where k_B = 1
+DX=1
+EB=`python -c "print (0.5 * (($OMEGA)*($DX))**2)"`
+
+# set up the theory
+E_FN="0.5*($OMEGA)**2 * (abs(x)-($DX))**2"
+BETA=`python -c "print (1.0 / (($KB)*($T)))"`
+# k ≅ ω²/2πγ √(πβE_b) e^(-βE_b)
+BETAEB=`python -c "print (($BETA)*($EB))"`
+THEORY_K="BETAEB(F) = ($BETA)*(($EB) - (0.5+1/($OMEGA)**2)*F**2 - F*($DX))
+K(F) = ($OMEGA)**2/(2*pi*($GAMMA)) * sqrt(pi*BETAEB(F)) * exp(-BETAEB(F))"
+
+# Get the diffusion coefficient from the Einstein-Smoluchowski relation
+# D = µ_p·k_B·T = k_B·T/γ
+# where µ_p=1/γ is the particle mobility
+# See, for example,
+# http://en.wikipedia.org/wiki/Einstein_relation_(kinetic_theory)
+D=`python -c "print (1/(($BETA)*($GAMMA)))"`
+# pick safe integration bounds
+XMIN=`python -c "print (-2*($DX))"`
+XMAX="$DX"
+# set up the cubic-spline approximation to a sharp cusp… get close anyway…
+SPLINE_PTS=`python -c "from sys import stdout
+dx = (($XMAX)-($XMIN))/100.0
+stdout.write('{')
+for i in range(100) :
+ x = ($XMIN)+i*dx
+ stdout.write('(%g,%g),' % (x, ($E_FN)))
+stdout.write('(%g,%g)}\n' % (($XMAX), ($E_FN)))"`
+
+# Run the simulation
+echo "$K_UTILS -kkramers_integ -K \"$D,$SPLINE_PTS,$XMIN,$XMAX\" -x$XMIN -X$XMAX -T$T -F$FMAX > $DATA"
+$K_UTILS -kkramers_integ -K "$D,$SPLINE_PTS,$XMIN,$XMAX" -x$XMIN -X$XMAX -T$T -F$FMAX > $DATA
+
+# fit the data
+FIT=`cat $DATA | ./fit_kramers`
+FITGAMMA=`echo "$FIT" | sed -n 's/gamma:[[:space:]]*//p'`
+FITOMEGA=`echo "$FIT" | sed -n 's/omega:[[:space:]]*//p'`
+FITT=`echo "$FIT" | sed -n 's/temp:[[:space:]]*//p'`
+FITDX=`echo "$FIT" | sed -n 's/delta x:[[:space:]]*//p'`
+echo -e "gamma:\t$FITGAMMA\t(expected $GAMMA)"
+echo -e "omega:\t$FITOMEGA\t(expected $OMEGA)"
+echo -e "temp:\t$FITT\t(expected $T)"
+echo -e "delta x:\t$FITDX\t(expected $DX)"
+
+FITEB=`python -c "print (0.5 * (($FITOMEGA)*($FITDX))**2)"`
+FITE_FN="0.5*($FITOMEGA)**2 * (abs(x)-($FITDX))**2"
+FITBETA=`python -c "print (1.0 / (($KB)*($FITT)))"`
+FITBETAEB=`python -c "print (($FITBETA)*($FITEB))"`
+FIT_K="FITBETAEB(F) = ($FITBETA)*(($FITEB) - (0.5+1/($FITOMEGA)**2)*F**2 - F*($FITDX))
+fitK(F) = ($FITOMEGA)**2/(2*pi*($FITGAMMA)) * sqrt(pi*FITBETAEB(F)) * exp(-FITBETAEB(F))"
+
+# generate gnuplot script
+echo "set terminal png" > $GPFILE
+echo "set output '$PNGFILE'" >> $GPFILE
+echo "set logscale y" >> $GPFILE
+echo "set xlabel 'Time'" >> $GPFILE
+echo "set ylabel 'Deaths'" >> $GPFILE
+echo "$THEORY_K" >> $GPFILE
+echo "$FIT_K" >> $GPFILE
+echo "plot K(x), fitK(x), '$DATA' with points" >> $GPFILE
+
+gnuplot $GPFILE