OSDN > Buscar Software > Scientific/Engineering > Physics > libtetrabz > SCM > git > python > Commit

libtetrabz
Fork
python
libtetrabz

R/O
HTTP
SSH
HTTPS

python: Commit

libtetrabz python package

Commit MetaInfo

Revisión	0fcf47ce73a8a84f5cd88a53885ff8e38cba82d9 (tree)
Tiempo	2021-10-22 01:55:50
Autor	Mitsuaki Kawamura <kawamitsuaki@gmai...>
Commiter	Mitsuaki Kawamura

Log Message

Port remaining functions

Cambiar Resumen

modified: LICENSE (diff)
modified: README.md (diff)
modified: pyproject.toml (diff)
modified: src/libtetrabz/__init__.py (diff)
add: src/libtetrabz/libtetrabz_dbldelta.py (diff)
add: src/libtetrabz/libtetrabz_dblstep.py (diff)
add: src/libtetrabz/libtetrabz_fermigr.py (diff)
add: src/libtetrabz/libtetrabz_polcmplx.py (diff)
delete: src/libtetrabz/{libtetrabz_dos_mod.py => libtetrabz_dos.py}
delete: src/libtetrabz/{libtetrabz_occ_mod.py => libtetrabz_occ.py}
delete: src/libtetrabz/{libtetrabz_polstat_mod.py => libtetrabz_polstat.py}

Diferencia incremental

--- a/LICENSE

+++ b/LICENSE

		@@ -1,4 +1,4 @@
1		-Copyright (c) 2018 The Python Packaging Authority
	1	+Copyright (c) Copyright (c) 2014 Mitsuaki Kawamura
2	2
3	3	Permission is hereby granted, free of charge, to any person obtaining a copy
4	4	of this software and associated documentation files (the "Software"), to deal

--- a/README.md

+++ b/README.md

		@@ -1,4 +1,4 @@
1		-# Example Package
	1	+# libtetrabz
2	2
3	3	This is a simple example package. You can use
4	4	[Github-flavored Markdown](https://guides.github.com/features/mastering-markdown/)

--- a/pyproject.toml

+++ b/pyproject.toml

		@@ -2,5 +2,6 @@
2	2	requires = [
3	3	"setuptools>=42",
4	4	"wheel"
	5	+ "numpy"
5	6	]
6	7	build-backend = "setuptools.build_meta"

--- a/src/libtetrabz/__init__.py

+++ b/src/libtetrabz/__init__.py

		@@ -0,0 +1,31 @@
	1	+#
	2	+# Copyright (c) 2014 Mitsuaki Kawamura
	3	+#
	4	+# Permission is hereby granted, free of charge, to any person obtaining a
	5	+# copy of this software and associated documentation files (the
	6	+# "Software"), to deal in the Software without restriction, including
	7	+# without limitation the rights to use, copy, modify, merge, publish,
	8	+# distribute, sublicense, and/or sell copies of the Software, and to
	9	+# permit persons to whom the Software is furnished to do so, subject to
	10	+# the following conditions:
	11	+#
	12	+# The above copyright notice and this permission notice shall be included
	13	+# in all copies or substantial portions of the Software.
	14	+#
	15	+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	16	+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
	17	+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
	18	+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
	19	+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
	20	+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
	21	+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
	22	+#
	23	+from libtetrabz_dos import dos
	24	+from libtetrabz_dos import intdos
	25	+from libtetrabz_occ import occ
	26	+from libtetrabz_occ import fermieng
	27	+from libtetrabz_polstat import polstat
	28	+from libtetrabz_fermigr import fermigr
	29	+from libtetrabz_dbldelta import dbldelta
	30	+from libtetrabz_dblstep import dblstep
	31	+from libtetrabz_polcmplx import polcmplx

--- /dev/null

+++ b/src/libtetrabz/libtetrabz_dbldelta.py

		@@ -0,0 +1,144 @@
	1	+#
	2	+# Copyright (c) 2014 Mitsuaki Kawamura
	3	+#
	4	+# Permission is hereby granted, free of charge, to any person obtaining a
	5	+# copy of this software and associated documentation files (the
	6	+# "Software"), to deal in the Software without restriction, including
	7	+# without limitation the rights to use, copy, modify, merge, publish,
	8	+# distribute, sublicense, and/or sell copies of the Software, and to
	9	+# permit persons to whom the Software is furnished to do so, subject to
	10	+# the following conditions:
	11	+#
	12	+# The above copyright notice and this permission notice shall be included
	13	+# in all copies or substantial portions of the Software.
	14	+#
	15	+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	16	+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
	17	+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
	18	+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
	19	+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
	20	+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
	21	+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
	22	+#
	23	+import numpy
	24	+import libtetrabz_common
	25	+
	26	+
	27	+def dbldelta(bvec, eig1, eig2):
	28	+ """
	29	+ Main SUBROUTINE for Delta(E1) * Delta(E2)
	30	+ :param bvec:
	31	+ :param eig1:
	32	+ :param eig2:
	33	+ :return:
	34	+ """
	35	+ thr = 1.0e-10
	36	+ ng = numpy.array(eig1.shape[0:3])
	37	+ nk = ng.prod(0)
	38	+ nb = eig1.shape[3]
	39	+ wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
	40	+
	41	+ wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)
	42	+
	43	+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
	44	+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
	45	+
	46	+ for it in range(6 * nk):
	47	+ #
	48	+ for ii in range(20):
	49	+ eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	50	+ eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	51	+
	52	+ ei1 = wlsm.dot(eig1t)
	53	+ ej1 = wlsm.dot(eig2t)
	54	+ for ib in range(nb):
	55	+ #
	56	+ w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
	57	+ e = ei1[0:4, ib]
	58	+ indx = e.argsort(0)
	59	+ e.sort(0)
	60	+ #
	61	+ if e[0] < 0.0 <= e[1]:
	62	+ #
	63	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_a1(e)
	64	+ #
	65	+ if v > thr:
	66	+ #
	67	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	68	+ w2 = libtetrabz_dbldelta2(ej2)
	69	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	70	+ #
	71	+ elif e[1] < 0.0 <= e[2]:
	72	+ #
	73	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_b1(e)
	74	+ #
	75	+ if v > thr:
	76	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	77	+ w2 = libtetrabz_dbldelta2(ej2)
	78	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	79	+ #
	80	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e)
	81	+ #
	82	+ if v > thr:
	83	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	84	+ w2 = libtetrabz_dbldelta2(ej2)
	85	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	86	+ #
	87	+ elif e[2] < 0.0 < e[3]:
	88	+ #
	89	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_c1(e)
	90	+ #
	91	+ if v > thr:
	92	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	93	+ w2 = libtetrabz_dbldelta2(ej2)
	94	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	95	+ #
	96	+ else:
	97	+ continue
	98	+ #
	99	+ for ii in range(20):
	100	+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb] += w1.dot(wlsm)
	101	+ #
	102	+ wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
	103	+
	104	+
	105	+def libtetrabz_dbldelta2(ej):
	106	+ """
	107	+ 2nd step of tetrahedra method.
	108	+ :param ej:
	109	+ :return:
	110	+ """
	111	+ nb = ej.shape[1]
	112	+ w = numpy.zeros([nb, 3], dtype=numpy.float_)
	113	+ a = numpy.empty([3, 3], dtype=numpy.float_)
	114	+ for ib in range(nb):
	115	+ #
	116	+ if abs(ej[0:3, ib]).max() < 1.0e-10:
	117	+ raise ValueError("Nesting ##")
	118	+ #
	119	+ e = ej[0:3, ib]
	120	+ indx = e.argsort(0)
	121	+ e.sort(0)
	122	+ #
	123	+ for ii in range(3):
	124	+ a[0:3, ii] = (0.0 - e[ii]) / (e[0:3] - e[ii])
	125	+ #
	126	+ if e[0] < 0.0 <= e[1] or e[0] <= 0.0 < e[1]:
	127	+ #
	128	+ # v = a[1, 0] * a[2, 0] / (0.0 - e[0])
	129	+ v = a[1, 0] / (e[2] - e[0])
	130	+ #
	131	+ w[ib, indx[0]] = v * (a[0, 1] + a[0, 2])
	132	+ w[ib, indx[1]] = v * a[1, 0]
	133	+ w[ib, indx[2]] = v * a[2, 0]
	134	+ #
	135	+ elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
	136	+ #
	137	+ # v = a[1, 2] * a[1, 2] / (e[2] - 0.0)
	138	+ v = a[1, 2] / (e[2] - e[0])
	139	+ #
	140	+ w[ib, indx[0]] = v * a[0, 2]
	141	+ w[ib, indx[1]] = v * a[1, 2]
	142	+ w[ib, indx[2]] = v * (a[2, 0] + a[2, 1])
	143	+ #
	144	+ return w

