Skip to content

Commit

Permalink
docs: more in-code documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
tfrederiksen committed Mar 11, 2018
1 parent eb2f8d6 commit 4f19628
Showing 1 changed file with 52 additions and 20 deletions.
72 changes: 52 additions & 20 deletions Inelastica/math/misc.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,10 @@
import numpy.linalg as LA


##############################################################

def outerAdd(* args):
"""
Returns :math:`A_{ijk}=B_i+C_j+D_k`.
"""
# A_ijk=B_i+C_j+D_k
tmp = args[0].copy()
for ii in range(1, len(args)):
Expand All @@ -13,39 +14,57 @@ def outerAdd(* args):


def dist(x):
"""
Euclidian norm :math:`||x||` of a vector :math:`x`.
"""
return N.sqrt(N.dot(x, x))


def mysqrt(x):
# Square root of matrix
"Square root of matrix."
ev, U = LA.eig(x)
U = N.transpose(U)

tmp = N.zeros((len(ev), len(ev)), N.complex)
for ii in range(len(ev)):
tmp[ii, ii] = N.sqrt(ev[ii])
for ii,evi in enumerate(ev):
tmp[ii, ii] = N.sqrt(evi)

return mm(LA.inv(U), tmp, U)

##############################################################


def fermi(mu, E, kT):
"""
Fermi function :math:`n_F(\mu) = 1/[exp((E-\mu)/k_BT)+1]`.
"""
return 1/(N.exp(N.clip((E-mu)/kT, -70.0, 70.0))+1)


def box(mu1, mu2, grid, kT):
"""
Window (box) function defined as
:math:`n_F(\mu_2)-n_F(\mu_1)`
"""
# f2-f1 (Box!)
return fermi(mu2, grid, kT)-fermi(mu1, grid, kT)

##############################################################


def trapez(x, f, equidistant=False):
'''
"""
Integration of vector f on grid x using the 3rd degree polynomial.
The grid x does not have to be equidistant.
'''
Parameters
----------
x : ndarray
f : ndarray
equidistant : bool
`False` = 3rd degree polynomial method.
`True` = Linear trapez method.
Returns
-------
res : complex number
Result of the integration
"""
if equidistant:
# Trapez method!
d = N.array((x[1]-x[0])*N.ones(len(x)), N.complex)
Expand All @@ -54,11 +73,11 @@ def trapez(x, f, equidistant=False):
return N.dot(d, f)
else:
# 3rd degree polynomial except for 1st and last bins
sum = (x[1]-x[0])*(f[0]+f[1])/2+(x[-1]-x[-2])*(f[-1]+f[-2])/2
res = (x[1]-x[0])*(f[0]+f[1])/2+(x[-1]-x[-2])*(f[-1]+f[-2])/2
for ii in range(1, len(x)-2):
x0, x1, x2, x3=x[ii-1], x[ii], x[ii+1], x[ii+2]
y0, y1, y2, y3=f[ii-1], f[ii], f[ii+1], f[ii+2]
sum+=((x1-x2)*(-6*x0**2*(x0-x3)*x3**2*(y1+y2)+4*x0*x2*(x0-x3)*x3*(x0+x3)*(2*y1+y2)+\
x0, x1, x2, x3 = x[ii-1], x[ii], x[ii+1], x[ii+2]
y0, y1, y2, y3 = f[ii-1], f[ii], f[ii+1], f[ii+2]
res += ((x1-x2)*(-6*x0**2*(x0-x3)*x3**2*(y1+y2)+4*x0*x2*(x0-x3)*x3*(x0+x3)*(2*y1+y2)+\
3*x2**3*(x3**2*(y0-y1)+x0**2*(y1-y3))+\
x1**3*(-2*x2*(x3*(y0-y2)+x0*(y2-y3))+3*(x3**2*(y0-y2)+x0**2*(y2-y3))+\
x2**2*(-y0+y3))+x2**4*(x3*(-y0+y1)+x0*(-y1+y3))+\
Expand All @@ -73,14 +92,28 @@ def trapez(x, f, equidistant=False):
3*x2*(x3**2*(y0+y1+2*y2)-x0**2*(y1+2*y2+y3))+\
3*x2**2*(x3*(2*y0+y1+y2)-x0*(y1+y2+2*y3)))))/\
(12.*(x0-x1)*(x0-x2)*(x0-x3)*(x1-x3)*(x2-x3))
return sum
##############################################################
return res


def interpolate(nx, x, y):
"""
Interpolate f(x)=y to find f(nx)
Makes no checks for nx inside x region!!!
Interpolate :math:`f(x)=y` to find :math:`f(nx)`. Extrapolation allowed.
NB: Makes no checks for nx inside x region!!!
Parameters
----------
nx : ndarray
Abscissa values to interpolate the function.
x : ndarray
Known abscissa values.
y : ndarray
Known function values.
Returns
-------
ny : ndarray
Interpolated function values on nx.
"""
ny = N.array([0.0]*len(nx))
Lpos = N.searchsorted(x, nx)
Expand All @@ -91,4 +124,3 @@ def interpolate(nx, x, y):
#TF: NB EXTRAPOLATION condition added!
ny[ix] = y[-1]+(y[-1]-y[-2])/(x[-1]-x[-2])*(nx[ix]-x[-1])
return ny

0 comments on commit 4f19628

Please sign in to comment.