test/woods-saxon/plot-woods-saxon.py

0001 #!/usr/bin/env python3
0002
0003 """
0004 Compile the C++ Woods-Saxon generator, run it, pipe the results into this
0005 script, and plot.
0006 """
0007
0008 import os.path
0009 import subprocess
0010
0011 import matplotlib.pyplot as plt
0012 import numpy as np
0013 from scipy import integrate
0014
0015
0016 def woods_saxon(r, R, a):
0017     return r*r/(1. + np.exp((r-R)/a))
0018
0019
0020 def main():
0021     basename = 'generate-woods-saxon'
0022     cxxname = '{}.cxx'.format(basename)
0023     if not os.path.exists(basename) or (
0024             os.path.getmtime(cxxname) > os.path.getmtime(basename)):
0025         print('compiling C++ Woods-Saxon generator')
0026         subprocess.check_call(
0027             ['g++', '-std=c++11', '-Wall', '-Wextra', '-O3', '-march=native',
0028              cxxname, '-o'+basename])
0029
0030     R = 6.5
0031     a = 0.5
0032     N = int(1e7)
0033
0034     print('generating Woods-Saxon numbers')
0035     with subprocess.Popen(['./'+basename, str(R), str(a), str(N)],
0036                           stdout=subprocess.PIPE) as proc:
0037         samples = np.fromiter(proc.stdout, dtype=float, count=N)
0038
0039     print('plotting')
0040     plt.hist(samples, bins=200, histtype='step', normed=True)
0041
0042     rmax = samples.max()
0043     r = np.linspace(0, rmax, 1000)
0044     norm = integrate.quad(woods_saxon, 0, rmax, args=(R, a))[0]
0045     plt.plot(r, woods_saxon(r, R, a)/norm)
0046
0047     plt.xlim(0, rmax)
0048     plt.show()
0049
0050
0051 if __name__ == "__main__":
0052     main()