--- /dev/null

+++ b/src/libtetrabz/libtetrabz_dblstep.py

		@@ -0,0 +1,203 @@
	1	+#
	2	+# Copyright (c) 2014 Mitsuaki Kawamura
	3	+#
	4	+# Permission is hereby granted, free of charge, to any person obtaining a
	5	+# copy of this software and associated documentation files (the
	6	+# "Software"), to deal in the Software without restriction, including
	7	+# without limitation the rights to use, copy, modify, merge, publish,
	8	+# distribute, sublicense, and/or sell copies of the Software, and to
	9	+# permit persons to whom the Software is furnished to do so, subject to
	10	+# the following conditions:
	11	+#
	12	+# The above copyright notice and this permission notice shall be included
	13	+# in all copies or substantial portions of the Software.
	14	+#
	15	+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	16	+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
	17	+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
	18	+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
	19	+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
	20	+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
	21	+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
	22	+#
	23	+import numpy
	24	+import libtetrabz_common
	25	+
	26	+
	27	+def dblstep(bvec, eig1, eig2):
	28	+ """
	29	+ Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
	30	+ :param bvec:
	31	+ :param eig1:
	32	+ :param eig2:
	33	+ :return:
	34	+ """
	35	+ ng = numpy.array(eig1.shape[0:3])
	36	+ nk = ng.prod(0)
	37	+ nb = eig1.shape[3]
	38	+ wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
	39	+
	40	+ wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)
	41	+
	42	+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
	43	+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
	44	+
	45	+ thr = 1.0e-8
	46	+ #
	47	+ for it in range(6 * nk):
	48	+ #
	49	+ for ii in range(20):
	50	+ eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	51	+ eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	52	+
	53	+ ei1 = wlsm.dot(eig1t)
	54	+ ej1 = wlsm.dot(eig2t)
	55	+ #
	56	+ for ib in range(nb):
	57	+ #
	58	+ w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
	59	+ e = ei1[0:4, ib]
	60	+ indx = e.argsort(0)
	61	+ e.sort(0)
	62	+ #
	63	+ if e(1) <= 0.0 < e(2):
	64	+ #
	65	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
	66	+ #
	67	+ if v > thr:
	68	+ #
	69	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	70	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	71	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	72	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	73	+ #
	74	+ elif e(2) <= 0.0 < e(3):
	75	+ #
	76	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
	77	+ #
	78	+ if v > thr:
	79	+ #
	80	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	81	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	82	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	83	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	84	+ #
	85	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
	86	+ #
	87	+ if v > thr:
	88	+ #
	89	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	90	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	91	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	92	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	93	+ #
	94	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
	95	+ #
	96	+ if v > thr:
	97	+ #
	98	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	99	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	100	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	101	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	102	+ #
	103	+ elif e(3) <= 0.0 < e(4):
	104	+ #
	105	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
	106	+ #
	107	+ if v > thr:
	108	+ #
	109	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	110	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	111	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	112	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	113	+ #
	114	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
	115	+ #
	116	+ if v > thr:
	117	+ #
	118	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	119	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	120	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	121	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	122	+ #
	123	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
	124	+ #
	125	+ if v > thr:
	126	+ #
	127	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	128	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	129	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	130	+ w1[0:nb, indx[0:4]] += v * w2.dot(tsmall)
	131	+ #
	132	+ elif e(4) <= 0.0:
	133	+ #
	134	+ ei2 = ei1[0:4, ib]
	135	+ ej2 = ej1[0:4, 0:nb]
	136	+ w2 = libtetrabz_dblstep2(ei2, ej2)
	137	+ w1[0:nb, 0:4] += w2[0:nb, 0:4]
	138	+ #
	139	+ else:
	140	+ continue
	141	+ #
	142	+ for ii in range(20):
	143	+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb] += w1.dot(wlsm)
	144	+
	145	+ wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb] /= (6.0 * nk)
	146	+
	147	+
	148	+def libtetrabz_dblstep2(ei1, ej1):
	149	+ """
	150	+ Tetrahedra method for theta( - de)
	151	+ :param ei1:
	152	+ :param ej1:
	153	+ :return:
	154	+ """
	155	+ thr = 1.0e-8
	156	+ nb = ej1.shape[1]
	157	+ w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
	158	+
	159	+ for ib in range(nb):
	160	+ #
	161	+ e = - ei1[0:4] + ej1[0:4, ib]
	162	+ indx = e.argsort(0)
	163	+ e.sort(0)
	164	+ #
	165	+ if abs(e(1)) < thr and abs(e(4)) < thr:
	166	+ #
	167	+ # Theta(0) = 0.5
	168	+ #
	169	+ v = 0.5
	170	+ w1[ib, 0:4] += v * 0.25
	171	+ #
	172	+ elif e(1) <= 0.0 < e(2) or e(1) < 0.0 <= e(2):
	173	+ #
	174	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
	175	+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
	176	+ #
	177	+ elif e(2) <= 0.0 < e(3) or e(2) < 0.0 <= e(3):
	178	+ #
	179	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
	180	+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
	181	+ #
	182	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
	183	+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
	184	+ #
	185	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
	186	+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
	187	+ #
	188	+ elif e(3) <= 0.0 < e(4) or e(3) < 0.0 <= e(4):
	189	+ #
	190	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
	191	+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
	192	+ #
	193	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
	194	+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
	195	+ #
	196	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
	197	+ w1[ib, indx[0:4]] += v * tsmall.sum(0) * 0.25
	198	+ #
	199	+ elif e(4) <= 0.0:
	200	+ #
	201	+ w1[ib, 0:4] += 0.25
	202	+ #
	203	+ return w1

