# Author: G.-H. Sam Gweon
# Date: 06-04-2014
#
# This is the code, using which homework solutions for problem 7.1 were
# computed.  This module uses many of my private custom modules (data, gmath,
# reCreateable, util) and so this python module will not be useable right
# away on a random computer.  However, this code should give a good sense of
# how things were calculated.

from data import DataXY
from gmath import brent
from Numeric import exp
from reCreateable import asDefinedReCreateable
from scipy import infty
from scipy.integrate import quad
from scipy.special import gamma
from util import toArray

def _fugacity (nl3inv, N = 1.e5): # <<<
	assert N >= 100
	assert nl3inv >= 0
	gs_val = 1 / (1 + 1. / N)
	if nl3inv == 0:
		return gs_val
	# brent is a zero finding algorithm.
	# TODO: The tolerance (xtol) for brent might need tweaking if higher
	#       precision is required.  This code is OK up to N ~ 1e9.
	return brent (lambda x, nl3inv = nl3inv, N = N:
			1. - g (x, m = 1.5, dxy = False) * nl3inv - x / N / (1. - x),
			0, gs_val)
# >>>
def _g (z, m = 1.5): # <<<
	if z == 0:
		return 0.0
	# TODO: The tolerance (epsabs, epsrel) for quad might need tweaking when
	#       higher precision is required.  This code is OK up to N ~ 1e9.
	return quad (lambda x, z = z, m = m: x ** (m - 1.) / (exp (x) / z - 1.)
			/ gamma (m), 0,  infty) [0]
# >>>

def BECfrac (nl3inv, N = 1.e5, dxy = True): # <<<
	"""
	This returns the condensate fraction n(0)/N.
	"""
	dargs = locals ().copy () if dxy else None
	z = toArray (fugacity (nl3inv, N = N, dxy = False), 'd')
	y = z / (1. - z) / N
	if type (nl3inv) in (int, float):
		return y
	else:
		return asDefinedReCreateable (DataXY (nl3inv, y,
			name = 'BECfrac (N=%g)' % N), None, BECfrac, dargs) if dxy else y
# >>>
def fugacity (nl3inv, N = 1.e5, dxy = True): # <<<
	dargs = locals ().copy () if dxy else None
	func = _fugacity
	if type (nl3inv) in (int, float):
		return func (nl3inv, N = N)
	else:
		y = [func (nl3invval, N = N) for nl3invval in nl3inv]
		return asDefinedReCreateable (DataXY (nl3inv, y, name =
			'fugacity (N=%g)' % N), None, fugacity, dargs) if dxy else y
# >>>
def g (z, m = 1.5, dxy = True): # <<<
	"""
	The "g function" for the Bose gas, which is well-defined for any
	0 <= z < 1 for any value of m, but diverges for z = 1 and m <= 1.
	"""
	dargs = locals ().copy () if dxy else None
	func = _g
	if type (z) in (int, float):
		return func (z, m = m)
	else:
		y = [func (zval, m = m) for zval in z]
		return asDefinedReCreateable (DataXY (z, y,
			name = 'g_{%g} function' % m), None, g, dargs) if dxy else y
# >>>