Merge branch 'prerelease-3.2.1'

[unres.git] / ctest / matplotlib_fit_hist.py

diff --git a/ctest/matplotlib_fit_hist.py b/ctest/matplotlib_fit_hist.py
new file mode 100755 (executable)
index 0000000..754540c
--- /dev/null
+++ b/ctest/matplotlib_fit_hist.py
@@ -0,0 +1,50 @@
+#! /usr/bin/env python
+
+import matplotlib
+#matplotlib.use('GTK')
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.optimize import curve_fit
+from math import sqrt
+import sys
+
+def prob_T(x,a):
+    gg=np.float128(g)
+    aa=np.float128(10**(-gg-2)*a)
+    Tr=np.float128(300.)
+    return np.exp( np.log(aa) + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
+#    return aa * (  x**((gg-2)/2) * np.exp( -gg*x/(2*Tr) ) )
+
+#x,y= np.loadtxt('1L2Y_L_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
+#x,y= np.loadtxt('1L2Y_NH_GB000.stat',usecols=(0,11),skiprows=10000,unpack=True)
+#x,y= np.loadtxt('1L2Y_B_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
+g=int(sys.argv[2])
+if (sys.argv[1] == '1L2Y_NH_GB000.stat' or sys.argv[1] == '1L2Y_NH_GB.stat'):
+ x,y= np.loadtxt(sys.argv[1],usecols=(0,11),skiprows=10000,unpack=True)
+else: 
+ x,y= np.loadtxt(sys.argv[1],usecols=(0,10),skiprows=10000,unpack=True)
+
+h,bin=np.histogram(y,bins=50,density=True)
+
+plt.bar(bin[:-1], h, width = bin[2]-bin[1])
+plt.xlim(min(bin), max(bin))
+plt.ylim(0,max(h))
+plt.ylabel('probality')
+plt.xlabel('temperature')
+
+center = (bin[:-1] + bin[1:]) / 2
+#print bin
+#print center
+popt, pcov = curve_fit(prob_T, center, h)
+xfine = np.linspace(min(bin), max(bin), 100)
+#print popt
+
+chi_squared = np.sum((prob_T(center, *popt)-h)**2)
+print '%15.10f' % chi_squared
+
+#print xfine
+#print prob_T(xfine,popt[0])
+plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
+plt.savefig(sys.argv[1]+'.png')
+#plt.show()