--- a/src/libtetrabz/libtetrabz_dos_mod.py

+++ b/src/libtetrabz/libtetrabz_dos.py

		@@ -24,7 +24,7 @@ import numpy
24	24	import libtetrabz_common
25	25
26	26
27		-def libtetrabz_dos(bvec, eig, e0):
	27	+def dos(bvec, eig, e0):
28	28	"""
29	29	Compute Dos : Delta(E - E1)
30	30	:param bvec:

		@@ -87,7 +87,7 @@ def libtetrabz_dos(bvec, eig, e0):
87	87	return wght
88	88
89	89
90		-def libtetrabz_intdos(bvec, eig, e0):
	90	+def intdos(bvec, eig, e0):
91	91	"""
92	92	Compute integrated Dos : theta(E - E1)
93	93	:param bvec:

--- /dev/null

+++ b/src/libtetrabz/libtetrabz_fermigr.py

		@@ -0,0 +1,275 @@
	1	+#
	2	+# Copyright (c) 2014 Mitsuaki Kawamura
	3	+#
	4	+# Permission is hereby granted, free of charge, to any person obtaining a
	5	+# copy of this software and associated documentation files (the
	6	+# "Software"), to deal in the Software without restriction, including
	7	+# without limitation the rights to use, copy, modify, merge, publish,
	8	+# distribute, sublicense, and/or sell copies of the Software, and to
	9	+# permit persons to whom the Software is furnished to do so, subject to
	10	+# the following conditions:
	11	+#
	12	+# The above copyright notice and this permission notice shall be included
	13	+# in all copies or substantial portions of the Software.
	14	+#
	15	+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	16	+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
	17	+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
	18	+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
	19	+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
	20	+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
	21	+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
	22	+#
	23	+import numpy
	24	+import libtetrabz_common
	25	+
	26	+
	27	+def fermigr(bvec, eig1, eig2, e0):
	28	+ """
	29	+ Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
	30	+ :param bvec:
	31	+ :param eig1:
	32	+ :param eig2:
	33	+ :param e0:
	34	+ :return:
	35	+ """
	36	+ ng = numpy.array(eig1.shape[0:3])
	37	+ nk = ng.prod(0)
	38	+ nb = eig1.shape[3]
	39	+ ne = e0.shape[0]
	40	+ wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
	41	+
	42	+ wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb, ne], dtype=numpy.float_)
	43	+
	44	+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
	45	+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
	46	+ thr = 1.0e-10
	47	+ #
	48	+ for it in range(6 * nk):
	49	+ #
	50	+ for ii in range(20):
	51	+ eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	52	+ eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	53	+
	54	+ ei1 = wlsm.dot(eig1t)
	55	+ ej1 = wlsm.dot(eig2t)
	56	+ #
	57	+ for ib in range(nb):
	58	+ #
	59	+ w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
	60	+ e = ei1[0:4, ib]
	61	+ indx = e.argsort(0)
	62	+ e.sort(0)
	63	+ #
	64	+ if e[0] <= 0.0 < e[1]:
	65	+ #
	66	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
	67	+ #
	68	+ if v > thr:
	69	+ #
	70	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	71	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	72	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	73	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	74	+ #
	75	+ elif e[1] <= 0.0 < e[2]:
	76	+ #
	77	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
	78	+ #
	79	+ if v > thr:
	80	+ #
	81	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	82	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	83	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	84	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	85	+ #
	86	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
	87	+ #
	88	+ if v > thr:
	89	+ #
	90	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	91	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	92	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	93	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	94	+ #
	95	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
	96	+ #
	97	+ if v > thr:
	98	+ #
	99	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	100	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	101	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	102	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	103	+ #
	104	+ elif e[2] <= 0.0 < e[3]:
	105	+ #
	106	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
	107	+ #
	108	+ if v > thr:
	109	+ #
	110	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	111	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	112	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	113	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	114	+ #
	115	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
	116	+ #
	117	+ if v > thr:
	118	+ #
	119	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	120	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	121	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	122	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	123	+ #
	124	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
	125	+ #
	126	+ if v > thr:
	127	+ #
	128	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	129	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	130	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	131	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	132	+ #
	133	+ elif e[3] <= 0.0:
	134	+ #
	135	+ ei2 = ei1[0:4, ib]
	136	+ ej2 = ej1[0:4, 0:nb]
	137	+ w2 = libtetrabz_fermigr2(e0, ei2, ej2)
	138	+ w1[0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
	139	+ #
	140	+ else:
	141	+ continue
	142	+ #
	143	+ for ii in range(20):
	144	+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb, 0:ne] += w1.dot(wlsm)
	145	+ #
	146	+ wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb, 0:ne] /= (6.0 * nk)
	147	+
	148	+
	149	+def libtetrabz_fermigr2(e0, ei1, ej1):
	150	+ """
	151	+ Tetrahedra method for theta( - E2)
	152	+ :param e0:
	153	+ :param ei1:
	154	+ :param ej1:
	155	+ :return:
	156	+ """
	157	+ #
	158	+ ne = e0.shape[0]
	159	+ nb = ej1.shape[1]
	160	+ w1 = numpy.zeros([nb, ne, 4], dtype=numpy.float_)
	161	+ thr = 1.0e-8
	162	+ #
	163	+ for ib in range(nb):
	164	+ #
	165	+ e = -ej1[0:4, ib]
	166	+ indx = e.argsort(0)
	167	+ e.sort(0)
	168	+ #
	169	+ if e[0] <= 0.0 < e[1] or e[0] < 0.0 <= e[1]:
	170	+ #
	171	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
	172	+ #
	173	+ if v > thr:
	174	+ #
	175	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	176	+ w2 = libtetrabz_fermigr3(e0, de)
	177	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	178	+ #
	179	+ #
	180	+ elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
	181	+ #
	182	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
	183	+ #
	184	+ if v > thr:
	185	+ #
	186	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	187	+ w2 = libtetrabz_fermigr3(e0, de)
	188	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	189	+ #
	190	+ #
	191	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
	192	+ #
	193	+ if v > thr:
	194	+ #
	195	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	196	+ w2 = libtetrabz_fermigr3(e0, de)
	197	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	198	+ #
	199	+ #
	200	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
	201	+ #
	202	+ if v > thr:
	203	+ #
	204	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	205	+ w2 = libtetrabz_fermigr3(e0, de)
	206	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	207	+ #
	208	+ #
	209	+ elif e[2] <= 0.0 < e[3] or e[2] < 0.0 <= e[3]:
	210	+ #
	211	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
	212	+ #
	213	+ if v > thr:
	214	+ #
	215	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	216	+ w2 = libtetrabz_fermigr3(e0, de)
	217	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	218	+ #
	219	+ #
	220	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
	221	+ #
	222	+ if v > thr:
	223	+ #
	224	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	225	+ w2 = libtetrabz_fermigr3(e0, de)
	226	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	227	+ #
	228	+ #
	229	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
	230	+ #
	231	+ if v > thr:
	232	+ #
	233	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	234	+ w2 = libtetrabz_fermigr3(e0, de)
	235	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	236	+ #
	237	+ #
	238	+ elif e[3] <= 0.0:
	239	+ #
	240	+ de = ej1[0:4, ib] - ei1[0:4]
	241	+ w2 = libtetrabz_fermigr3(e0, de)
	242	+ w1[ib, 0:ne, 0:4] += w2
	243	+ #
	244	+ return w1
	245	+
	246	+
	247	+def libtetrabz_fermigr3(e0, de):
	248	+ #
	249	+ ne = e0.shape[0]
	250	+ e = de[0:4]
	251	+ indx = e.argsort(0)
	252	+ e.sort(0)
	253	+ w1 = numpy.zeros([ne, 4], dtype=numpy.float_)
	254	+ #
	255	+ for ie in range(ne):
	256	+ #
	257	+ if e[0] < e0[ie] .AND. e0[ie] <= e[1]:
	258	+ #
	259	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_a1(e[0:4] - e0[ie])
	260	+ w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
	261	+ #
	262	+ elif e[1] < e0[ie] .AND. e0[ie] <= e[2]:
	263	+ #
	264	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_b1(e[0:4] - e0[ie])
	265	+ w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
	266	+ #
	267	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_b2(e[0:4] - e0[ie])
	268	+ w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
	269	+ #
	270	+ elif e[2] < e0[ie] .AND. e0[ie] < e[3]:
	271	+ #
	272	+ v, tsmall = libtetrabz_common.libtetrabz_triangle_c1(e[0:4] - e0[ie])
	273	+ w1[ie, indx[0:4]] += v * tsmall.sum(0) / 3.0
	274	+ #
	275	+ return w1

