text: 
  """
  This file provides a way to obtain thermodynamic quantities from an 
  interpolation of available NLCE solutions 
  """
  import ldaconf
  import numpy as np
  import glob
  
  import matplotlib.pyplot as plt
  import matplotlib
  from matplotlib import rc
  rc('font', **{'family':'serif'})
  rc('text', usetex=True)
  
  from scipy.spatial import Delaunay
  from scipy.interpolate import CloughTocher2DInterpolator, LinearNDInterpolator
  from scipy.interpolate.interpnd import _ndim_coords_from_arrays
  
  import logging
  # create logger
  logger = logging.getLogger(__name__)
  logger.addHandler(logging.NullHandler())
  #logger.disabled = True
  
  from qmc import get_qty_mu
  
  def find_closest_nlce( U=8, T=0.67, mu=4., qty='dens', **kwargs):
      """
      This function finds the closest values of U and T in the NLCE data 
      that straddle the values U and T given as arguments.
      """
  
      msg0 = "U={:0.2f}, T={:0.2f}".format( U, T) 
     
      nUs = 3 
      us = [ float(u.split('/U')[-1]) for u in \
               glob.glob( ldaconf.basedir + 'NLCE8_FinalK/U*' ) ] 
      du = [ np.abs( U-u ) for u in us ] 
      index = np.argsort( du ) 
      Ulist0 = range( nUs ) 
      Upts = [ us[index[i]] for i in Ulist0 ]  
   
      
      # The T points are not uniformly spaced so we find the two closest ones
      # We start with a list of available T points: 
      Ts = np.array( sorted( [ float(g.split('/T')[1].split('.dat')[0]) for g in \
               glob.glob(ldaconf.basedir + 'NLCE8_FinalK/U00/T*') ] ) )
      
      
      diff = T-Ts
  
      error = False
  
      if np.all( diff < 0 ):
          print "Available temperatures do not make it this low:"
          print " T = ", T
          error = True  
  
      if not error:    
          order_pos =  np.argsort(np.abs( diff[diff>0] ))
          order_neg =  np.argsort(np.abs( diff[diff<0] ))
          Tpts = sorted( [  Ts[diff>0][ order_pos[0] ] ,  Ts[diff<0][ order_neg[0]] ] )  
      else:
          order = np.argsort( np.abs( diff ) ) 
          Tpts = sorted( [ Ts[order[0]], Ts[order[1]], Ts[order[2]] ] ) 
          #Ta = min(Ts[order[0]], Ts[order[1]])
          #Tb = max(Ts[order[0]], Ts[order[1]])
          #print "T in ", Ta, Tb
      
      datadir = ldaconf.basedir
      datfiles = []
      for Uval in Upts:
          for Tval in Tpts:
              fname =  datadir + \
                  'NLCE8_FinalK/U{:02d}/T{:0.2f}.dat'.format(int(Uval),Tval)
              datfiles.append([ fname, Uval, Tval ])
              
      if qty == 'dens':
          COL = 1 
      elif qty == 'entr':
          COL = 2 
      elif qty == 'spi':
          COL = 3 
      elif qty == 'kappa':
          COL = 4 
      else:
          raise "Qty not defined:", qty
              
      MUCOL = 0 
      basedat = [] 
  
      qtyinterp = kwargs.get( 'qtyinterp', 'nearest' )
  
      for f in datfiles:
          msg = 'U={:0.2f}, T={:0.2f}'.format(U,T) + \
             ' mu={:0.2f}, Upt={:0.3f}, Tpt={:0.3f}'.\
             format(mu, f[1], f[2])
  
          try:
              dat = np.loadtxt(f[0])
              # careful = False, for NLCE data it is ok to not worry about
              # chemical potentials that are outside the data range
              qtyresult = get_qty_mu( dat, mu, MUCOL, COL, msg=msg, careful=False) 
              if qtyresult == 'out-of-bounds':
                  print msg
                  print "out-of-bounds"
                  continue
   
              basedat.append( [f[1], f[2], qtyresult] )
                  
          except Exception as e:
              print "Failed to get data from file = ", f
  
              print msg
  
              fig = plt.figure( figsize=(3.5,3.5))
              gs = matplotlib.gridspec.GridSpec( 1,1 ,\
                      left=0.15, right=0.96, bottom=0.12, top=0.88)
              ax = fig.add_subplot( gs[0] )
              ax.grid(alpha=0.5)
              ax.plot( dat[:,MUCOL], dat[:,COL], '.-')
              ax.axvline( mu )
  
              ax.text( 0.5, 1.05, msg, ha='center', va='bottom', \
                  transform=ax.transAxes)
              if matplotlib.get_backend() == 'agg':
                  fig.savefig('err_mu.png', dpi=200) 
              else:
                  plt.show()
  
              raise e 
          
      basedat =   np.array(basedat)
      #print "Closest dat = ", basedat
      points = _ndim_coords_from_arrays(( basedat[:,0] , basedat[:,1]))
      
      
      
      #finterp = CloughTocher2DInterpolator(points, basedat[:,2])
      finterp = LinearNDInterpolator( points, basedat[:,2] )
      
      
      try:
          result = finterp( U,T )
          if np.isnan(result):
              if U >= 30 and U<=32.5:
                  result = finterp( 29.99, T ) 
                  logger.warning(" nlce: U={:0.1f} replaced to U=29.99 ".format(U) )
          if np.isnan(result):
              raise Exception("!!!! nlce: Invalid result !!!!\n" + msg0)
          
      except Exception as e:
          print e
          error = True
          
      if error or kwargs.get('showinterp',False):
          #print "Interp points:"
          #print basedat
          
          tri = Delaunay(points)
          fig = plt.figure( figsize=(3.5,3.5))
          gs = matplotlib.gridspec.GridSpec( 1,1 ,\
                  left=0.15, right=0.96, bottom=0.12, top=0.88)
          ax = fig.add_subplot( gs[0] )
          ax.grid(alpha=0.5)
          ax.triplot(points[:,0], points[:,1], tri.simplices.copy())
          ax.plot(points[:,0], points[:,1], 'o')
          ax.plot( U, T, 'o', ms=6., color='red')
          xlim = ax.get_xlim()
          dx = (xlim[1]-xlim[0])/10.
          ax.set_xlim( xlim[0]-dx, xlim[1]+dx )
          ylim = ax.get_ylim()
          dy = (ylim[1]-ylim[0])/10.
          ax.set_ylim( ylim[0]-dy, ylim[1]+dy )
          ax.set_xlabel('$U/t$')
          ax.set_ylabel('$T/t$',rotation=0,labelpad=8)
          
          tt = kwargs.get('title_text','')
          ax.set_title( tt + '$U/t={:.2f}$'.format(U) + ',\ \ ' + '$T/t={:.2f}$'.format(T), \
                  ha='center', va='bottom', fontsize=10)
          save_err = kwargs.get('save_err',None) 
          if save_err is not None:
              print "saving png to ", save_err 
              fig.savefig( save_err, dpi=300)
  
          if matplotlib.get_backend() == 'agg':
              print "saving png to err.png"
              fig.savefig('err.png', dpi=200) 
          else:
              plt.show()
  
          raise
      
      return result
  
  
  QTYINTERP = 'linear' 
  
  def nlce_dens( T, t, mu, U, ignoreLowT=False, verbose=True):
      U_ = U/t 
      T_ = T/t 
      mu_ = mu/t   
  
      result = np.empty_like(mu) 
      for i in range( len(mu_)):
          result[i] = find_closest_nlce( U=U_[i], T=T_[i], mu=mu_[i], \
                       qty='dens',  qtyinterp=QTYINTERP  ) 
      return result    
   
  def nlce_entr( T, t, mu, U, ignoreLowT=False, verbose=True):
      U_ = U/t 
      T_ = T/t 
      mu_ = mu/t   
  
      result = np.empty_like(mu) 
      for i in range( len(mu_)):
          result[i] = find_closest_nlce( U=U_[i], T=T_[i], mu=mu_[i], \
                       qty='entr', qtyinterp=QTYINTERP ) 
      return result     
  
  def nlce_spi( T, t, mu, U, ignoreLowT=False, verbose=True):
      U_ = U/t 
      T_ = T/t 
      mu_ = mu/t   
  
      result = np.empty_like(mu) 
      for i in range( len(mu_)):
          result[i] = find_closest_nlce( U=U_[i], T=T_[i], mu=mu_[i], \
                       qty='spi', qtyinterp=QTYINTERP ) 
      return result     
   
  def nlce_cmpr( T, t, mu, U, ignoreLowT=False, verbose=True):
      U_ = U/t 
      T_ = T/t 
      mu_ = mu/t   
  
      result = np.empty_like(mu) 
      for i in range( len(mu_)):
          result[i] = find_closest_nlce( U=U_[i], T=T_[i], mu=mu_[i], \
                       qty='kappa',  qtyinterp=QTYINTERP  ) 
      return result     
  
  