testing/kramers_rate/fit_kramers

   1 #!/usr/bin/python
   2 # -*- coding: utf-8 -*-
   3 #
   4 # fit data coming in on stdin to Kramers' k(F) over a cusp.
   5 #   k ≅ ω²/2πγ √(πβE_b) e^(-βE_b)
   6 #   E_b(F) = ½ω²Δx² - FΔx - (½+1/ω²)F²
   7 # fitting gamma, omega, T, and Δx
   8 # output:
   9 #     <convergence information>
  10 #   gamma:\t<some number>
  11 #   omega:\t<another number>
  12 #   temp:\t<yet another number>
  13 #   delta x:\t<a final number>
  14 #
  15 # usage: fit_kramers
  16 # for example:
  17 #   $ cat datafile | fit_kramers
  18
  19 import sys
  20 sys.path.append('../common/')
  21 import fit
  22 from scipy import pi, sqrt, exp
  23
  24 kB=1.3806503e-23 # m²·kg/s²·K
  25
  26 class Model (fit.Model) :
  27     def printModel(self) :
  28         print "k ≅ ω²/2πγ √(πβE_b) e^(-βE_b) with E_b(F) = ½ω²Δx² - FΔx - (½+1/ω²)F²"
  29     def printParams(self, xopt) :
  30         print "gamma:\t%g" % xopt[0]
  31         print "omega:\t%g" % xopt[1]
  32         print "temp:\t%g" % xopt[2]
  33         print "delta x:\t%g" % xopt[3]
  34     def model(self, x, params) :
  35         F = x
  36         gamma = params[0]
  37         omega = params[1]
  38         T = params[2]
  39         dX = params[3]
  40         betaEb = (0.5*(omega**2)*(dX**2) - F*dX - (0.5 + 1.0/omega**2)*F**2)/(kB*T)
  41         return (omega**2/(2*pi*gamma) * sqrt(pi*betaEb) * exp(-betaEb))
  42     def guessInitialParams(self) :
  43         return [1, 1, 1.0/kB, 1.0]
  44
  45 if __name__ == "__main__" :
  46     fit.run(4, Model)