--- a/src/libtetrabz/libtetrabz_occ_mod.py

+++ b/src/libtetrabz/libtetrabz_occ.py

		@@ -24,7 +24,7 @@ import numpy
24	24	import libtetrabz_common
25	25
26	26
27		-def libtetrabz_fermieng(bvec, eig, nelec):
	27	+def fermieng(bvec, eig, nelec):
28	28	"""
29	29	Calculate Fermi energy
30	30	:param bvec:

		@@ -47,7 +47,7 @@ def libtetrabz_fermieng(bvec, eig, nelec):
47	47	#
48	48	# Calc. # of electrons
49	49	#
50		- wght = libtetrabz_occ(bvec, eig-ef)
	50	+ wght = occ(bvec, eig-ef)
51	51	sumkmid = wght.sum()
52	52	#
53	53	# convergence check

		@@ -62,7 +62,7 @@ def libtetrabz_fermieng(bvec, eig, nelec):
62	62	raise ValueError("libtetrabz_fermieng")
63	63
64	64
65		-def libtetrabz_occ(bvec, eig):
	65	+def occ(bvec, eig):
66	66	"""
67	67	Main SUBROUTINE for occupation : Theta(EF - E1)
68	68	:param bvec:

--- /dev/null

+++ b/src/libtetrabz/libtetrabz_polcmplx.py

		@@ -0,0 +1,622 @@
	1	+#
	2	+# Copyright (c) 2014 Mitsuaki Kawamura
	3	+#
	4	+# Permission is hereby granted, free of charge, to any person obtaining a
	5	+# copy of this software and associated documentation files (the
	6	+# "Software"), to deal in the Software without restriction, including
	7	+# without limitation the rights to use, copy, modify, merge, publish,
	8	+# distribute, sublicense, and/or sell copies of the Software, and to
	9	+# permit persons to whom the Software is furnished to do so, subject to
	10	+# the following conditions:
	11	+#
	12	+# The above copyright notice and this permission notice shall be included
	13	+# in all copies or substantial portions of the Software.
	14	+#
	15	+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	16	+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
	17	+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
	18	+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
	19	+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
	20	+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
	21	+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
	22	+#
	23	+import math
	24	+import numpy
	25	+import libtetrabz_common
	26	+
	27	+
	28	+def polcmplx(bvec, eig1, eig2, e0):
	29	+ """
	30	+ Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
	31	+ :param bvec:
	32	+ :param eig1:
	33	+ :param eig2:
	34	+ :param e0:
	35	+ :return:
	36	+ """
	37	+ ng = numpy.array(eig1.shape[0:3])
	38	+ nk = ng.prod(0)
	39	+ nb = eig1.shape[3]
	40	+ ne = e0.shape[0]
	41	+ wlsm, ikv = libtetrabz_common.libtetrabz_initialize(ng, bvec)
	42	+
	43	+ wght = numpy.zeros([ng[0], ng[1], ng[2], nb, nb, ne], dtype=numpy.complex_)
	44	+
	45	+ eig1t = numpy.empty([20, nb], dtype=numpy.float_)
	46	+ eig2t = numpy.empty([20, nb], dtype=numpy.float_)
	47	+
	48	+ #
	49	+ thr = 1.0e-8
	50	+ for it in range(6 * nk):
	51	+ #
	52	+ for ii in range(20):
	53	+ eig1t[ii, 0:nb] = eig1[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	54	+ eig2t[ii, 0:nb] = eig2[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], 0:nb]
	55	+
	56	+ ei1 = wlsm.dot(eig1t)
	57	+ ej1 = wlsm.dot(eig2t)
	58	+ #
	59	+ for ib in range(nb):
	60	+ #
	61	+ w1 = numpy.zeros([nb, 4], dtype=numpy.float_)
	62	+ e = ei1[0:4, ib]
	63	+ indx = e.argsort(0)
	64	+ e.sort(0)
	65	+ #
	66	+ if e[0] <= 0.0 < e(2):
	67	+ #
	68	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
	69	+ #
	70	+ if v > thr:
	71	+ #
	72	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	73	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	74	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	75	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	76	+ #
	77	+ #
	78	+ elif e(2) <= 0.0 < e[2]:
	79	+ #
	80	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
	81	+ #
	82	+ if v > thr:
	83	+ #
	84	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	85	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	86	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	87	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	88	+ #
	89	+ #
	90	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
	91	+ #
	92	+ if v > thr:
	93	+ #
	94	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	95	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	96	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	97	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	98	+ #
	99	+ #
	100	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
	101	+ #
	102	+ if v > thr:
	103	+ #
	104	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	105	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	106	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	107	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	108	+ #
	109	+ #
	110	+ elif e[2] <= 0.0 < e[3]:
	111	+ #
	112	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
	113	+ #
	114	+ if v > thr:
	115	+ #
	116	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	117	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	118	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	119	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	120	+ #
	121	+ #
	122	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
	123	+ #
	124	+ if v > thr:
	125	+ #
	126	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	127	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	128	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	129	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	130	+ #
	131	+ #
	132	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
	133	+ #
	134	+ if v > thr:
	135	+ #
	136	+ ei2 = tsmall.dot(ei1[indx[0:4], ib])
	137	+ ej2 = tsmall.dot(ej1[indx[0:4], 0:nb])
	138	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	139	+ w1[0:nb, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	140	+ #
	141	+ #
	142	+ elif e[3] <= 0.0:
	143	+ #
	144	+ ei2 = ei1[0:4, ib]
	145	+ ej2 = ej1[0:4, 0:nb]
	146	+ w2 = libtetrabz_polcmplx2(e0, ei2, ej2)
	147	+ w1[0:nb, 0:ne, 0:4] += w2[0:nb, 0:ne, 0:4]
	148	+ #
	149	+ else:
	150	+ continue
	151	+ #
	152	+ for ii in range(20):
	153	+ wght[ikv[it, ii, 0], ikv[it, ii, 1], ikv[it, ii, 2], ib, 0:nb, 0:ne] += w1.dot(wlsm)
	154	+ #
	155	+ wght[0:ng[0], 0:ng[1], 0:ng[2], 0:nb, 0:nb, 0:ne] /= (6.0 * nk)
	156	+
	157	+
	158	+def libtetrabz_polcmplx2(e0, ei1, ej1):
	159	+ """
	160	+ Tetrahedra method for theta( - E2)
	161	+ :param e0:
	162	+ :param ei1:
	163	+ :param ej1:
	164	+ :return:
	165	+ """
	166	+ ne = e0.shape[0]
	167	+ nb = ej1.shape[1]
	168	+ thr = 1.0e-8
	169	+ w1 = numpy.zeros([nb, ne, 4], dtype=numpy.complex_)
	170	+
	171	+ for ib in range(nb):
	172	+ #
	173	+ e = -ej1[0:4, ib]
	174	+ indx = e.argsort(0)
	175	+ e.sort(0)
	176	+ #
	177	+ if e[0] <= 0.0 < e(2) or e[0] < 0.0 <= e[1]:
	178	+ #
	179	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_a1(e)
	180	+ #
	181	+ if v > thr:
	182	+ #
	183	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	184	+ w2 = libtetrabz_polcmplx3(e0, de)
	185	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	186	+ #
	187	+ #
	188	+ elif e[1] <= 0.0 < e[2] or e[1] < 0.0 <= e[2]:
	189	+ #
	190	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b1(e)
	191	+ #
	192	+ if v > thr:
	193	+ #
	194	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	195	+ w2 = libtetrabz_polcmplx3(e0, de)
	196	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	197	+ #
	198	+ #
	199	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b2(e)
	200	+ #
	201	+ if v > thr:
	202	+ #
	203	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	204	+ w2 = libtetrabz_polcmplx3(e0, de)
	205	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	206	+ #
	207	+ #
	208	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_b3(e)
	209	+ #
	210	+ if v > thr:
	211	+ #
	212	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	213	+ w2 = libtetrabz_polcmplx3(e0, de)
	214	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	215	+ #
	216	+ #
	217	+ elif e[2] <= 0.0 < e[3] or e[2] < 0.0 <= e[3]:
	218	+ #
	219	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c1(e)
	220	+ #
	221	+ if v > thr:
	222	+ #
	223	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	224	+ w2 = libtetrabz_polcmplx3(e0, de)
	225	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	226	+ #
	227	+ #
	228	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c2(e)
	229	+ #
	230	+ if v > thr:
	231	+ #
	232	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	233	+ w2 = libtetrabz_polcmplx3(e0, de)
	234	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	235	+ #
	236	+ #
	237	+ v, tsmall = libtetrabz_common.libtetrabz_tsmall_c3(e)
	238	+ #
	239	+ if v > thr:
	240	+ #
	241	+ de = tsmall.dot(ej1[indx[0:4], ib] - ei1[indx[0:4]])
	242	+ w2 = libtetrabz_polcmplx3(e0, de)
	243	+ w1[ib, 0:ne, indx[0:4]] += v * w2.dot(tsmall)
	244	+ #
	245	+ #
	246	+ elif e[3] <= 0.0:
	247	+ #
	248	+ de = ej1[0:4, ib] - ei1[0:4]
	249	+ w2 = libtetrabz_polcmplx3(e0, de)
	250	+ w1[ib, 0:ne, 0:4] += w2
	251	+
	252	+ return w1
	253	+
	254	+
	255	+def libtetrabz_polcmplx3(e0, de):
	256	+ """
	257	+ Tetrahedron method for delta(om - ep + e)
	258	+ :param e0:
	259	+ :param de:
	260	+ :return:
	261	+ """
	262	+ ne = e0.shape[0]
	263	+ e = de[0:4]
	264	+ indx = e.argsort(0)
	265	+ e.sort(0)
	266	+ w1 = numpy.zeros([ne, 4], dtype=numpy.complex_)
	267	+ w2 = numpy.empty([2, 4], dtype=numpy.float_)
	268	+ #
	269	+ for ie in range(ne):
	270	+ #
	271	+ # I don't know which one is better.
	272	+ # The former is more stable.
	273	+ # The latter is more accurate ?
	274	+ #
	275	+ w1[ie, 0:4] = 0.25 / (de[0:4] + e0[ie])
	276	+ #
	277	+ continue
	278	+ #
	279	+ x = (e[0:4] + e0[ie].real) / e0[ie].imag
	280	+ # thr = maxval(de(1:4)) * 1d-3
	281	+ thr = max(1.0e-3, x.max() * 1.0e-2)
	282	+ #
	283	+ if abs(x[3] - x(3)) < thr:
	284	+ if abs(x[3] - x(2)) < thr:
	285	+ if abs(x[3] - x(1)) < thr:
	286	+ #
	287	+ # e[3] = e[2] = e[1] = e[0]
	288	+ #
	289	+ w2[0, 3] = 0.25 * x[3] / (1.0 + x[3]**2)
	290	+ w2[1, 3] = 0.25 / (1.0 + x[3]**2)
	291	+ w2[0:2, 2] = w2[0:2, 3]
	292	+ w2[0:2, 1] = w2[0:2, 3]
	293	+ w2[0:2, 0] = w2[0:2, 3]
	294	+ #
	295	+ else:
	296	+ #
	297	+ # e[3] = e[2] = e[1]
	298	+ #
	299	+ w2[0:2, 3] = libtetrabz_polcmplx_1211(x[3], x(1))
	300	+ w2[0:2, 2] = w2[0:2, 3]
	301	+ w2[0:2, 1] = w2[0:2, 3]
	302	+ w2[0:2, 0] = libtetrabz_polcmplx_1222(x[0], x(4))
	303	+ #
	304	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	305	+ # WRITE(,) ie
	306	+ # WRITE(*,'(100e15.5)') x[0:4]
	307	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	308	+ # STOP "weighting 4=3=2"
	309	+ #
	310	+ elif abs(x[1] - x(1)) < thr:
	311	+ #
	312	+ # e[3] = e[2], e[1] = e[0]
	313	+ #
	314	+ w2[0:2, 3] = libtetrabz_polcmplx_1221(x[3], x(2))
	315	+ w2[0:2, 2] = w2[0:2, 3]
	316	+ w2[0:2, 1] = libtetrabz_polcmplx_1221(x[1], x(4))
	317	+ w2[0:2, 0] = w2[0:2, 1]
	318	+ #
	319	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	320	+ # WRITE(,) ie
	321	+ # WRITE(*,'(100e15.5)') x[0:4]
	322	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	323	+ # STOP "weighting 4=3 2=1"
	324	+ #
	325	+ else:
	326	+ #
	327	+ # e[3] = e[2]
	328	+ #
	329	+ w2[0:2, 3] = libtetrabz_polcmplx_1231(x[3], x[0], x(2))
	330	+ w2[0:2, 2] = w2[0:2, 3]
	331	+ w2[0:2, 1] = libtetrabz_polcmplx_1233(x[1], x[0], x(4))
	332	+ w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[1], x(4))
	333	+ #
	334	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	335	+ # WRITE(,) ie
	336	+ # WRITE(*,'(100e15.5)') x[0:4]
	337	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	338	+ # STOP "weighting 4=3"
	339	+ #
	340	+ elif abs(x[2] - x(2)) < thr:
	341	+ if abs(x[2] - x(1)) < thr:
	342	+ #
	343	+ # e[2] = e[1] = e[0]
	344	+ #
	345	+ w2[0:2, 3] = libtetrabz_polcmplx_1222(x[3], x(3))
	346	+ w2[0:2, 2] = libtetrabz_polcmplx_1211(x[2], x(4))
	347	+ w2[0:2, 1] = w2[0:2, 2]
	348	+ w2[0:2, 0] = w2[0:2, 2]
	349	+ #
	350	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	351	+ # WRITE(,) ie
	352	+ # WRITE(*,'(100e15.5)') x[0:4]
	353	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	354	+ # STOP "weighting 3=2=1"
	355	+ #
	356	+ else:
	357	+ #
	358	+ # e[2] = e[1]
	359	+ #
	360	+ w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[0], x(3))
	361	+ w2[0:2, 2] = libtetrabz_polcmplx_1231(x[2], x[0], x(4))
	362	+ w2[0:2, 1] = w2[0:2, 2]
	363	+ w2[0:2, 0] = libtetrabz_polcmplx_1233(x[0], x[3], x(3))
	364	+ #
	365	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	366	+ # WRITE(,) ie
	367	+ # WRITE(*,'(100e15.5)') x[0:4]
	368	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	369	+ # STOP "weighting 3=2"
	370	+ #
	371	+ elif abs(x[1] - x(1)) < thr:
	372	+ #
	373	+ # e[1] = e[0]
	374	+ #
	375	+ w2[0:2, 3] = libtetrabz_polcmplx_1233(x[3], x[2], x(2))
	376	+ w2[0:2, 2] = libtetrabz_polcmplx_1233(x[2], x[3], x(2))
	377	+ w2[0:2, 1] = libtetrabz_polcmplx_1231(x[1], x[2], x(4))
	378	+ w2[0:2, 0] = w2[0:2, 1]
	379	+ #
	380	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	381	+ # WRITE(,) ie
	382	+ # WRITE(*,'(100e15.5)') x[0:4]
	383	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	384	+ # STOP "weighting 2=1"
	385	+ #
	386	+ else:
	387	+ #
	388	+ # Different each other.
	389	+ #
	390	+ w2[0:2, 3] = libtetrabz_polcmplx_1234(x[3], x[0], x[1], x(3))
	391	+ w2[0:2, 2] = libtetrabz_polcmplx_1234(x[2], x[0], x[1], x(4))
	392	+ w2[0:2, 1] = libtetrabz_polcmplx_1234(x[1], x[0], x[2], x(4))
	393	+ w2[0:2, 0] = libtetrabz_polcmplx_1234(x[0], x[1], x[2], x(4))
	394	+ #
	395	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	396	+ # WRITE(,) ie
	397	+ # WRITE(*,'(100e15.5)') x[0:4]
	398	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	399	+ # STOP "weighting"
	400	+ #
	401	+ #
	402	+ w1[ie, indx[0:4]] = complex(w2[0, 0:4] / e0[ie].imag, w2[1, 0:4] / (- e0[ie].imag))
	403	+ #
	404	+ return w1
	405	+#
	406	+# Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
	407	+# for 0<x<1, 0<y<1-x, 0<z<1-x-y
	408	+
	409	+
	410	+def libtetrabz_polcmplx_1234(g1, g2, g3, g4):
	411	+ """
	412	+ 1, Different each other
	413	+ :param g1:
	414	+ :param g2:
	415	+ :param g3:
	416	+ :param g4:
	417	+ :return:
	418	+ """
	419	+ w = numpy.empty(2, dtype=numpy.float_)
	420	+ #
	421	+ # Real
	422	+ #
	423	+ w2 = 2.0(3.0g2*2 - 1.0)(math.atan(g2) - math.atan(g1)) + (g2*2 - 3.0)g2math.log((1.0 + g22)/(1.0 + g1*2))
	424	+ w2 = -2.0(g2*2 - 1.0) + w2/(g2 - g1)
	425	+ w2 = w2/(g2 - g1)
	426	+ w3 = 2.0(3.0g3*2 - 1.0)(math.atan(g3) - math.atan(g1)) + (g3*2 - 3.0)g3math.log((1.0 + g32)/(1.0 + g1*2))
	427	+ w3 = -2.0(g3*2 - 1.0) + w3/(g3 - g1)
	428	+ w3 = w3/(g3 - g1)
	429	+ w4 = 2.0(3.0g4*2 - 1.0)(math.atan(g4) - math.atan(g1)) + (g4*2 - 3.0)g4math.log((1.0 + g42)/(1.0 + g1*2))
	430	+ w4 = -2.0(g4*2 - 1.0) + w4/(g4 - g1)
	431	+ w4 = w4/(g4 - g1)
	432	+ w2 = (w2 - w3)/(g2 - g3)
	433	+ w4 = (w4 - w3)/(g4 - g3)
	434	+ w[0] = (w4 - w2)/(2.0*(g4 - g2))
	435	+ #
	436	+ # Imaginal
	437	+ #
	438	+ w2 = 2.0(3.0 - g22)g2(math.atan(g2) - math.atan(g1)) + (3.0g2*2 - 1.0)math.log((1.0 + g22)/(1.0 + g12))
	439	+ w2 = 4.0*g2 - w2/(g2 - g1)
	440	+ w2 = w2/(g2 - g1)
	441	+ w3 = 2.0(3.0 - g32)g3(math.atan(g3) - math.atan(g1)) + (3.0g3*2 - 1.0)math.log((1.0 + g32)/(1.0 + g12))
	442	+ w3 = 4.0*g3 - w3/(g3 - g1)
	443	+ w3 = w3/(g3 - g1)
	444	+ w4 = 2.0(3.0 - g42)g4(math.atan(g4) - math.atan(g1)) + (3.0g4*2 - 1.0)math.log((1.0 + g42)/(1.0 + g12))
	445	+ w4 = 4.0*g4 - w4/(g4 - g1)
	446	+ w4 = w4/(g4 - g1)
	447	+ w2 = (w2 - w3)/(g2 - g3)
	448	+ w4 = (w4 - w3)/(g4 - g3)
	449	+ w[1] = (w4 - w2)/(2.0*(g4 - g2))
	450	+ #
	451	+ return w
	452	+
	453	+
	454	+def libtetrabz_polcmplx_1231(g1, g2, g3):
	455	+ """
	456	+ 2, g4 = g1
	457	+ :param g1:
	458	+ :param g2:
	459	+ :param g3:
	460	+ :return:
	461	+ """
	462	+ w = numpy.empty(2, dtype=numpy.float_)
	463	+ #
	464	+ # Real
	465	+ #
	466	+ w2 = 2.0(-1.0 + 3.0g2*2)(math.atan(g2) - math.atan(g1))\
	467	+ + g2(-3.0 + g22)math.log((1.0 + g22)/(1.0 + g12))
	468	+ w2 = 2.0(1.0 - g2*2) + w2/(g2 - g1)
	469	+ w2 = -g1 + w2/(g2 - g1)
	470	+ w2 = w2/(g2 - g1)
	471	+ w3 = 2.0(-1.0 + 3.0g3*2)(math.atan(g3) - math.atan(g1))\
	472	+ + g3(-3.0 + g32)math.log((1.0 + g32)/(1.0 + g12))
	473	+ w3 = 2.0(1 - g3*2) + w3/(g3 - g1)
	474	+ w3 = -g1 + w3/(g3 - g1)
	475	+ w3 = w3/(g3 - g1)
	476	+ w[0] = (w3 - w2)/(2.0*(g3 - g2))
	477	+ #
	478	+ # Imaginal
	479	+ #
	480	+ w2 = 2.0g2(3.0 - g2*2)(math.atan(g2) - math.atan(g1)) + (-1.0 + 3.0g22)math.log((1.0 + g22)/(1.0 + g12))
	481	+ w2 = 4.0*g2 - w2/(g2 - g1)
	482	+ w2 = 1 + w2/(g2 - g1)
	483	+ w2 = w2/(g2 - g1)
	484	+ w3 = 2.0g3(3.0 - g3*2)(math.atan(g3) - math.atan(g1)) + (-1.0 + 3.0g32)math.log((1.0 + g32)/(1.0 + g12))
	485	+ w3 = 4.0*g3 - w3/(g3 - g1)
	486	+ w3 = 1 + w3/(g3 - g1)
	487	+ w3 = w3/(g3 - g1)
	488	+ w[1] = (w3 - w2)/(2.0*(g3 - g2))
	489	+ #
	490	+ return w
	491	+
	492	+
	493	+def libtetrabz_polcmplx_1233(g1, g2, g3):
	494	+ """
	495	+ 3, g4 = g3
	496	+ :param g1:
	497	+ :param g2:
	498	+ :param g3:
	499	+ :return:
	500	+ """
	501	+ w = numpy.empty(2, dtype=numpy.float_)
	502	+ #
	503	+ # Real
	504	+ #
	505	+ w2 = 2.0(1.0 - 3.0g2*2)(math.atan(g2) - math.atan(g1)) + g2(3.0 - g22)math.log((1.0 + g22)/(1.0 + g12))
	506	+ w2 = 2.0(1 - g2*2) - w2/(g2 - g1)
	507	+ w2 = w2/(g2 - g1)
	508	+ w3 = 2.0(1.0 - 3.0g3*2)(math.atan(g3) - math.atan(g1)) + g3(3.0 - g32)math.log((1.0 + g32)/(1.0 + g12))
	509	+ w3 = 2.0(1 - g3*2) - w3/(g3 - g1)
	510	+ w3 = w3/(g3 - g1)
	511	+ w2 = (w3 - w2)/(g3 - g2)
	512	+ w3 = 4.0(1.0 - 3.0g1g3)(math.atan(g3) - math.atan(g1))\
	513	+ + (3.0g1 + 3.0g3 - 3.0g1g32 + g33) * math.log((1.0 + g32)/(1.0 + g12))
	514	+ w3 = -4.0(1.0 - g1*2) + w3/(g3 - g1)
	515	+ w3 = 4.0*g1 + w3/(g3 - g1)
	516	+ w3 = w3/(g3 - g1)
	517	+ w[0] = (w3 - w2)/(2.0*(g3 - g2))
	518	+ #
	519	+ # Imaginal
	520	+ #
	521	+ w2 = 2.0g2(3.0 - g2*2)(math.atan(g2) - math.atan(g1)) + (-1.0 + 3.0g22)math.log((1.0 + g22)/(1.0 + g12))
	522	+ w2 = 4.0*g2 - w2/(g2 - g1)
	523	+ w2 = w2/(g2 - g1)
	524	+ w3 = 2.0g3(3.0 - g3*2)(math.atan(g3) - math.atan(g1)) + (-1.0 + 3.0g32)math.log((1.0 + g32)/(1.0 + g12))
	525	+ w3 = 4.0*g3 - w3/(g3 - g1)
	526	+ w3 = w3/(g3 - g1)
	527	+ w2 = (w3 - w2)/(g3 - g2)
	528	+ w3 = (3.0g1 - 3.0g1g32 + 3.0g3 + g3*3)(math.atan(g3) - math.atan(g1)) \
	529	+ + (3.0g1g3 - 1.0)math.log((1.0 + g32)/(1.0 + g1*2))
	530	+ w3 = w3/(g3 - g1) - 4.0*g1
	531	+ w3 = w3/(g3 - g1) - 2.0
	532	+ w3 = (2.0*w3)/(g3 - g1)
	533	+ w[1] = (w3 - w2)/(2.0*(g3 - g2))
	534	+ #
	535	+ return w
	536	+
	537	+
	538	+def libtetrabz_polcmplx_1221(g1, g2):
	539	+ """
	540	+ 4, g4 = g1 and g3 = g2
	541	+ :param g1:
	542	+ :param g2:
	543	+ :return:
	544	+ """
	545	+ w = numpy.empty(2, dtype=numpy.float_)
	546	+ #
	547	+ # Real
	548	+ #
	549	+ w[0] = -2.0(-1.0 + 2.0g1g2 + g22)(math.atan(g2) - math.atan(g1)) \
	550	+ + (g1 + 2.0g2 - g1g2*2)math.log((1.0 + g22)/(1.0 + g12))
	551	+ w[0] = 2.0(-1.0 + g1*2) + w[0]/(g2 - g1)
	552	+ w[0] = 3.0*g1 + w[0]/(g2 - g1)
	553	+ w[0] = 2.0 + (3.0*w[0])/(g2 - g1)
	554	+ w[0] = w[0]/(2.0*(g2 - g1))
	555	+ #
	556	+ # Imaginal
	557	+ #
	558	+ w[1] = 2.0(g1 + 2.0g2 - g1g22)(math.atan(g2) - math.atan(g1)) \
	559	+ + (-1.0 + 2.0g1g2 + g2*2)math.log((1 + g22)/(1 + g12))
	560	+ w[1] = -4.0*g1 + w[1]/(g2 - g1)
	561	+ w[1] = -3.0 + w[1]/(g2 - g1)
	562	+ w[1] = (3.0w[1])/(2.0(g2 - g1)**2)
	563	+ #
	564	+ return w
	565	+
	566	+
	567	+def libtetrabz_polcmplx_1222(g1, g2):
	568	+ """
	569	+ 5, g4 = g3 = g2
	570	+ :param g1:
	571	+ :param g2:
	572	+ :return:
	573	+ """
	574	+ w = numpy.empty(2, dtype=numpy.float_)
	575	+ #
	576	+ # Real
	577	+ #
	578	+ w[0] = 2.0(-1.0 + g12 + 2.0g1g2)(math.atan(g2) - math.atan(g1)) \
	579	+ + (-2.0g1 - g2 + g12g2) * math.log((1.0 + g22)/(1.0 + g12))
	580	+ w[0] = 2.0(1.0 - g1*2) + w[0]/(g2 - g1)
	581	+ w[0] = g1 - w[0]/(g2 - g1)
	582	+ w[0] = 1.0 - (3.0*w[0])/(g2 - g1)
	583	+ w[0] = w[0]/(2.0*(g2 - g1))
	584	+ #
	585	+ # Imaginal
	586	+ #
	587	+ w[1] = 2.0(-2.0g1 - g2 + g1*2g2)*(math.atan(g2) - math.atan(g1)) \
	588	+ + (1.0 - g1*2 - 2.0g1g2) math.log((1.0 + g22)/(1.0 + g12))
	589	+ w[1] = 4.0*g1 + w[1]/(g2 - g1)
	590	+ w[1] = 1.0 + w[1]/(g2 - g1)
	591	+ w[1] = (3.0w[1])/(2.0(g2 - g1)**2)
	592	+ #
	593	+ return w
	594	+
	595	+
	596	+def libtetrabz_polcmplx_1211(g1, g2):
	597	+ """
	598	+ 6, g4 = g3 = g1
	599	+ :param g1:
	600	+ :param g2:
	601	+ :return:
	602	+ """
	603	+ w = numpy.empty(2, dtype=numpy.float_)
	604	+ #
	605	+ # Real
	606	+ #
	607	+ w[0] = 2.0(3.0g2*2 - 1.0)(math.atan(g2) - math.atan(g1)) \
	608	+ + g2(g22 - 3.0)math.log((1.0 + g22)/(1.0 + g12))
	609	+ w[0] = 2.0(1.0 - g1*2) + w[0]/(g2 - g1)
	610	+ w[0] = -5.0*g1 + w[0]/(g2 - g1)
	611	+ w[0] = -11.0 + (3.0*w[0])/(g2 - g1)
	612	+ w[0] = w[0]/(6.0*(g2 - g1))
	613	+ #
	614	+ # Imaginal
	615	+ #
	616	+ w[1] = 2.0g2(-3.0 + g2*2)(math.atan(g2) - math.atan(g1)) \
	617	+ + (1.0 - 3.0g22)math.log((1.0 + g22)/(1.0 + g12))
	618	+ w[1] = 4.0*g2 + w[1]/(g2 - g1)
	619	+ w[1] = 1.0 + w[1]/(g2 - g1)
	620	+ w[1] = w[1]/(2.0(g2 - g1)*2)
	621	+ #
	622	+ return w

