#!/usr/bin/python
# -*- coding: utf-8 -*-
#
# fit data coming in on stdin to Kramers' k(F) over a cusp.
#   k ≅ ω²/2πγ √(πβE_b) e^(-βE_b)
#   E_b(F) = ½ω²Δx² - FΔx - (½+1/ω²)F²
# fitting gamma, omega, T, and Δx
# output:
#     <convergence information>
#   gamma:\t<some number>
#   omega:\t<another number>
#   temp:\t<yet another number>
#   delta x:\t<a final number>
#
# usage: fit_kramers
# for example:
#   $ cat datafile | fit_kramers

import sys
sys.path.append('../common/')
import fit
from scipy import pi, sqrt, exp

kB=1.3806503e-23 # m²·kg/s²·K

class Model (fit.Model) :
    def printModel(self) :
        print "k ≅ ω²/2πγ √(πβE_b) e^(-βE_b) with E_b(F) = ½ω²Δx² - FΔx - (½+1/ω²)F²"
    def printParams(self, xopt) :
        print "gamma:\t%g" % xopt[0]
        print "omega:\t%g" % xopt[1]
        print "temp:\t%g" % xopt[2]
        print "delta x:\t%g" % xopt[3]
    def model(self, x, params) :
        F = x
        gamma = params[0]
        omega = params[1]
        T = params[2]
        dX = params[3]
        betaEb = (0.5*(omega**2)*(dX**2) - F*dX - (0.5 + 1.0/omega**2)*F**2)/(kB*T)
        return (omega**2/(2*pi*gamma) * sqrt(pi*betaEb) * exp(-betaEb))
    def guessInitialParams(self) :
        return [1, 1, 1.0/kB, 1.0]

if __name__ == "__main__" :
    fit.run(4, Model)