# Path: lmatools/lasso/empirical_charge_density.py
# class rho_retrieve(object):
#     def __init__(self, area, d, zinit, separation, constant=False, arbitrary_rho=0.4e-9):
#     def rho_br(self, z, separation):
#     def energy_estimation(self, rho, d, area):
#     def calculate(self):


import numpy as np
import logging
import re
from scipy.spatial import Delaunay, ConvexHull
from scipy.special import factorial
from scipy.spatial.qhull import QhullError
from lmatools.lasso import empirical_charge_density as cd
        
    tri = Delaunay(xyz[:,0:3])
    vertices = tri.points[tri.vertices]
    
    # This is the volume formula in 
    # https://github.com/scipy/scipy/blob/master/scipy/spatial/tests/test_qhull.py#L106
    # Except the formula needs to be divided by ndim! to get the volume, cf., 
    # http://en.wikipedia.org/wiki/Simplex#Geometric_properties
    # Credit Pauli Virtanen, Oct 14, 2012, scipy-user list
    
    q = vertices[:,:-1,:] - vertices[:,-1,None,:]
    simplex_volumes = (1.0 / factorial(q.shape[-1])) * np.fromiter(
           (np.linalg.det(q[k,:,:]) for k in range(tri.nsimplex)) , dtype=float)
    # print vertices.shape # number of simplices, points per simplex, coords
    # print q.shape
    
    # The simplex volumes have negative values since they are oriented 
    # (think surface normal direction for a triangle
    volume=np.sum(np.abs(simplex_volumes))
    return volume, vertices, simplex_volumes

##############ADDED 01/05/2017 ###############
def energy(area, separation, zinit, constant, eta):
    #Charge separation computed from 27th and 73rd percentiles of 
    #flash altitude source locations - marks where the most sources are typically
    #found from synthetic flashes generated in the NSSL COMMAS.
    #
    #eta = 0.01 is recommended and is a ballpark neutrlization efficiency as found in Salinas et al. [In Progress - 060220]

    distance = separation #np.abs(random)

# Path: lmatools/lasso/empirical_charge_density.py
# class rho_retrieve(object):
#     def __init__(self, area, d, zinit, separation, constant=False, arbitrary_rho=0.4e-9):
#     def rho_br(self, z, separation):
#     def energy_estimation(self, rho, d, area):
#     def calculate(self):


# Path: lmatools/flashsort/flash_stats.py
import numpy as np
import logging
import re
from scipy.spatial import Delaunay, ConvexHull
from scipy.special import factorial
from scipy.spatial.qhull import QhullError
from lmatools.lasso import empirical_charge_density as cd

        
    tri = Delaunay(xyz[:,0:3])
    vertices = tri.points[tri.vertices]
    
    # This is the volume formula in 
    # https://github.com/scipy/scipy/blob/master/scipy/spatial/tests/test_qhull.py#L106
    # Except the formula needs to be divided by ndim! to get the volume, cf., 
    # http://en.wikipedia.org/wiki/Simplex#Geometric_properties
    # Credit Pauli Virtanen, Oct 14, 2012, scipy-user list
    
    q = vertices[:,:-1,:] - vertices[:,-1,None,:]
    simplex_volumes = (1.0 / factorial(q.shape[-1])) * np.fromiter(
           (np.linalg.det(q[k,:,:]) for k in range(tri.nsimplex)) , dtype=float)
    # print vertices.shape # number of simplices, points per simplex, coords
    # print q.shape
    
    # The simplex volumes have negative values since they are oriented 
    # (think surface normal direction for a triangle
    volume=np.sum(np.abs(simplex_volumes))
    return volume, vertices, simplex_volumes

##############ADDED 01/05/2017 ###############
def energy(area, separation, zinit, constant, eta):
    #Charge separation computed from 27th and 73rd percentiles of 
    #flash altitude source locations - marks where the most sources are typically
    #found from synthetic flashes generated in the NSSL COMMAS.
    #
    #eta = 0.01 is recommended and is a ballpark neutrlization efficiency as found in Salinas et al. [In Progress - 060220]

    distance = separation #np.abs(random)

    density  = cd.rho_retrieve(area, distance, zinit, separation, False, None) #None - No constant charge density specified