--- a/src/libtetrabz/libtetrabz_polstat_mod.py

+++ b/src/libtetrabz/libtetrabz_polstat.py

		@@ -25,7 +25,7 @@ import numpy
25	25	import libtetrabz_common
26	26
27	27
28		-def libtetrabz_polstat(bvec, eig1, eig2):
	28	+def polstat(bvec, eig1, eig2):
29	29	"""
30	30	Compute Static polarization function
31	31	:param bvec:

		@@ -134,7 +134,7 @@ def libtetrabz_polstat(bvec, eig1, eig2):
134	134	ei2 = ei1[0:4, ib]
135	135	ej2 = ej1[0:4, 0:nb]
136	136	w2 = libtetrabz_polstat2(ei2, ej2)
137		- w1[0:nb, 0:4] += w2[1:nb, 0:4]
	137	+ w1[0:nb, 0:4] += w2[0:nb, 0:4]
138	138	#
139	139	else:
140	140	continue

		@@ -228,7 +228,7 @@ def libtetrabz_polstat2(ei1, ej1):
228	228	#
229	229	de = ej1[0:4, ib] - ei1[0:4]
230	230	w2 = libtetrabz_polstat3(de)
231		- w1[ib, indx[0:4]] += w2
	231	+ w1[ib, 0:4] += w2
232	232	#
233	233	return w1
234	234

Show on old repository browser

libtetrabz Forkpythonlibtetrabz

python: Commit

Commit MetaInfo

Log Message

Cambiar Resumen

Diferencia incremental

libtetrabz
Fork
python
libtetrabz