from scipy import exp, infty
from scipy.integrate import quad
from scipy.special import gamma
from scipy.optimize import brentq

# GHG, 06-08-2014:
#
# This file is very incomplete.  It shows some code snippets which can be
# helpful in dealing with the BEC problem.
# 
# Assuming that you have scipy installed, running this file as
#
#          python BEC-hint.py
#
# at the command line should give you results that look like the following.
#
# ... some possible warning message about exp ...
# Evaluating:  g (z=1, m=1.5) = 2.61237534869
# Evaluating:  g (z=0.9, m=1.5) = 1.61443852857
# Evaluating:  g (z=0., m=1.5) = 0.0
# Evaluating:  fugacity (3, N = 1000) = 0.296001870486
# Evaluating:  fugacity (1. / g(z=1), N = 1000) = 0.991688094173
# Evaluating:  fugacity (1. / g(z=1), N = 1e5) = 0.999619385047

def fugacity (nl3inv, N = 1.e5):
	assert N >= 100
	assert nl3inv >= 0
	gs_val = 1 / (1 + 1. / N)
	if nl3inv == 0:
		return gs_val
	# The tolerance for brentq might need tweaking if high precision is
	# required.
	return brentq (lambda x, nl3inv = nl3inv, N = N:
			1. - g (x, m = 1.5) * nl3inv - x / N / (1. - x),
			0, gs_val)

# This function calculates the g_m function, as required for the analysis
# of the Bose gas problem.
def g (z, m = 1.5):
	if z == 0:
		return 0.0
	# The tolerance (epsabs, epsrel) for quad might need tweaking when high
	# precision is required.
	return quad (lambda x, z = z, m = m: x ** (m - 1.) / (exp (x) / z - 1.)
			/ gamma (m), 0,  infty) [0]

if __name__ == '__main__':

	def eval_str (s):
		print 'Evaluating: ', s,
		try:
			ans = eval (s)
			print '=', ans
		except:
			print
			print '=' * 80
			print '** An error occured---reporting.'
			print '=' * 80
			import traceback
			traceback.print_exc ()
			print '=' * 80
	
	for s in ['g (z=1, m=1.5)',
			'g (z=0.9, m=1.5)',
			'g (z=0., m=1.5)',
			'fugacity (3, N = 1000)',
			'fugacity (1. / g(z=1), N = 1000)',
			'fugacity (1. / g(z=1), N = 1e5)']:
		eval_str (s)