OSDN > ソフトウェアを探す > 科学/工学 > 物理学 > libtetrabz > SCM > git > python > コミット

libtetrabz
Fork
python
libtetrabz

R/O
HTTP
SSH
HTTPS

python: コミット

libtetrabz python package

コミットメタ情報

リビジョン	5ce1d1d3343db3cc752bd38baae5f5488ad29662 (tree)
日時	2022-03-31 15:54:18
作者	Mitsuaki Kawamura <kawamitsuaki@gmai...>
コミッター	Mitsuaki Kawamura

ログメッセージ

Forgot to add C source

変更サマリ

add: src/libtetrabz/libtetrabz.c (diff)

差分

--- /dev/null

+++ b/src/libtetrabz/libtetrabz.c

		@@ -0,0 +1,4116 @@
	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	+#define PY_SSIZE_T_CLEAN
	24	+#include <Python.h>
	25	+#include <stdio.h>
	26	+#include <stdlib.h>
	27	+#include <math.h>
	28	+
	29	+void libtetrabz_initialize(
	30	+ int ng[3],
	31	+ double bvec[3][3],
	32	+ double wlsm[20][4],
	33	+ int **ikv
	34	+ )
	35	+ {
	36	+ /*
	37	+ define shortest diagonal line & define type of tetragonal
	38	+ */
	39	+ int ii, i0, i1, i2, i3, it, itype, ivvec0[4], divvec[4][4], ivvec[6][20][3], ikv0[3], nt;
	40	+ double bvec2[3][3], bvec3[4][3], bnorm[4];
	41	+ double wlsm1[4][4] = { {1440.0, 0.0, 30.0, 0.0},
	42	+ {0.0, 1440.0, 0.0, 30.0},
	43	+ {30.0, 0.0, 1440.0, 0.0},
	44	+ {0.0, 30.0, 0.0, 1440.0} },
	45	+ wlsm2[4][4] = {{-38.0, 7.0, 17.0, -28.0},
	46	+ {-28.0, -38.0, 7.0, 17.0},
	47	+ {17.0, -28.0, -38.0, 7.0},
	48	+ {7.0, 17.0, -28.0, -38.0}},
	49	+ wlsm3[4][4] = {{-56.0, 9.0, -46.0, 9.0},
	50	+ {9.0, -56.0, 9.0, -46.0},
	51	+ {-46.0, 9.0, -56.0, 9.0},
	52	+ {9.0, -46.0, 9.0, -56.0}},
	53	+ wlsm4[4][4] = {{-38.0, -28.0, 17.0, 7.0},
	54	+ {7.0, -38.0, -28.0, 17.0},
	55	+ {17.0, 7.0, -38.0, -28.0},
	56	+ {-28.0, 17.0, 7.0, -38.0}},
	57	+ wlsm5[4][4] = {{-18.0, -18.0, 12.0, -18.0},
	58	+ {-18.0, -18.0, -18.0, 12.0},
	59	+ {12.0, -18.0, -18.0, -18.0},
	60	+ {-18.0, 12.0, -18.0, -18.0}};
	61	+
	62	+ for (i1 = 0; i1 < 3; i1++)
	63	+ for (i2 = 0; i2 < 3; i2++)
	64	+ bvec2[i1][i2] = bvec[i1][i2] / (double) ng[i1];
	65	+
	66	+ for (i1 = 0; i1 < 3; i1++) {
	67	+ bvec3[0][i1] = -bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
	68	+ bvec3[1][i1] = bvec2[0][i1] - bvec2[1][i1] + bvec2[2][i1];
	69	+ bvec3[2][i1] = bvec2[0][i1] + bvec2[1][i1] - bvec2[2][i1];
	70	+ bvec3[3][i1] = bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
	71	+ }
	72	+ /*
	73	+ length of delta bvec
	74	+ */
	75	+ for (i1 = 0; i1 < 4; i1++) {
	76	+ bnorm[i1] = 0.0;
	77	+ for (i2 = 0; i2 < 3; i2++)
	78	+ bnorm[i1] += bvec3[i1][i2] * bvec3[i1][i2];
	79	+ }
	80	+ itype = 0;
	81	+ for (i1 = 1; i1 < 4; i1++)
	82	+ if (bnorm[i1] < bnorm[itype]) itype = i1;
	83	+ /*
	84	+ start & last
	85	+ */
	86	+ for (i0 = 0; i0 < 4; i0++) {
	87	+ ivvec0[i0] = 0;
	88	+ for (i1 = 0; i1 < 4; i1++)divvec[i0][i1] = 0;
	89	+ divvec[i0][i0] = 1;
	90	+ }
	91	+ ivvec0[itype] = 1;
	92	+ divvec[itype][itype] = -1;
	93	+ /*
	94	+ Corners of tetrahedron
	95	+ */
	96	+ it = 0;
	97	+ for (i0 = 0; i0 < 3; i0++) {
	98	+ for (i1 = 0; i1 < 3; i1++) {
	99	+ if (i1 == i0) continue;
	100	+ for (i2 = 0; i2 < 3; i2++) {
	101	+ if (i2 == i1 \|\| i2 == i0) continue;
	102	+
	103	+ for (i3 = 0; i3 < 3; i3++) {
	104	+ ivvec[it][0][i3] = ivvec0[i3];
	105	+ ivvec[it][1][i3] = ivvec[it][0][i3] + divvec[i0][i3];
	106	+ ivvec[it][2][i3] = ivvec[it][1][i3] + divvec[i1][i3];
	107	+ ivvec[it][3][i3] = ivvec[it][2][i3] + divvec[i2][i3];
	108	+ }
	109	+
	110	+ it += 1;
	111	+ }
	112	+ }
	113	+ }
	114	+ /*
	115	+ Additional points
	116	+ */
	117	+ for (i1 = 0; i1 < 6; i1++) {
	118	+ for (i2 = 0; i2 < 3; i2++) {
	119	+ ivvec[i1][4][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][1][i2];
	120	+ ivvec[i1][5][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][2][i2];
	121	+ ivvec[i1][6][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][3][i2];
	122	+ ivvec[i1][7][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][0][i2];
	123	+
	124	+ ivvec[i1][8][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][2][i2];
	125	+ ivvec[i1][9][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][3][i2];
	126	+ ivvec[i1][10][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][0][i2];
	127	+ ivvec[i1][11][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][1][i2];
	128	+
	129	+ ivvec[i1][12][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][3][i2];
	130	+ ivvec[i1][13][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][0][i2];
	131	+ ivvec[i1][14][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][1][i2];
	132	+ ivvec[i1][15][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][2][i2];
	133	+
	134	+ ivvec[i1][16][i2] = ivvec[i1][3][i2] - ivvec[i1][0][i2] + ivvec[i1][1][i2];
	135	+ ivvec[i1][17][i2] = ivvec[i1][0][i2] - ivvec[i1][1][i2] + ivvec[i1][2][i2];
	136	+ ivvec[i1][18][i2] = ivvec[i1][1][i2] - ivvec[i1][2][i2] + ivvec[i1][3][i2];
	137	+ ivvec[i1][19][i2] = ivvec[i1][2][i2] - ivvec[i1][3][i2] + ivvec[i1][0][i2];
	138	+ }
	139	+ }
	140	+
	141	+ for (i1 = 0; i1 < 4; i1++) {
	142	+ for (i2 = 0; i2 < 4; i2++) {
	143	+ wlsm[i2][i1] = wlsm1[i1][i2] /= 1260.0;
	144	+ wlsm[i2 + 4][i1] = wlsm2[i1][i2] /= 1260.0;
	145	+ wlsm[i2 + 8][i1] = wlsm3[i1][i2] /= 1260.0;
	146	+ wlsm[i2 + 12][i1] = wlsm4[i1][i2] /= 1260.0;
	147	+ wlsm[i2 + 16][i1] = wlsm5[i1][i2] /= 1260.0;
	148	+ }
	149	+ }
	150	+ /*
	151	+ k-index for energy
	152	+ */
	153	+ nt = 0;
	154	+ for (i2 = 0; i2 < ng[2]; i2++) {
	155	+ for (i1 = 0; i1 < ng[1]; i1++) {
	156	+ for (i0 = 0; i0 < ng[0]; i0++) {
	157	+
	158	+ for (it = 0; it < 6; it++) {
	159	+
	160	+ for (ii = 0; ii < 20; ii++) {
	161	+ ikv0[0] = (i0 + ivvec[it][ii][0]) % ng[0];
	162	+ ikv0[1] = (i1 + ivvec[it][ii][1]) % ng[1];
	163	+ ikv0[2] = (i2 + ivvec[it][ii][2]) % ng[2];
	164	+ for (i3 = 0; i3 < 3; i3++) if (ikv0[i3] < 0) ikv0[i3] += ng[i3];
	165	+ ikv[nt][ii] = ikv0[2] + ng[2] * ikv0[1] + ng[2] * ng[1] * ikv0[0];
	166	+ }
	167	+ nt += 1;
	168	+ }
	169	+ }
	170	+ }
	171	+ }
	172	+}
	173	+/*
	174	+Cut small tetrahedron A1
	175	+*/
	176	+void libtetrabz_tsmall_a1(
	177	+ double *e,
	178	+ double e0,
	179	+ double *v,
	180	+ double tsmall[4][4]
	181	+) {
	182	+ double a10, a20, a30;
	183	+ a10 = (e0 - e[0]) / (e[1] - e[0]);
	184	+ a20 = (e0 - e[0]) / (e[2] - e[0]);
	185	+ a30 = (e0 - e[0]) / (e[3] - e[0]);
	186	+
	187	+ v = a10 a20 * a30;
	188	+
	189	+ tsmall[0][0] = 1.0;
	190	+ tsmall[0][1] = 1.0 - a10;
	191	+ tsmall[0][2] = 1.0 - a20;
	192	+ tsmall[0][3] = 1.0 - a30;
	193	+ tsmall[1][0] = 0.0;
	194	+ tsmall[1][1] = a10;
	195	+ tsmall[1][2] = 0.0;
	196	+ tsmall[1][3] = 0.0;
	197	+ tsmall[2][0] = 0.0;
	198	+ tsmall[2][1] = 0.0;
	199	+ tsmall[2][2] = a20;
	200	+ tsmall[2][3] = 0.0;
	201	+ tsmall[3][0] = 0.0;
	202	+ tsmall[3][1] = 0.0;
	203	+ tsmall[3][2] = 0.0;
	204	+ tsmall[3][3] = a30;
	205	+}
	206	+/*
	207	+Cut small tetrahedron B1
	208	+*/
	209	+void libtetrabz_tsmall_b1(
	210	+ double *e,
	211	+ double e0,
	212	+ double *v,
	213	+ double tsmall[4][4]
	214	+)
	215	+{
	216	+ double a13, a20, a30;
	217	+ a13 = (e0 - e[3]) / (e[1] - e[3]);
	218	+ a20 = (e0 - e[0]) / (e[2] - e[0]);
	219	+ a30 = (e0 - e[0]) / (e[3] - e[0]);
	220	+
	221	+ v = a20 a30 * a13;
	222	+
	223	+ tsmall[0][0] = 1.0;
	224	+ tsmall[0][1] = 1.0 - a20;
	225	+ tsmall[0][2] = 1.0 - a30;
	226	+ tsmall[0][3] = 0.0;
	227	+ tsmall[1][0] = 0.0;
	228	+ tsmall[1][1] = 0.0;
	229	+ tsmall[1][2] = 0.0;
	230	+ tsmall[1][3] = a13;
	231	+ tsmall[2][0] = 0.0;
	232	+ tsmall[2][1] = a20;
	233	+ tsmall[2][2] = 0.0;
	234	+ tsmall[2][3] = 0.0;
	235	+ tsmall[3][0] = 0.0;
	236	+ tsmall[3][1] = 0.0;
	237	+ tsmall[3][2] = a30;
	238	+ tsmall[3][3] = 1.0 - a13;
	239	+}
	240	+/*
	241	+Cut small tetrahedron B2
	242	+*/
	243	+void libtetrabz_tsmall_b2(
	244	+ double *e,
	245	+ double e0,
	246	+ double *v,
	247	+ double tsmall[4][4]
	248	+)
	249	+{
	250	+ double a21, a31;
	251	+ a21 = (e0 - e[1]) / (e[2] - e[1]);
	252	+ a31 = (e0 - e[1]) / (e[3] - e[1]);
	253	+
	254	+ v = a21 a31;
	255	+
	256	+ tsmall[0][0] = 1.0;
	257	+ tsmall[0][1] = 0.0;
	258	+ tsmall[0][2] = 0.0;
	259	+ tsmall[0][3] = 0.0;
	260	+ tsmall[1][0] = 0.0;
	261	+ tsmall[1][1] = 1.0;
	262	+ tsmall[1][2] = 1.0 - a21;
	263	+ tsmall[1][3] = 1.0 - a31;
	264	+ tsmall[2][0] = 0.0;
	265	+ tsmall[2][1] = 0.0;
	266	+ tsmall[2][2] = a21;
	267	+ tsmall[2][3] = 0.0;
	268	+ tsmall[3][0] = 0.0;
	269	+ tsmall[3][1] = 0.0;
	270	+ tsmall[3][2] = 0.0;
	271	+ tsmall[3][3] = a31;
	272	+}
	273	+/*
	274	+Cut small tetrahedron B3
	275	+*/
	276	+void libtetrabz_tsmall_b3(
	277	+ double *e,
	278	+ double e0,
	279	+ double *v,
	280	+ double tsmall[4][4]
	281	+)
	282	+{
	283	+ double a12, a20, a31;
	284	+ a12 = (e0 - e[2]) / (e[1] - e[2]);
	285	+ a20 = (e0 - e[0]) / (e[2] - e[0]);
	286	+ a31 = (e0 - e[1]) / (e[3] - e[1]);
	287	+
	288	+ v = a12 a20 * a31;
	289	+
	290	+ tsmall[0][0] = 1.0;
	291	+ tsmall[0][1] = 1.0 - a20;
	292	+ tsmall[0][2] = 0.0;
	293	+ tsmall[0][3] = 0.0;
	294	+ tsmall[1][0] = 0.0;
	295	+ tsmall[1][1] = 0.0;
	296	+ tsmall[1][2] = a12;
	297	+ tsmall[1][3] = 1.0 - a31;
	298	+ tsmall[2][0] = 0.0;
	299	+ tsmall[2][1] = a20;
	300	+ tsmall[2][2] = 1.0 - a12;
	301	+ tsmall[2][3] = 0.0;
	302	+ tsmall[3][0] = 0.0;
	303	+ tsmall[3][1] = 0.0;
	304	+ tsmall[3][2] = 0.0;
	305	+ tsmall[3][3] = a31;
	306	+}
	307	+/*
	308	+Cut small tetrahedron C1
	309	+*/
	310	+void libtetrabz_tsmall_c1(
	311	+ double *e,
	312	+ double e0,
	313	+ double *v,
	314	+ double tsmall[4][4]
	315	+)
	316	+{
	317	+ double a32;
	318	+ a32 = (e0 - e[2]) / (e[3] - e[2]);
	319	+
	320	+ *v = a32;
	321	+
	322	+ tsmall[0][0] = 1.0;
	323	+ tsmall[0][1] = 0.0;
	324	+ tsmall[0][2] = 0.0;
	325	+ tsmall[0][3] = 0.0;
	326	+ tsmall[1][0] = 0.0;
	327	+ tsmall[1][1] = 1.0;
	328	+ tsmall[1][2] = 0.0;
	329	+ tsmall[1][3] = 0.0;
	330	+ tsmall[2][0] = 0.0;
	331	+ tsmall[2][1] = 0.0;
	332	+ tsmall[2][2] = 1.0;
	333	+ tsmall[2][3] = 1.0 - a32;
	334	+ tsmall[3][0] = 0.0;
	335	+ tsmall[3][1] = 0.0;
	336	+ tsmall[3][2] = 0.0;
	337	+ tsmall[3][3] = a32;
	338	+}
	339	+/*
	340	+Cut small tetrahedron C2
	341	+*/
	342	+void libtetrabz_tsmall_c2(
	343	+ double *e,
	344	+ double e0,
	345	+ double *v,
	346	+ double tsmall[4][4]
	347	+)
	348	+{
	349	+ double a23, a31;
	350	+ a23 = (e0 - e[3]) / (e[2] - e[3]);
	351	+ a31 = (e0 - e[1]) / (e[3] - e[1]);
	352	+
	353	+ v = a23 a31;
	354	+
	355	+ tsmall[0][0] = 1.0;
	356	+ tsmall[0][1] = 0.0;
	357	+ tsmall[0][2] = 0.0;
	358	+ tsmall[0][3] = 0.0;
	359	+ tsmall[1][0] = 0.0;
	360	+ tsmall[1][1] = 1.0;
	361	+ tsmall[1][2] = 1.0 - a31;
	362	+ tsmall[1][3] = 0.0;
	363	+ tsmall[2][0] = 0.0;
	364	+ tsmall[2][1] = 0.0;
	365	+ tsmall[2][2] = 0.0;
	366	+ tsmall[2][3] = a23;
	367	+ tsmall[3][0] = 0.0;
	368	+ tsmall[3][1] = 0.0;
	369	+ tsmall[3][2] = a31;
	370	+ tsmall[3][3] = 1.0 - a23;
	371	+}
	372	+/*
	373	+Cut small tetrahedron C3
	374	+*/
	375	+void libtetrabz_tsmall_c3(
	376	+ double *e,
	377	+ double e0,
	378	+ double *v,
	379	+ double tsmall[4][4]
	380	+)
	381	+{
	382	+ double a23, a13, a30;
	383	+ a23 = (e0 - e[3]) / (e[2] - e[3]);
	384	+ a13 = (e0 - e[3]) / (e[1] - e[3]);
	385	+ a30 = (e0 - e[0]) / (e[3] - e[0]);
	386	+
	387	+ v = a23 a13 * a30;
	388	+
	389	+ tsmall[0][0] = 1.0;
	390	+ tsmall[0][1] = 1.0 - a30;
	391	+ tsmall[0][2] = 0.0;
	392	+ tsmall[0][3] = 0.0;
	393	+ tsmall[1][0] = 0.0;
	394	+ tsmall[1][1] = 0.0;
	395	+ tsmall[1][2] = a13;
	396	+ tsmall[1][3] = 0.0;
	397	+ tsmall[2][0] = 0.0;
	398	+ tsmall[2][1] = 0.0;
	399	+ tsmall[2][2] = 0.0;
	400	+ tsmall[2][3] = a23;
	401	+ tsmall[3][0] = 0.0;
	402	+ tsmall[3][1] = a30;
	403	+ tsmall[3][2] = 1.0 - a13;
	404	+ tsmall[3][3] = 1.0 - a23;
	405	+}
	406	+/*
	407	+Cut triangle A1
	408	+*/
	409	+void libtetrabz_triangle_a1(
	410	+ double *e,
	411	+ double e0,
	412	+ double *v,
	413	+ double tsmall[4][3]
	414	+)
	415	+{
	416	+ double a10, a20, a30;
	417	+ a10 = (e0 - e[0]) / (e[1] - e[0]);
	418	+ a20 = (e0 - e[0]) / (e[2] - e[0]);
	419	+ a30 = (e0 - e[0]) / (e[3] - e[0]);
	420	+
	421	+ v = 3.0 a10 * a20 / (e[3] - e[0]);
	422	+
	423	+ tsmall[0][0] = 1.0 - a10;
	424	+ tsmall[0][1] = 1.0 - a20;
	425	+ tsmall[0][2] = 1.0 - a30;
	426	+ tsmall[1][0] = a10;
	427	+ tsmall[1][1] = 0.0;
	428	+ tsmall[1][2] = 0.0;
	429	+ tsmall[2][0] = 0.0;
	430	+ tsmall[2][1] = a20;
	431	+ tsmall[2][2] = 0.0;
	432	+ tsmall[3][0] = 0.0;
	433	+ tsmall[3][1] = 0.0;
	434	+ tsmall[3][2] = a30;
	435	+}
	436	+/*
	437	+Cut triangle B1
	438	+*/
	439	+void libtetrabz_triangle_b1(
	440	+ double *e,
	441	+ double e0,
	442	+ double *v,
	443	+ double tsmall[4][3]
	444	+) {
	445	+ double a30, a13, a20;
	446	+ a30 = (e0 - e[0]) / (e[3] - e[0]);
	447	+ a13 = (e0 - e[3]) / (e[1] - e[3]);
	448	+ a20 = (e0 - e[0]) / (e[2] - e[0]);
	449	+
	450	+ v = 3.0 a30 * a13 / (e[2] - e[0]);
	451	+
	452	+ tsmall[0][0] = 1.0 - a20;
	453	+ tsmall[0][1] = 1.0 - a30;
	454	+ tsmall[0][2] = 0.0;
	455	+ tsmall[1][0] = 0.0;
	456	+ tsmall[1][1] = 0.0;
	457	+ tsmall[1][2] = a13;
	458	+ tsmall[2][0] = a20;
	459	+ tsmall[2][1] = 0.0;
	460	+ tsmall[2][2] = 0.0;
	461	+ tsmall[3][0] = 0.0;
	462	+ tsmall[3][1] = a30;
	463	+ tsmall[3][2] = 1.0 - a13;
	464	+}
	465	+/*
	466	+Cut triangle B2
	467	+*/
	468	+void libtetrabz_triangle_b2(
	469	+ double *e,
	470	+ double e0,
	471	+ double *v,
	472	+ double tsmall[4][3]
	473	+)
	474	+{
	475	+ double a12, a31, a20;
	476	+ a12 = (e0 - e[2]) / (e[1] - e[2]);
	477	+ a31 = (e0 - e[1]) / (e[3] - e[1]);
	478	+ a20 = (e0 - e[0]) / (e[2] - e[0]);
	479	+
	480	+ v = 3.0 a12 * a31 / (e[2] - e[0]);
	481	+
	482	+ tsmall[0][0] = 1.0 - a20;
	483	+ tsmall[0][1] = 0.0;
	484	+ tsmall[0][2] = 0.0;
	485	+ tsmall[1][0] = 0.0;
	486	+ tsmall[1][1] = a12;
	487	+ tsmall[1][2] = 1.0 - a31;
	488	+ tsmall[2][0] = a20;
	489	+ tsmall[2][1] = 1.0 - a12;
	490	+ tsmall[2][2] = 0.0;
	491	+ tsmall[3][0] = 0.0;
	492	+ tsmall[3][1] = 0.0;
	493	+ tsmall[3][2]= a31;
	494	+}
	495	+/*
	496	+Cut triangle C1
	497	+*/
	498	+void libtetrabz_triangle_c1(
	499	+ double *e,
	500	+ double e0,
	501	+ double *v,
	502	+ double tsmall[4][3]
	503	+)
	504	+{
	505	+ double a03, a13, a23;
	506	+ a03 = (e0 - e[3]) / (e[0] - e[3]);
	507	+ a13 = (e0 - e[3]) / (e[1] - e[3]);
	508	+ a23 = (e0 - e[3]) / (e[2] - e[3]);
	509	+
	510	+ v = 3.0 a03 * a13 / (e[3] - e[2]);
	511	+
	512	+ tsmall[0][0] = a03;
	513	+ tsmall[0][1] = 0.0;
	514	+ tsmall[0][2] = 0.0;
	515	+ tsmall[1][0] = 0.0;
	516	+ tsmall[1][1] = a13;
	517	+ tsmall[1][2] = 0.0;
	518	+ tsmall[2][0] = 0.0;
	519	+ tsmall[2][1] = 0.0;
	520	+ tsmall[2][2] = a23;
	521	+ tsmall[3][0] = 1.0 - a03;
	522	+ tsmall[3][1] = 1.0 - a13;
	523	+ tsmall[3][2] = 1.0 - a23;
	524	+}
	525	+/*
	526	+Sort eigenvalues
	527	+*/
	528	+void eig_sort(
	529	+ int n, //!< [in] the number of components
	530	+ double *key, //!< [in] Variables to be sorted [n].
	531	+ int *swap //!< [out] Order of index (sorted)
	532	+)
	533	+{
	534	+ int i, j, k, min_loc;
	535	+ double min_val;
	536	+
	537	+ for (i = 0; i < n; ++i) swap[i] = i;
	538	+
	539	+ for (i = 0; i < n - 1; ++i) {
	540	+ min_val = key[i];
	541	+ min_loc = i;
	542	+ for (j = i + 1; j < n; ++j) {
	543	+ if (min_val > key[j]) {
	544	+ min_val = key[j];
	545	+ min_loc = j;
	546	+ }
	547	+ }
	548	+ if (key[i] > min_val) {
	549	+ /*
	550	+ Swap
	551	+ */
	552	+ key[min_loc] = key[i];
	553	+ key[i] = min_val;
	554	+
	555	+ k = swap[min_loc];
	556	+ swap[min_loc] = swap[i];
	557	+ swap[i] = k;
	558	+ }
	559	+ }/for (i = 0; i < n - 1; ++i)/
	560	+}/eig_sort/
	561	+/*
	562	+2nd step of tetrahedron method.
	563	+*/
	564	+static void libtetrabz_dbldelta2(
	565	+ int nb,
	566	+ double **ej,
	567	+ double **w
	568	+) {
	569	+ int i3, ib, indx[3];
	570	+ double a10, a20, a02, a12, v, e[3], e_abs;
	571	+
	572	+ for (ib = 0; ib < nb; ib++) {
	573	+
	574	+ e_abs = 0.0;
	575	+ for (i3 = 0; i3 < 3; i3++) {
	576	+ e[i3] = ej[ib][i3];
	577	+ if (e_abs < fabs(e[i3])) e_abs = fabs(e[i3]);
	578	+ }
	579	+ eig_sort(3, e, indx);
	580	+
	581	+ if (e_abs < 1.0e-10) {
	582	+ printf("Nesting ##\n");
	583	+ }
	584	+
	585	+ if ((e[0] < 0.0 && 0.0 <= e[1]) \|\| (e[0] <= 0.0 && 0.0 < e[1])) {
	586	+
	587	+ a10 = (0.0 - e[0]) / (e[1] - e[0]);
	588	+ a20 = (0.0 - e[0]) / (e[2] - e[0]);
	589	+
	590	+ v = a10 / (e[2] - e[0]);
	591	+
	592	+ w[ib][indx[0]] = v * (2.0 - a10 - a20);
	593	+ w[ib][indx[1]] = v * a10;
	594	+ w[ib][indx[2]] = v * a20;
	595	+ }
	596	+ else if ((e[1] <= 0.0 && 0.0 < e[2]) \|\| (e[1] < 0.0 && 0.0 <= e[2])) {
	597	+
	598	+ a02 = (0.0 - e[2]) / (e[0] - e[2]);
	599	+ a12 = (0.0 - e[2]) / (e[1] - e[2]);
	600	+
	601	+ v = a12 / (e[2] - e[0]);
	602	+
	603	+ w[ib][indx[0]] = v * a02;
	604	+ w[ib][indx[1]] = v * a12;
	605	+ w[ib][indx[2]] = v * (2.0 - a02 - a12);
	606	+ }
	607	+ else {
	608	+ for (i3 = 0; i3 < 3; i3++)
	609	+ w[ib][i3] = 0.0;
	610	+ }
	611	+ }
	612	+}
	613	+/*
	614	+Main SUBROUTINE for Delta(E1) * Delta(E2)
	615	+*/
	616	+static PyObject* dbldelta_c(PyObject* self, PyObject* args)
	617	+{
	618	+ int it, ik, ib, i20, i4, j3, jb, ** ikv, indx[4], ierr, ng[3], nk, nb;
	619	+ double wlsm[20][4], ei1, ej1, ej2, e[4], * w1, **w2, v, tsmall[4][3], thr,
	620	+ bvec[3][3], eig1, eig2, *** wght;
	621	+ PyObject* eig1_po, * eig2_po, * wght_po;
	622	+ /*
	623	+ Read input from python object
	624	+ */
	625	+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOO",
	626	+ &ng[0], &ng[1], &ng[2], &nk, &nb,
	627	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	628	+ &eig1_po, &eig2_po))
	629	+ return NULL;
	630	+ /*
	631	+ convert python list object to array
	632	+ */
	633	+ eig1 = (double*)malloc(nk sizeof(double*));
	634	+ eig1[0] = (double)malloc(nk nb * sizeof(double));
	635	+ eig2 = (double*)malloc(nk sizeof(double*));
	636	+ eig2[0] = (double)malloc(nk nb * sizeof(double));
	637	+ wght = (double**)malloc(nk sizeof(double**));
	638	+ wght[0] = (double*)malloc(nk nb * sizeof(double*));
	639	+ wght[0][0] = (double)malloc(nk nb * nb * sizeof(double));
	640	+ for (ik = 0; ik < nk; ik++) {
	641	+ eig1[ik] = eig1[0] + ik * nb;
	642	+ eig2[ik] = eig2[0] + ik * nb;
	643	+ wght[ik] = wght[0] + ik * nb;
	644	+ for (ib = 0; ib < nb; ib++) {
	645	+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
	646	+ }
	647	+ }
	648	+
	649	+ for (ik = 0; ik < nk; ik++)
	650	+ for (ib = 0; ib < nb; ib++) {
	651	+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
	652	+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
	653	+ }
	654	+ /*
	655	+ Start main calculation
	656	+ */
	657	+ thr = 1.0e-10;
	658	+
	659	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	660	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	661	+ for (ik = 0; ik < 6 * nk; ik++) {
	662	+ ikv[ik] = ikv[0] + ik * 20;
	663	+ }
	664	+
	665	+ ei1 = (double*)malloc(4 sizeof(double*));
	666	+ ej1 = (double*)malloc(4 sizeof(double*));
	667	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	668	+ ej1[0] = (double)malloc(4 nb * sizeof(double));
	669	+ for (i4 = 0; i4 < 4; i4++) {
	670	+ ei1[i4] = ei1[0] + i4 * nb;
	671	+ ej1[i4] = ej1[0] + i4 * nb;
	672	+ }
	673	+
	674	+ w1 = (double**)malloc(nb sizeof(double**));
	675	+ w1[0] = (double*)malloc(nb 4 * sizeof(double*));
	676	+ w1[0][0] = (double)malloc(nb 4 * nb * sizeof(double));
	677	+ for (ib = 0; ib < nb; ib++) {
	678	+ w1[ib] = w1[0] + ib * 4;
	679	+ for (i4 = 0; i4 < 4; i4++) {
	680	+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
	681	+ }
	682	+ }
	683	+
	684	+ ej2 = (double*)malloc(nb sizeof(double*));
	685	+ ej2[0] = (double)malloc(nb 3 * sizeof(double));
	686	+ for (ib = 0; ib < nb; ib++) {
	687	+ ej2[ib] = ej2[0] + ib * 3;
	688	+ }
	689	+
	690	+ w2 = (double*)malloc(nb sizeof(double*));
	691	+ w2[0] = (double)malloc(nb 3 * sizeof(double));
	692	+ for (ib = 0; ib < nb; ib++) {
	693	+ w2[ib] = w2[0] + ib * 3;
	694	+ }
	695	+
	696	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	697	+
	698	+ for (ik = 0; ik < nk; ik++)
	699	+ for (ib = 0; ib < nb; ib++)
	700	+ for (jb = 0; jb < nb; jb++)
	701	+ wght[ik][ib][jb] = 0.0;
	702	+
	703	+ for (it = 0; it < 6 * nk; it++) {
	704	+
	705	+ for (i4 = 0; i4 < 4; i4++)
	706	+ for (ib = 0; ib < nb; ib++) {
	707	+ ei1[i4][ib] = 0.0;
	708	+ ej1[i4][ib] = 0.0;
	709	+ }
	710	+ for (i20 = 0; i20 < 20; i20++) {
	711	+ for (i4 = 0; i4 < 4; i4++) {
	712	+ for (ib = 0; ib < nb; ib++) {
	713	+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
	714	+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
	715	+ }
	716	+ }
	717	+ }
	718	+
	719	+ for (ib = 0; ib < nb; ib++) {
	720	+
	721	+ for (i4 = 0; i4 < 4; i4++)
	722	+ for (jb = 0; jb < nb; jb++)
	723	+ w1[ib][i4][jb] = 0.0;
	724	+
	725	+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
	726	+ eig_sort(4, e, indx);
	727	+
	728	+ if (e[0] < 0.0 && 0.0 <= e[1]) {
	729	+
	730	+ libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
	731	+
	732	+ if (v > thr) {
	733	+ for (jb = 0; jb < nb; jb++)
	734	+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
	735	+ for (i4 = 0; i4 < 4; i4++)
	736	+ for (jb = 0; jb < nb; jb++)
	737	+ for (j3 = 0; j3 < 3; j3++)
	738	+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
	739	+ libtetrabz_dbldelta2(nb, ej2, w2);
	740	+ for (i4 = 0; i4 < 4; i4++)
	741	+ for (jb = 0; jb < nb; jb++)
	742	+ for (j3 = 0; j3 < 3; j3++)
	743	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
	744	+ }
	745	+ }
	746	+ else if (e[1] < 0.0 && 0.0 <= e[2]) {
	747	+
	748	+ libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
	749	+
	750	+ if (v > thr) {
	751	+ for (jb = 0; jb < nb; jb++)
	752	+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
	753	+ for (i4 = 0; i4 < 4; i4++)
	754	+ for (jb = 0; jb < nb; jb++)
	755	+ for (j3 = 0; j3 < 3; j3++)
	756	+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
	757	+ libtetrabz_dbldelta2(nb, ej2, w2);
	758	+ for (i4 = 0; i4 < 4; i4++)
	759	+ for (jb = 0; jb < nb; jb++)
	760	+ for (j3 = 0; j3 < 3; j3++)
	761	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
	762	+ }
	763	+
	764	+ libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
	765	+
	766	+ if (v > thr) {
	767	+ for (jb = 0; jb < nb; jb++)
	768	+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
	769	+ for (i4 = 0; i4 < 4; i4++)
	770	+ for (jb = 0; jb < nb; jb++)
	771	+ for (j3 = 0; j3 < 3; j3++)
	772	+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
	773	+ libtetrabz_dbldelta2(nb, ej2, w2);
	774	+ for (i4 = 0; i4 < 4; i4++)
	775	+ for (jb = 0; jb < nb; jb++)
	776	+ for (j3 = 0; j3 < 3; j3++)
	777	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
	778	+ }
	779	+ }
	780	+ else if (e[2] < 0.0 && 0.0 < e[3]) {
	781	+
	782	+ libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
	783	+
	784	+ if (v > thr) {
	785	+ for (jb = 0; jb < nb; jb++)
	786	+ for (j3 = 0; j3 < 3; j3++) ej2[jb][j3] = 0.0;
	787	+ for (i4 = 0; i4 < 4; i4++)
	788	+ for (jb = 0; jb < nb; jb++)
	789	+ for (j3 = 0; j3 < 3; j3++)
	790	+ ej2[jb][j3] += ej1[indx[i4]][jb] * tsmall[i4][j3];
	791	+ libtetrabz_dbldelta2(nb, ej2, w2);
	792	+ for (i4 = 0; i4 < 4; i4++)
	793	+ for (jb = 0; jb < nb; jb++)
	794	+ for (j3 = 0; j3 < 3; j3++)
	795	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j3] * w2[jb][j3];
	796	+ }
	797	+ }
	798	+ else {
	799	+ continue;
	800	+ }
	801	+ }
	802	+ for (i20 = 0; i20 < 20; i20++)
	803	+ for (ib = 0; ib < nb; ib++)
	804	+ for (i4 = 0; i4 < 4; i4++)
	805	+ for (jb = 0; jb < nb; jb++)
	806	+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
	807	+ }
	808	+ for (ik = 0; ik < nk; ik++)
	809	+ for (ib = 0; ib < nb; ib++)
	810	+ for (jb = 0; jb < nb; jb++)
	811	+ wght[ik][ib][jb] /= (6.0 * (double)nk);
	812	+
	813	+ free(ikv[0]);
	814	+ free(ikv);
	815	+ free(ei1[0]);
	816	+ free(ei1);
	817	+ free(ej1[0]);
	818	+ free(ej1);
	819	+ free(ej2[0]);
	820	+ free(ej2);
	821	+ free(w1[0][0]);
	822	+ free(w1[0]);
	823	+ free(w1);
	824	+ free(w2[0]);
	825	+ free(w2);
	826	+ /*
	827	+ Convert weight to python list object
	828	+ */
	829	+ wght_po = PyList_New(nk * nb * nb);
	830	+ ierr = 0;
	831	+ for (ik = 0; ik < nk; ik++)
	832	+ for (ib = 0; ib < nb; ib++)
	833	+ for (jb = 0; jb < nb; jb++)
	834	+ ierr = PyList_SetItem(wght_po, ik * nb * nb + ib * nb + jb, PyFloat_FromDouble(wght[ik][ib][jb]));
	835	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	836	+
	837	+ free(eig1[0]);
	838	+ free(eig1);
	839	+ free(eig2[0]);
	840	+ free(eig2);
	841	+ free(wght[0][0]);
	842	+ free(wght[0]);
	843	+ free(wght);
	844	+
	845	+ return wght_po;
	846	+}
	847	+/*
	848	+Tetrahedron method for theta( - de)
	849	+*/
	850	+static void libtetrabz_dblstep2(
	851	+ int nb,
	852	+ double ei[4],
	853	+ double **ej,
	854	+ double **w1
	855	+) {
	856	+ int i4, j4, ib, indx[4];
	857	+ double v, thr, e[4], tsmall[4][4];
	858	+
	859	+ thr = 1.0e-8;
	860	+
	861	+ for (ib = 0; ib < nb; ib++) {
	862	+
	863	+ for (i4 = 0; i4 < 4; i4++)
	864	+ w1[ib][i4] = 0.0;
	865	+
	866	+ for (i4 = 0; i4 < 4; i4++) e[i4] = - ei[i4] + ej[ib][i4];
	867	+ eig_sort(4, e, indx);
	868	+
	869	+ if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
	870	+ /*
	871	+ Theta(0) = 0.5
	872	+ */
	873	+ v = 0.5;
	874	+ for (i4 = 0; i4 < 4; i4++)
	875	+ w1[ib][i4] += v * 0.25;
	876	+ }
	877	+ else if ((e[0] <= 0.0 && 0.0 < e[1]) \|\| (e[0] < 0.0 && 0.0 <= e[1])) {
	878	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	879	+ for (i4 = 0; i4 < 4; i4++)
	880	+ for (j4 = 0; j4 < 4; j4++)
	881	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
	882	+ }
	883	+ else if ((e[1] <= 0.0 && 0.0 < e[2]) \|\| (e[1] < 0.0 && 0.0 <= e[2])) {
	884	+
	885	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	886	+ for (i4 = 0; i4 < 4; i4++)
	887	+ for (j4 = 0; j4 < 4; j4++)
	888	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
	889	+
	890	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	891	+ for (i4 = 0; i4 < 4; i4++)
	892	+ for (j4 = 0; j4 < 4; j4++)
	893	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
	894	+
	895	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	896	+ for (i4 = 0; i4 < 4; i4++)
	897	+ for (j4 = 0; j4 < 4; j4++)
	898	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
	899	+ }
	900	+ else if ((e[2] <= 0.0 && 0.0 < e[3]) \|\| (e[2] < 0.0 && 0.0 <= e[3])) {
	901	+
	902	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	903	+ for (i4 = 0; i4 < 4; i4++)
	904	+ for (j4 = 0; j4 < 4; j4++)
	905	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
	906	+
	907	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	908	+ for (i4 = 0; i4 < 4; i4++)
	909	+ for (j4 = 0; j4 < 4; j4++)
	910	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
	911	+
	912	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	913	+ for (i4 = 0; i4 < 4; i4++)
	914	+ for (j4 = 0; j4 < 4; j4++)
	915	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * 0.25;
	916	+ }
	917	+ else if (e[3] <= 0.0) {
	918	+ for (i4 = 0; i4 < 4; i4++)
	919	+ w1[ib][i4] += 0.25;
	920	+ }
	921	+ }
	922	+}
	923	+/*
	924	+Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
	925	+*/
	926	+static PyObject* dblstep_c(PyObject* self, PyObject* args)
	927	+{
	928	+ int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4], ierr, ng[3], nk, nb;
	929	+ double wlsm[20][4], ei1, ej1, ei2[4], ej2, e[4], * w1, ** w2, v, tsmall[4][4], thr,
	930	+ bvec[3][3], eig1, eig2, *** wght;
	931	+ PyObject* eig1_po, * eig2_po, * wght_po;
	932	+
	933	+ /*
	934	+ Read input from python object
	935	+ */
	936	+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOO",
	937	+ &ng[0], &ng[1], &ng[2], &nk, &nb,
	938	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	939	+ &eig1_po, &eig2_po))
	940	+ return NULL;
	941	+ /*
	942	+ convert python list object to array
	943	+ */
	944	+ eig1 = (double*)malloc(nk sizeof(double*));
	945	+ eig1[0] = (double)malloc(nk nb * sizeof(double));
	946	+ eig2 = (double*)malloc(nk sizeof(double*));
	947	+ eig2[0] = (double)malloc(nk nb * sizeof(double));
	948	+ wght = (double**)malloc(nk sizeof(double**));
	949	+ wght[0] = (double*)malloc(nk nb * sizeof(double*));
	950	+ wght[0][0] = (double)malloc(nk nb * nb * sizeof(double));
	951	+ for (ik = 0; ik < nk; ik++) {
	952	+ eig1[ik] = eig1[0] + ik * nb;
	953	+ eig2[ik] = eig2[0] + ik * nb;
	954	+ wght[ik] = wght[0] + ik * nb;
	955	+ for (ib = 0; ib < nb; ib++) {
	956	+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
	957	+ }
	958	+ }
	959	+
	960	+ for (ik = 0; ik < nk; ik++)
	961	+ for (ib = 0; ib < nb; ib++) {
	962	+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
	963	+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
	964	+ }
	965	+ /*
	966	+ Start main calculation
	967	+ */
	968	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	969	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	970	+ for (ik = 0; ik < 6 * nk; ik++) {
	971	+ ikv[ik] = ikv[0] + ik * 20;
	972	+ }
	973	+
	974	+ ei1 = (double*)malloc(4 sizeof(double*));
	975	+ ej1 = (double*)malloc(4 sizeof(double*));
	976	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	977	+ ej1[0] = (double)malloc(4 nb * sizeof(double));
	978	+ for (i4 = 0; i4 < 4; i4++) {
	979	+ ei1[i4] = ei1[0] + i4 * nb;
	980	+ ej1[i4] = ej1[0] + i4 * nb;
	981	+ }
	982	+
	983	+ ej2 = (double*)malloc(nb sizeof(double*));
	984	+ ej2[0] = (double)malloc(nb 4 * sizeof(double));
	985	+ for (ib = 0; ib < nb; ib++) {
	986	+ ej2[ib] = ej2[0] + ib * 4;
	987	+ }
	988	+
	989	+ w1 = (double**)malloc(nb sizeof(double**));
	990	+ w1[0] = (double*)malloc(nb 4 * sizeof(double*));
	991	+ w1[0][0] = (double)malloc(nb 4 * nb * sizeof(double));
	992	+ for (ib = 0; ib < nb; ib++) {
	993	+ w1[ib] = w1[0] + ib * 4;
	994	+ for (i4 = 0; i4 < 4; i4++) {
	995	+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
	996	+ }
	997	+ }
	998	+
	999	+ w2 = (double*)malloc(nb sizeof(double*));
	1000	+ w2[0] = (double)malloc(nb 4 * sizeof(double));
	1001	+ for (ib = 0; ib < nb; ib++) {
	1002	+ w2[ib] = w2[0] + ib * 4;
	1003	+ }
	1004	+
	1005	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	1006	+
	1007	+ for (ik = 0; ik < nk; ik++)
	1008	+ for (ib = 0; ib < nb; ib++)
	1009	+ for (jb = 0; jb < nb; jb++)
	1010	+ wght[ik][ib][jb] = 0.0;
	1011	+
	1012	+ thr = 1.0e-8;
	1013	+
	1014	+ for(it = 0; it < 6*nk; it++) {
	1015	+
	1016	+ for (i4 = 0; i4 < 4; i4++)
	1017	+ for (ib = 0; ib < nb; ib++) {
	1018	+ ei1[i4][ib] = 0.0;
	1019	+ ej1[i4][ib] = 0.0;
	1020	+ }
	1021	+ for (i20 = 0; i20 < 20; i20++) {
	1022	+ for (i4 = 0; i4 < 4; i4++) {
	1023	+ for (ib = 0; ib < nb; ib++) {
	1024	+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
	1025	+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
	1026	+ }
	1027	+ }
	1028	+ }
	1029	+
	1030	+ for (ib = 0; ib < nb; ib++) {
	1031	+
	1032	+ for (i4 = 0; i4 < 4; i4++)
	1033	+ for (jb = 0; jb < nb; jb++)
	1034	+ w1[ib][i4][jb] = 0.0;
	1035	+
	1036	+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
	1037	+ eig_sort(4, e, indx);
	1038	+
	1039	+ if (e[0] <= 0.0 && 0.0 < e[1]) {
	1040	+
	1041	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	1042	+
	1043	+ if (v > thr) {
	1044	+ for (j4 = 0; j4 < 4; j4++) {
	1045	+ ei2[j4] = 0.0;
	1046	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1047	+ }
	1048	+ for (i4 = 0; i4 < 4; i4++)
	1049	+ for (j4 = 0; j4 < 4; j4++) {
	1050	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1051	+ for (jb = 0; jb < nb; jb++)
	1052	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1053	+ }
	1054	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1055	+ for (i4 = 0; i4 < 4; i4++)
	1056	+ for (jb = 0; jb < nb; jb++)
	1057	+ for (j4 = 0; j4 < 4; j4++)
	1058	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	1059	+ }
	1060	+ }
	1061	+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
	1062	+
	1063	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	1064	+
	1065	+ if (v > thr) {
	1066	+ for (j4 = 0; j4 < 4; j4++) {
	1067	+ ei2[j4] = 0.0;
	1068	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1069	+ }
	1070	+ for (i4 = 0; i4 < 4; i4++)
	1071	+ for (j4 = 0; j4 < 4; j4++) {
	1072	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1073	+ for (jb = 0; jb < nb; jb++)
	1074	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1075	+ }
	1076	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1077	+ for (i4 = 0; i4 < 4; i4++)
	1078	+ for (jb = 0; jb < nb; jb++)
	1079	+ for (j4 = 0; j4 < 4; j4++)
	1080	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	1081	+ }
	1082	+
	1083	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	1084	+
	1085	+ if (v > thr) {
	1086	+ for (j4 = 0; j4 < 4; j4++) {
	1087	+ ei2[j4] = 0.0;
	1088	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1089	+ }
	1090	+ for (i4 = 0; i4 < 4; i4++)
	1091	+ for (j4 = 0; j4 < 4; j4++) {
	1092	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1093	+ for (jb = 0; jb < nb; jb++)
	1094	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1095	+ }
	1096	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1097	+ for (i4 = 0; i4 < 4; i4++)
	1098	+ for (jb = 0; jb < nb; jb++)
	1099	+ for (j4 = 0; j4 < 4; j4++)
	1100	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	1101	+ }
	1102	+
	1103	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	1104	+
	1105	+ if (v > thr) {
	1106	+ for (j4 = 0; j4 < 4; j4++) {
	1107	+ ei2[j4] = 0.0;
	1108	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1109	+ }
	1110	+ for (i4 = 0; i4 < 4; i4++)
	1111	+ for (j4 = 0; j4 < 4; j4++) {
	1112	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1113	+ for (jb = 0; jb < nb; jb++)
	1114	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1115	+ }
	1116	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1117	+ for (i4 = 0; i4 < 4; i4++)
	1118	+ for (jb = 0; jb < nb; jb++)
	1119	+ for (j4 = 0; j4 < 4; j4++)
	1120	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	1121	+ }
	1122	+ }
	1123	+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
	1124	+
	1125	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	1126	+
	1127	+ if (v > thr) {
	1128	+ for (j4 = 0; j4 < 4; j4++) {
	1129	+ ei2[j4] = 0.0;
	1130	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1131	+ }
	1132	+ for (i4 = 0; i4 < 4; i4++)
	1133	+ for (j4 = 0; j4 < 4; j4++) {
	1134	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1135	+ for (jb = 0; jb < nb; jb++)
	1136	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1137	+ }
	1138	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1139	+ for (i4 = 0; i4 < 4; i4++)
	1140	+ for (jb = 0; jb < nb; jb++)
	1141	+ for (j4 = 0; j4 < 4; j4++)
	1142	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	1143	+ }
	1144	+
	1145	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	1146	+
	1147	+ if (v > thr) {
	1148	+ for (j4 = 0; j4 < 4; j4++) {
	1149	+ ei2[j4] = 0.0;
	1150	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1151	+ }
	1152	+ for (i4 = 0; i4 < 4; i4++)
	1153	+ for (j4 = 0; j4 < 4; j4++) {
	1154	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1155	+ for (jb = 0; jb < nb; jb++)
	1156	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1157	+ }
	1158	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1159	+ for (i4 = 0; i4 < 4; i4++)
	1160	+ for (jb = 0; jb < nb; jb++)
	1161	+ for (j4 = 0; j4 < 4; j4++)
	1162	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	1163	+ }
	1164	+
	1165	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	1166	+
	1167	+ if (v > thr) {
	1168	+ for (j4 = 0; j4 < 4; j4++) {
	1169	+ ei2[j4] = 0.0;
	1170	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1171	+ }
	1172	+ for (i4 = 0; i4 < 4; i4++)
	1173	+ for (j4 = 0; j4 < 4; j4++) {
	1174	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1175	+ for (jb = 0; jb < nb; jb++)
	1176	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1177	+ }
	1178	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1179	+ for (i4 = 0; i4 < 4; i4++)
	1180	+ for (jb = 0; jb < nb; jb++)
	1181	+ for (j4 = 0; j4 < 4; j4++)
	1182	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	1183	+ }
	1184	+ }
	1185	+ else if (e[3] <= 0.0) {
	1186	+ for (i4 = 0; i4 < 4; i4++) {
	1187	+ ei2[i4] = ei1[i4][ib];
	1188	+ for (jb = 0; jb < nb; jb++)
	1189	+ ej2[jb][i4] = ej1[i4][jb];
	1190	+ }
	1191	+ libtetrabz_dblstep2(nb, ei2, ej2, w2);
	1192	+ for (i4 = 0; i4 < 4; i4++)
	1193	+ for (jb = 0; jb < nb; jb++)
	1194	+ w1[ib][i4][jb] += w2[jb][i4];
	1195	+ }
	1196	+ else {
	1197	+ continue;
	1198	+ }
	1199	+ }
	1200	+ for (i20 = 0; i20 < 20; i20++)
	1201	+ for (ib = 0; ib < nb; ib++)
	1202	+ for (i4 = 0; i4 < 4; i4++)
	1203	+ for (jb = 0; jb < nb; jb++)
	1204	+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
	1205	+ }
	1206	+ for (ik = 0; ik < nk; ik++)
	1207	+ for (ib = 0; ib < nb; ib++)
	1208	+ for (jb = 0; jb < nb; jb++)
	1209	+ wght[ik][ib][jb] /= (6.0 * (double) nk);
	1210	+
	1211	+ free(ikv[0]);
	1212	+ free(ikv);
	1213	+ free(ei1[0]);
	1214	+ free(ei1);
	1215	+ free(ej1[0]);
	1216	+ free(ej1);
	1217	+ free(ej2[0]);
	1218	+ free(ej2);
	1219	+ free(w1[0][0]);
	1220	+ free(w1[0]);
	1221	+ free(w1);
	1222	+ free(w2[0]);
	1223	+ free(w2);
	1224	+ /*
	1225	+ Convert weight to python list object
	1226	+ */
	1227	+ wght_po = PyList_New(nk * nb * nb);
	1228	+ ierr = 0;
	1229	+ for (ik = 0; ik < nk; ik++)
	1230	+ for (ib = 0; ib < nb; ib++)
	1231	+ for (jb = 0; jb < nb; jb++)
	1232	+ ierr = PyList_SetItem(wght_po, ik * nb * nb + ib * nb + jb, PyFloat_FromDouble(wght[ik][ib][jb]));
	1233	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	1234	+
	1235	+ free(eig1[0]);
	1236	+ free(eig1);
	1237	+ free(eig2[0]);
	1238	+ free(eig2);
	1239	+ free(wght[0][0]);
	1240	+ free(wght[0]);
	1241	+ free(wght);
	1242	+
	1243	+ return wght_po;
	1244	+}
	1245	+/*
	1246	+Compute Dos : Delta(E - E1)
	1247	+*/
	1248	+static PyObject* dos_c(PyObject* self, PyObject* args) {
	1249	+ int it, ik, ib, ii, jj, ie, ** ikv, indx[4], ierr, ng[3], nk, nb, ne;
	1250	+ double wlsm[20][4], ei1, e[4], * w1, v, tsmall[4][3], bvec[3][3], ** eig, * e0, *** wght;
	1251	+ PyObject* eig_po, * e0_po, * wght_po;
	1252	+ /*
	1253	+ Read input from python object
	1254	+ */
	1255	+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOO",
	1256	+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
	1257	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	1258	+ &eig_po, &e0_po))
	1259	+ return NULL;
	1260	+ /*
	1261	+ convert python list object to array
	1262	+ */
	1263	+ eig = (double*)malloc(nk sizeof(double*));
	1264	+ eig[0] = (double)malloc(nk nb * sizeof(double));
	1265	+ wght = (double**)malloc(nk sizeof(double**));
	1266	+ wght[0] = (double*)malloc(nk nb * sizeof(double*));
	1267	+ wght[0][0] = (double)malloc(nk nb * ne * sizeof(double));
	1268	+ for (ik = 0; ik < nk; ik++) {
	1269	+ eig[ik] = eig[0] + ik * nb;
	1270	+ wght[ik] = wght[0] + ik * nb;
	1271	+ for (ib = 0; ib < nb; ib++) {
	1272	+ wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;
	1273	+ }
	1274	+ }
	1275	+ e0 = (double)malloc(ne sizeof(double));
	1276	+
	1277	+ for (ik = 0; ik < nk; ik++)
	1278	+ for (ib = 0; ib < nb; ib++)
	1279	+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
	1280	+ for (ie = 0; ie < ne; ie++)e0[ie] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie));
	1281	+ /*
	1282	+ Start main calculation
	1283	+ */
	1284	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	1285	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	1286	+ for (ik = 0; ik < 6 * nk; ik++) {
	1287	+ ikv[ik] = ikv[0] + ik * 20;
	1288	+ }
	1289	+
	1290	+ ei1 = (double*)malloc(4 sizeof(double*));
	1291	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	1292	+ for (ii = 0; ii < 4; ii++) {
	1293	+ ei1[ii] = ei1[0] + ii * nb;
	1294	+ }
	1295	+
	1296	+ w1 = (double**)malloc(nb sizeof(double**));
	1297	+ w1[0] = (double*)malloc(nb ne * sizeof(double*));
	1298	+ w1[0][0] = (double)malloc(nb ne * 4 * sizeof(double));
	1299	+ for (ib = 0; ib < nb; ib++) {
	1300	+ w1[ib] = w1[0] + ib * ne;
	1301	+ for (ie = 0; ie < ne; ie++) {
	1302	+ w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
	1303	+ }
	1304	+ }
	1305	+
	1306	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	1307	+
	1308	+ for (ik = 0; ik < nk; ik++)
	1309	+ for (ib = 0; ib < nb; ib++)
	1310	+ for (ie = 0; ie < ne; ie++)
	1311	+ wght[ik][ib][ie] = 0.0;
	1312	+
	1313	+ for (it = 0; it < 6 * nk; it++) {
	1314	+
	1315	+ for (ii = 0; ii < 4; ii++)
	1316	+ for (ib = 0; ib < nb; ib++)
	1317	+ ei1[ii][ib] = 0.0;
	1318	+ for (jj = 0; jj < 20; jj++)
	1319	+ for (ii = 0; ii < 4; ii++)
	1320	+ for (ib = 0; ib < nb; ib++)
	1321	+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
	1322	+
	1323	+ for (ib = 0; ib < nb; ib++) {
	1324	+
	1325	+ for (ie = 0; ie < ne; ie++)
	1326	+ for (ii = 0; ii < 4; ii++)
	1327	+ w1[ib][ie][ii] = 0.0;
	1328	+
	1329	+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
	1330	+ eig_sort(4, e, indx);
	1331	+
	1332	+ for (ie = 0; ie < ne; ie++) {
	1333	+
	1334	+ if ((e[0] < e0[ie] && e0[ie] <= e[1]) \|\| (e[0] <= e0[ie] && e0[ie] < e[1])) {
	1335	+
	1336	+ libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
	1337	+ for (ii = 0; ii < 4; ii++)
	1338	+ for (jj = 0; jj < 3; jj++)
	1339	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
	1340	+
	1341	+ }
	1342	+ else if ((e[1] < e0[ie] && e0[ie] <= e[2]) \|\| (e[1] <= e0[ie] && e0[ie] < e[2])) {
	1343	+
	1344	+ libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
	1345	+ for (ii = 0; ii < 4; ii++)
	1346	+ for (jj = 0; jj < 3; jj++)
	1347	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
	1348	+
	1349	+ libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
	1350	+ for (ii = 0; ii < 4; ii++)
	1351	+ for (jj = 0; jj < 3; jj++)
	1352	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
	1353	+ }
	1354	+ else if ((e[2] < e0[ie] && e0[ie] <= e[3]) \|\| (e[2] <= e0[ie] && e0[ie] < e[3])) {
	1355	+
	1356	+ libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
	1357	+ for (ii = 0; ii < 4; ii++)
	1358	+ for (jj = 0; jj < 3; jj++)
	1359	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] / 3.0;
	1360	+ }
	1361	+ else {
	1362	+ continue;
	1363	+ }
	1364	+ }
	1365	+ }
	1366	+ for (ii = 0; ii < 20; ii++)
	1367	+ for (ib = 0; ib < nb; ib++)
	1368	+ for (ie = 0; ie < ne; ie++)
	1369	+ for (jj = 0; jj < 4; jj++)
	1370	+ wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
	1371	+ }
	1372	+ for (ik = 0; ik < nk; ik++)
	1373	+ for (ib = 0; ib < nb; ib++)
	1374	+ for (ie = 0; ie < ne; ie++)
	1375	+ wght[ik][ib][ie] /= (6.0 * (double)nk);
	1376	+
	1377	+ free(ikv[0]);
	1378	+ free(ikv);
	1379	+ free(ei1[0]);
	1380	+ free(ei1);
	1381	+ free(w1[0][0]);
	1382	+ free(w1[0]);
	1383	+ free(w1);
	1384	+ /*
	1385	+ Convert weight to python list object
	1386	+ */
	1387	+ wght_po = PyList_New(nk * nb * ne);
	1388	+ ierr = 0;
	1389	+ for (ik = 0; ik < nk; ik++)
	1390	+ for (ib = 0; ib < nb; ib++)
	1391	+ for (ie = 0; ie < ne; ie++)
	1392	+ ierr = PyList_SetItem(wght_po, ik * nb * ne + ib * ne + ie, PyFloat_FromDouble(wght[ik][ib][ie]));
	1393	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	1394	+
	1395	+ free(eig[0]);
	1396	+ free(eig);
	1397	+ free(wght[0][0]);
	1398	+ free(wght[0]);
	1399	+ free(wght);
	1400	+ free(e0);
	1401	+
	1402	+ return wght_po;
	1403	+}
	1404	+/*
	1405	+Compute integrated Dos : theta(E - E1)
	1406	+*/
	1407	+static PyObject* intdos_c(PyObject* self, PyObject* args) {
	1408	+ int it, ik, ib, ii, jj, ie, ** ikv, indx[4], ierr, ng[3], nk, nb, ne;
	1409	+ double wlsm[20][4], ei1, e[4], * w1, v, tsmall[4][4], bvec[3][3], ** eig, * e0, *** wght;
	1410	+ PyObject* eig_po, * e0_po, * wght_po;
	1411	+ /*
	1412	+ Read input from python object
	1413	+ */
	1414	+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOO",
	1415	+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
	1416	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	1417	+ &eig_po, &e0_po))
	1418	+ return NULL;
	1419	+ /*
	1420	+ convert python list object to array
	1421	+ */
	1422	+ eig = (double*)malloc(nk sizeof(double*));
	1423	+ eig[0] = (double)malloc(nk nb * sizeof(double));
	1424	+ wght = (double**)malloc(nk sizeof(double**));
	1425	+ wght[0] = (double*)malloc(nk nb * sizeof(double*));
	1426	+ wght[0][0] = (double)malloc(nk nb * ne * sizeof(double));
	1427	+ for (ik = 0; ik < nk; ik++) {
	1428	+ eig[ik] = eig[0] + ik * nb;
	1429	+ wght[ik] = wght[0] + ik * nb;
	1430	+ for (ib = 0; ib < nb; ib++) {
	1431	+ wght[ik][ib] = wght[0][0] + ik * nb * ne + ib * ne;
	1432	+ }
	1433	+ }
	1434	+ e0 = (double)malloc(ne sizeof(double));
	1435	+
	1436	+ for (ik = 0; ik < nk; ik++)
	1437	+ for (ib = 0; ib < nb; ib++)
	1438	+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
	1439	+ for (ie = 0; ie < ne; ie++)e0[ie] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie));
	1440	+ /*
	1441	+ Start main calculation
	1442	+ */
	1443	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	1444	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	1445	+ for (ik = 0; ik < 6 * nk; ik++) {
	1446	+ ikv[ik] = ikv[0] + ik * 20;
	1447	+ }
	1448	+
	1449	+ ei1 = (double*)malloc(4 sizeof(double*));
	1450	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	1451	+ for (ii = 0; ii < 4; ii++) {
	1452	+ ei1[ii] = ei1[0] + ii * nb;
	1453	+ }
	1454	+
	1455	+ w1 = (double**)malloc(nb sizeof(double**));
	1456	+ w1[0] = (double*)malloc(nb ne * sizeof(double*));
	1457	+ w1[0][0] = (double)malloc(nb ne * 4 * sizeof(double));
	1458	+ for (ib = 0; ib < nb; ib++) {
	1459	+ w1[ib] = w1[0] + ib * ne;
	1460	+ for (ie = 0; ie < ne; ie++) {
	1461	+ w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
	1462	+ }
	1463	+ }
	1464	+
	1465	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	1466	+
	1467	+ for (ik = 0; ik < nk; ik++)
	1468	+ for (ib = 0; ib < nb; ib++)
	1469	+ for (ie = 0; ie < ne; ie++)
	1470	+ wght[ik][ib][ie] = 0.0;
	1471	+
	1472	+ for (it = 0; it < 6 * nk; it++) {
	1473	+
	1474	+ for (ii = 0; ii < 4; ii++)
	1475	+ for (ib = 0; ib < nb; ib++)
	1476	+ ei1[ii][ib] = 0.0;
	1477	+ for (jj = 0; jj < 20; jj++)
	1478	+ for (ii = 0; ii < 4; ii++)
	1479	+ for (ib = 0; ib < nb; ib++)
	1480	+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
	1481	+
	1482	+ for (ib = 0; ib < nb; ib++)
	1483	+ for (ie = 0; ie < ne; ie++)
	1484	+ for (ii = 0; ii < 4; ii++)
	1485	+ w1[ib][ie][ii] = 0.0;
	1486	+
	1487	+ for (ib = 0; ib < nb; ib++) {
	1488	+
	1489	+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
	1490	+ eig_sort(4, e, indx);
	1491	+
	1492	+ for (ie = 0; ie < ne; ie++) {
	1493	+
	1494	+ if ((e[0] <= e0[ie] && e0[ie] < e[1]) \|\| (e[0] < e0[ie] && e0[ie] <= e[1])) {
	1495	+ libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
	1496	+ for (ii = 0; ii < 4; ii++)
	1497	+ for (jj = 0; jj < 4; jj++)
	1498	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	1499	+ }
	1500	+ else if ((e[1] <= e0[ie] && e0[ie] < e[2]) \|\| (e[1] < e0[ie] && e0[ie] <= e[2])) {
	1501	+
	1502	+ libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
	1503	+ for (ii = 0; ii < 4; ii++)
	1504	+ for (jj = 0; jj < 4; jj++)
	1505	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	1506	+
	1507	+ libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
	1508	+ for (ii = 0; ii < 4; ii++)
	1509	+ for (jj = 0; jj < 4; jj++)
	1510	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	1511	+
	1512	+ libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
	1513	+ for (ii = 0; ii < 4; ii++)
	1514	+ for (jj = 0; jj < 4; jj++)
	1515	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	1516	+
	1517	+ }
	1518	+ else if ((e[2] <= e0[ie] && e0[ie] < e[3]) \|\| (e[2] < e0[ie] && e0[ie] <= e[3])) {
	1519	+
	1520	+ libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
	1521	+ for (ii = 0; ii < 4; ii++)
	1522	+ for (jj = 0; jj < 4; jj++)
	1523	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	1524	+
	1525	+ libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
	1526	+ for (ii = 0; ii < 4; ii++)
	1527	+ for (jj = 0; jj < 4; jj++)
	1528	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	1529	+
	1530	+ libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
	1531	+ for (ii = 0; ii < 4; ii++)
	1532	+ for (jj = 0; jj < 4; jj++)
	1533	+ w1[ib][ie][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	1534	+
	1535	+ }
	1536	+ else if (e[3] <= e0[ie]) {
	1537	+ for (ii = 0; ii < 4; ii++)
	1538	+ w1[ib][ie][ii] = 0.25;
	1539	+ }
	1540	+ else {
	1541	+ continue;
	1542	+ }
	1543	+ }
	1544	+ }
	1545	+
	1546	+ for (ii = 0; ii < 20; ii++)
	1547	+ for (ib = 0; ib < nb; ib++)
	1548	+ for (ie = 0; ie < ne; ie++)
	1549	+ for (jj = 0; jj < 4; jj++)
	1550	+ wght[ikv[it][ii]][ib][ie] += wlsm[ii][jj] * w1[ib][ie][jj];
	1551	+ }
	1552	+ for (ik = 0; ik < nk; ik++)
	1553	+ for (ib = 0; ib < nb; ib++)
	1554	+ for (ie = 0; ie < ne; ie++)
	1555	+ wght[ik][ib][ie] /= (6.0 * (double)nk);
	1556	+
	1557	+ free(ikv[0]);
	1558	+ free(ikv);
	1559	+ free(ei1[0]);
	1560	+ free(ei1);
	1561	+ free(w1[0][0]);
	1562	+ free(w1[0]);
	1563	+ free(w1);
	1564	+ /*
	1565	+ Convert weight to python list object
	1566	+ */
	1567	+ wght_po = PyList_New(nk * nb * ne);
	1568	+ ierr = 0;
	1569	+ for (ik = 0; ik < nk; ik++)
	1570	+ for (ib = 0; ib < nb; ib++)
	1571	+ for (ie = 0; ie < ne; ie++)
	1572	+ ierr = PyList_SetItem(wght_po, ik * nb * ne + ib * ne + ie, PyFloat_FromDouble(wght[ik][ib][ie]));
	1573	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	1574	+
	1575	+ free(eig[0]);
	1576	+ free(eig);
	1577	+ free(wght[0][0]);
	1578	+ free(wght[0]);
	1579	+ free(wght);
	1580	+ free(e0);
	1581	+
	1582	+ return wght_po; }
	1583	+/*
	1584	+ 3rd step for Fermi's golden rule
	1585	+*/
	1586	+static void libtetrabz_fermigr3(
	1587	+ int ne,
	1588	+ double *e0,
	1589	+ double de[4],
	1590	+ double **w1
	1591	+) {
	1592	+ int i4, j3, ie, indx[4];
	1593	+ double e[4], tsmall[4][3], v;
	1594	+
	1595	+ for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
	1596	+ eig_sort(4, e, indx);
	1597	+
	1598	+ for (ie = 0; ie < ne; ie++) {
	1599	+
	1600	+ for (i4 = 0; i4 < 4; i4++) w1[i4][ie] = 0.0;
	1601	+
	1602	+ if (e[0] < e0[ie] && e0[ie] <= e[1]) {
	1603	+
	1604	+ libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
	1605	+ for (i4 = 0; i4 < 4; i4++)
	1606	+ for (j3 = 0; j3 < 3; j3++)
	1607	+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
	1608	+ }
	1609	+ else if (e[1] < e0[ie] && e0[ie] <= e[2]) {
	1610	+
	1611	+ libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
	1612	+ for (i4 = 0; i4 < 4; i4++)
	1613	+ for (j3 = 0; j3 < 3; j3++)
	1614	+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
	1615	+
	1616	+ libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
	1617	+ for (i4 = 0; i4 < 4; i4++)
	1618	+ for (j3 = 0; j3 < 3; j3++)
	1619	+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
	1620	+ }
	1621	+ else if (e[2] < e0[ie] && e0[ie] < e[3]) {
	1622	+
	1623	+ libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
	1624	+ for (i4 = 0; i4 < 4; i4++)
	1625	+ for (j3 = 0; j3 < 3; j3++)
	1626	+ w1[indx[i4]][ie] += v * tsmall[i4][j3] / 3.0;
	1627	+ }
	1628	+ }
	1629	+}
	1630	+/*
	1631	+2nd step for Fermi's golden rule
	1632	+*/
	1633	+static void libtetrabz_fermigr2(
	1634	+ int nb,
	1635	+ int ne,
	1636	+ double *e0,
	1637	+ double *ei1,
	1638	+ double **ej1,
	1639	+ double ***w1
	1640	+) {
	1641	+ int ib, i4, j4, ie, indx[4];
	1642	+ double e[4], tsmall[4][4], v, de[4], thr, **w2;
	1643	+
	1644	+ w2 = (double*)malloc(4 sizeof(double*));
	1645	+ w2[0] = (double)malloc(4 ne * sizeof(double));
	1646	+ for (i4 = 0; i4 < 4; i4++) {
	1647	+ w2[i4] = w2[0] + i4 * ne;
	1648	+ }
	1649	+
	1650	+ thr = 1.0e-8;
	1651	+
	1652	+ for (ib = 0; ib < nb; ib++) {
	1653	+
	1654	+ for (i4 = 0; i4 < 4; i4++)
	1655	+ for (ie = 0; ie < ne; ie++)
	1656	+ w1[ib][i4][ie] = 0.0;
	1657	+
	1658	+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
	1659	+ eig_sort(4, e, indx);
	1660	+
	1661	+ if ((e[0] <= 0.0 && 0.0 < e[1]) \|\| (e[0] < 0.0 && 0.0 <= e[1])) {
	1662	+
	1663	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	1664	+
	1665	+ if (v > thr) {
	1666	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	1667	+ for (i4 = 0; i4 < 4; i4++)
	1668	+ for (j4 = 0; j4 < 4; j4++)
	1669	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	1670	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1671	+ for (i4 = 0; i4 < 4; i4++)
	1672	+ for (j4 = 0; j4 < 4; j4++)
	1673	+ for (ie = 0; ie < ne; ie++)
	1674	+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
	1675	+ }
	1676	+ }
	1677	+ else if ((e[1] <= 0.0 && 0.0 < e[2]) \|\| (e[1] < 0.0 && 0.0 <= e[2])) {
	1678	+
	1679	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	1680	+
	1681	+ if (v > thr) {
	1682	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	1683	+ for (i4 = 0; i4 < 4; i4++)
	1684	+ for (j4 = 0; j4 < 4; j4++)
	1685	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	1686	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1687	+ for (i4 = 0; i4 < 4; i4++)
	1688	+ for (j4 = 0; j4 < 4; j4++)
	1689	+ for (ie = 0; ie < ne; ie++)
	1690	+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
	1691	+ }
	1692	+
	1693	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	1694	+
	1695	+ if (v > thr) {
	1696	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	1697	+ for (i4 = 0; i4 < 4; i4++)
	1698	+ for (j4 = 0; j4 < 4; j4++)
	1699	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	1700	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1701	+ for (i4 = 0; i4 < 4; i4++)
	1702	+ for (j4 = 0; j4 < 4; j4++)
	1703	+ for (ie = 0; ie < ne; ie++)
	1704	+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
	1705	+ }
	1706	+
	1707	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	1708	+
	1709	+ if (v > thr) {
	1710	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	1711	+ for (i4 = 0; i4 < 4; i4++)
	1712	+ for (j4 = 0; j4 < 4; j4++)
	1713	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	1714	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1715	+ for (i4 = 0; i4 < 4; i4++)
	1716	+ for (j4 = 0; j4 < 4; j4++)
	1717	+ for (ie = 0; ie < ne; ie++)
	1718	+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
	1719	+ }
	1720	+ }
	1721	+ else if ((e[2] <= 0.0 && 0.0 < e[3]) \|\| (e[2] < 0.0 && 0.0 <= e[3])) {
	1722	+
	1723	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	1724	+
	1725	+ if (v > thr) {
	1726	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	1727	+ for (i4 = 0; i4 < 4; i4++)
	1728	+ for (j4 = 0; j4 < 4; j4++)
	1729	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	1730	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1731	+ for (i4 = 0; i4 < 4; i4++)
	1732	+ for (j4 = 0; j4 < 4; j4++)
	1733	+ for (ie = 0; ie < ne; ie++)
	1734	+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
	1735	+ }
	1736	+
	1737	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	1738	+
	1739	+ if (v > thr) {
	1740	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	1741	+ for (i4 = 0; i4 < 4; i4++)
	1742	+ for (j4 = 0; j4 < 4; j4++)
	1743	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	1744	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1745	+ for (i4 = 0; i4 < 4; i4++)
	1746	+ for (j4 = 0; j4 < 4; j4++)
	1747	+ for (ie = 0; ie < ne; ie++)
	1748	+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
	1749	+ }
	1750	+
	1751	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	1752	+
	1753	+ if (v > thr) {
	1754	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	1755	+ for (i4 = 0; i4 < 4; i4++)
	1756	+ for (j4 = 0; j4 < 4; j4++)
	1757	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	1758	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1759	+ for (i4 = 0; i4 < 4; i4++)
	1760	+ for (j4 = 0; j4 < 4; j4++)
	1761	+ for (ie = 0; ie < ne; ie++)
	1762	+ w1[ib][indx[i4]][ie] += v * tsmall[i4][j4] * w2[j4][ie];
	1763	+ }
	1764	+ }
	1765	+ else if (e[3] <= 0.0) {
	1766	+ for (i4 = 0; i4 < 4; i4++)
	1767	+ de[i4] = ej1[ib][i4] - ei1[i4];
	1768	+ libtetrabz_fermigr3(ne, e0, de, w2);
	1769	+ for (i4 = 0; i4 < 4; i4++)
	1770	+ for (ie = 0; ie < ne; ie++)
	1771	+ w1[ib][i4][ie] += w2[i4][ie];
	1772	+ }
	1773	+ }
	1774	+
	1775	+ free(w2[0]);
	1776	+ free(w2);
	1777	+}
	1778	+/*
	1779	+Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
	1780	+*/
	1781	+static PyObject* fermigr_c(PyObject* self, PyObject* args)
	1782	+{
	1783	+ int it, ik, ib, i20, i4, j4, jb, ie, **ikv, indx[4], ierr, ng[3], nk, nb, ne;
	1784	+ double wlsm[20][4], ei1, ej1, ei2[4], ej2, e[4], w1, * w2, v, tsmall[4][4], thr,
	1785	+ bvec[3][3], eig1, eig2, e0, ***wght;
	1786	+ PyObject* eig1_po, * eig2_po, e0_po, wght_po;
	1787	+ /*
	1788	+ Read input from python object
	1789	+ */
	1790	+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOOO",
	1791	+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
	1792	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	1793	+ &eig1_po, &eig2_po, &e0_po))
	1794	+ return NULL;
	1795	+ /*
	1796	+ convert python list object to array
	1797	+ */
	1798	+ eig1 = (double*)malloc(nk sizeof(double*));
	1799	+ eig1[0] = (double)malloc(nk nb * sizeof(double));
	1800	+ eig2 = (double*)malloc(nk sizeof(double*));
	1801	+ eig2[0] = (double)malloc(nk nb * sizeof(double));
	1802	+ wght = (double***)malloc(nk sizeof(double***));
	1803	+ wght[0] = (double**)malloc(nk nb * sizeof(double**));
	1804	+ wght[0][0] = (double*)malloc(nk nb * nb * sizeof(double*));
	1805	+ wght[0][0][0] = (double)malloc(nk nb * nb ne sizeof(double));
	1806	+ for (ik = 0; ik < nk; ik++) {
	1807	+ eig1[ik] = eig1[0] + ik * nb;
	1808	+ eig2[ik] = eig2[0] + ik * nb;
	1809	+ wght[ik] = wght[0] + ik * nb;
	1810	+ for (ib = 0; ib < nb; ib++) {
	1811	+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
	1812	+ for (jb = 0; jb < nb; jb++) {
	1813	+ wght[ik][ib][jb] = wght[0][0][0]
	1814	+ + ik * nb * nb * ne + ib * nb * ne + jb * ne;
	1815	+ }
	1816	+ }
	1817	+ }
	1818	+ e0 = (double)malloc(ne sizeof(double));
	1819	+
	1820	+ for (ik = 0; ik < nk; ik++)
	1821	+ for (ib = 0; ib < nb; ib++) {
	1822	+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
	1823	+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
	1824	+ }
	1825	+ for (ie = 0; ie < ne; ie++)e0[ie] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie));
	1826	+ /*
	1827	+ Start main calculation
	1828	+ */
	1829	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	1830	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	1831	+ for (it = 0; it < 6 * nk; it++) {
	1832	+ ikv[it] = ikv[0] + it * 20;
	1833	+ }
	1834	+
	1835	+ ei1 = (double*)malloc(4 sizeof(double*));
	1836	+ ej1 = (double*)malloc(4 sizeof(double*));
	1837	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	1838	+ ej1[0] = (double)malloc(4 nb * sizeof(double));
	1839	+ for (i4 = 0; i4 < 4; i4++) {
	1840	+ ei1[i4] = ei1[0] + i4 * nb;
	1841	+ ej1[i4] = ej1[0] + i4 * nb;
	1842	+ }
	1843	+
	1844	+ ej2 = (double*)malloc(nb sizeof(double*));
	1845	+ ej2[0] = (double)malloc(nb 4 * sizeof(double));
	1846	+ for (ib = 0; ib < nb; ib++)
	1847	+ ej2[ib] = ej2[0] + ib * 4;
	1848	+
	1849	+ w1 = (double***)malloc(nb sizeof(double***));
	1850	+ w1[0] = (double**)malloc(nb 4 * sizeof(double**));
	1851	+ w1[0][0] = (double*)malloc(nb 4 * nb * sizeof(double*));
	1852	+ w1[0][0][0] = (double)malloc(nb 4 * nb * ne * sizeof(double));
	1853	+ for (ib = 0; ib < nb; ib++) {
	1854	+ w1[ib] = w1[0] + ib * 4;
	1855	+ for (i4 = 0; i4 < 4; i4++) {
	1856	+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
	1857	+ for (jb = 0; jb < nb; jb++) {
	1858	+ w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
	1859	+ }
	1860	+ }
	1861	+ }
	1862	+
	1863	+ w2 = (double**)malloc(nb sizeof(double**));
	1864	+ w2[0] = (double*)malloc(nb 4 * sizeof(double*));
	1865	+ w2[0][0] = (double)malloc(nb 4 * ne * sizeof(double));
	1866	+ for (ib = 0; ib < nb; ib++) {
	1867	+ w2[ib] = w2[0] + ib * 4;
	1868	+ for (i4 = 0; i4 < 4; i4++) {
	1869	+ w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
	1870	+ }
	1871	+ }
	1872	+
	1873	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	1874	+
	1875	+ for (ik = 0; ik < nk; ik++)
	1876	+ for (ib = 0; ib < nb; ib++)
	1877	+ for (jb = 0; jb < nb; jb++)
	1878	+ for (ie = 0; ie < ne; ie++)
	1879	+ wght[ik][ib][jb][ie] = 0.0;
	1880	+
	1881	+ thr = 1.0e-10;
	1882	+
	1883	+ for (it = 0; it < 6*nk; it++) {
	1884	+
	1885	+ for (i4 = 0; i4 < 4; i4++)
	1886	+ for (ib = 0; ib < nb; ib++) {
	1887	+ ei1[i4][ib] = 0.0;
	1888	+ ej1[i4][ib] = 0.0;
	1889	+ }
	1890	+ for (i20 = 0; i20 < 20; i20++) {
	1891	+ for (i4 = 0; i4 < 4; i4++) {
	1892	+ for (ib = 0; ib < nb; ib++) {
	1893	+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
	1894	+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
	1895	+ }
	1896	+ }
	1897	+ }
	1898	+
	1899	+ for (ib = 0; ib < nb; ib++) {
	1900	+
	1901	+ for (i4 = 0; i4 < 4; i4++)
	1902	+ for (jb = 0; jb < nb; jb++)
	1903	+ for (ie = 0; ie < ne; ie++)
	1904	+ w1[ib][i4][jb][ie] = 0.0;
	1905	+
	1906	+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
	1907	+ eig_sort(4, e, indx);
	1908	+
	1909	+ if (e[0] <= 0.0 && 0.0 < e[1]) {
	1910	+
	1911	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	1912	+
	1913	+ if (v > thr) {
	1914	+ for (j4 = 0; j4 < 4; j4++) {
	1915	+ ei2[j4] = 0.0;
	1916	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1917	+ }
	1918	+ for (i4 = 0; i4 < 4; i4++)
	1919	+ for (j4 = 0; j4 < 4; j4++) {
	1920	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1921	+ for (jb = 0; jb < nb; jb++)
	1922	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1923	+ }
	1924	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	1925	+ for (i4 = 0; i4 < 4; i4++)
	1926	+ for (jb = 0; jb < nb; jb++)
	1927	+ for (j4 = 0; j4 < 4; j4++)
	1928	+ for (ie = 0; ie < ne; ie++)
	1929	+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
	1930	+ }
	1931	+ }
	1932	+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
	1933	+
	1934	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	1935	+
	1936	+ if (v > thr) {
	1937	+ for (j4 = 0; j4 < 4; j4++) {
	1938	+ ei2[j4] = 0.0;
	1939	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1940	+ }
	1941	+ for (i4 = 0; i4 < 4; i4++)
	1942	+ for (j4 = 0; j4 < 4; j4++) {
	1943	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1944	+ for (jb = 0; jb < nb; jb++)
	1945	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1946	+ }
	1947	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	1948	+ for (i4 = 0; i4 < 4; i4++)
	1949	+ for (jb = 0; jb < nb; jb++)
	1950	+ for (j4 = 0; j4 < 4; j4++)
	1951	+ for (ie = 0; ie < ne; ie++)
	1952	+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
	1953	+ }
	1954	+
	1955	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	1956	+
	1957	+ if (v > thr) {
	1958	+ for (j4 = 0; j4 < 4; j4++) {
	1959	+ ei2[j4] = 0.0;
	1960	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1961	+ }
	1962	+ for (i4 = 0; i4 < 4; i4++)
	1963	+ for (j4 = 0; j4 < 4; j4++) {
	1964	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1965	+ for (jb = 0; jb < nb; jb++)
	1966	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1967	+ }
	1968	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	1969	+ for (i4 = 0; i4 < 4; i4++)
	1970	+ for (jb = 0; jb < nb; jb++)
	1971	+ for (j4 = 0; j4 < 4; j4++)
	1972	+ for (ie = 0; ie < ne; ie++)
	1973	+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
	1974	+ }
	1975	+
	1976	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	1977	+
	1978	+ if (v > thr) {
	1979	+ for (j4 = 0; j4 < 4; j4++) {
	1980	+ ei2[j4] = 0.0;
	1981	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	1982	+ }
	1983	+ for (i4 = 0; i4 < 4; i4++)
	1984	+ for (j4 = 0; j4 < 4; j4++) {
	1985	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	1986	+ for (jb = 0; jb < nb; jb++)
	1987	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	1988	+ }
	1989	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	1990	+ for (i4 = 0; i4 < 4; i4++)
	1991	+ for (jb = 0; jb < nb; jb++)
	1992	+ for (j4 = 0; j4 < 4; j4++)
	1993	+ for (ie = 0; ie < ne; ie++)
	1994	+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
	1995	+ }
	1996	+ }
	1997	+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
	1998	+
	1999	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	2000	+
	2001	+ if (v > thr) {
	2002	+ for (j4 = 0; j4 < 4; j4++) {
	2003	+ ei2[j4] = 0.0;
	2004	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	2005	+ }
	2006	+ for (i4 = 0; i4 < 4; i4++)
	2007	+ for (j4 = 0; j4 < 4; j4++) {
	2008	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	2009	+ for (jb = 0; jb < nb; jb++)
	2010	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	2011	+ }
	2012	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	2013	+ for (i4 = 0; i4 < 4; i4++)
	2014	+ for (jb = 0; jb < nb; jb++)
	2015	+ for (j4 = 0; j4 < 4; j4++)
	2016	+ for (ie = 0; ie < ne; ie++)
	2017	+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
	2018	+ }
	2019	+
	2020	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	2021	+
	2022	+ if (v > thr) {
	2023	+ for (j4 = 0; j4 < 4; j4++) {
	2024	+ ei2[j4] = 0.0;
	2025	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	2026	+ }
	2027	+ for (i4 = 0; i4 < 4; i4++)
	2028	+ for (j4 = 0; j4 < 4; j4++) {
	2029	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	2030	+ for (jb = 0; jb < nb; jb++)
	2031	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	2032	+ }
	2033	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	2034	+ for (i4 = 0; i4 < 4; i4++)
	2035	+ for (jb = 0; jb < nb; jb++)
	2036	+ for (j4 = 0; j4 < 4; j4++)
	2037	+ for (ie = 0; ie < ne; ie++)
	2038	+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
	2039	+ }
	2040	+
	2041	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	2042	+
	2043	+ if (v > thr) {
	2044	+ for (j4 = 0; j4 < 4; j4++) {
	2045	+ ei2[j4] = 0.0;
	2046	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	2047	+ }
	2048	+ for (i4 = 0; i4 < 4; i4++)
	2049	+ for (j4 = 0; j4 < 4; j4++) {
	2050	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	2051	+ for (jb = 0; jb < nb; jb++)
	2052	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	2053	+ }
	2054	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	2055	+ for (i4 = 0; i4 < 4; i4++)
	2056	+ for (jb = 0; jb < nb; jb++)
	2057	+ for (j4 = 0; j4 < 4; j4++)
	2058	+ for (ie = 0; ie < ne; ie++)
	2059	+ w1[ib][indx[i4]][jb][ie] += v * tsmall[i4][j4] * w2[jb][j4][ie];
	2060	+ }
	2061	+ }
	2062	+ else if (e[3] <= 0.0) {
	2063	+ for (i4 = 0; i4 < 4; i4++) {
	2064	+ ei2[i4] = ei1[i4][ib];
	2065	+ for (jb = 0; jb < nb; jb++)
	2066	+ ej2[jb][i4] = ej1[i4][jb];
	2067	+ }
	2068	+ libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
	2069	+ for (i4 = 0; i4 < 4; i4++)
	2070	+ for (jb = 0; jb < nb; jb++)
	2071	+ for (ie = 0; ie < ne; ie++)
	2072	+ w1[ib][i4][jb][ie] += w2[jb][i4][ie];
	2073	+ }
	2074	+ else {
	2075	+ continue;
	2076	+ }
	2077	+ }
	2078	+ for (i20 = 0; i20 < 20; i20++)
	2079	+ for (ib = 0; ib < nb; ib++)
	2080	+ for (i4 = 0; i4 < 4; i4++)
	2081	+ for (jb = 0; jb < nb; jb++)
	2082	+ for (ie = 0; ie < ne; ie++)
	2083	+ wght[ikv[it][i20]][ib][jb][ie] += wlsm[i20][i4] * w1[ib][i4][jb][ie];
	2084	+ }
	2085	+ for (ik = 0; ik < nk; ik++)
	2086	+ for (ib = 0; ib < nb; ib++)
	2087	+ for (jb = 0; jb < nb; jb++)
	2088	+ for (ie = 0; ie < ne; ie++)
	2089	+ wght[ik][ib][jb][ie] /= (6.0 * (double) nk);
	2090	+
	2091	+ free(ikv[0]);
	2092	+ free(ikv);
	2093	+ free(ei1[0]);
	2094	+ free(ei1);
	2095	+ free(ej1[0]);
	2096	+ free(ej1);
	2097	+ free(ej2[0]);
	2098	+ free(ej2);
	2099	+ free(w1[0][0][0]);
	2100	+ free(w1[0][0]);
	2101	+ free(w1[0]);
	2102	+ free(w1);
	2103	+ free(w2[0][0]);
	2104	+ free(w2[0]);
	2105	+ free(w2);
	2106	+ /*
	2107	+ Convert weight to python list object
	2108	+ */
	2109	+ wght_po = PyList_New(nk * nb * nb * ne);
	2110	+ ierr = 0;
	2111	+ for (ik = 0; ik < nk; ik++)
	2112	+ for (ib = 0; ib < nb; ib++)
	2113	+ for (jb = 0; jb < nb; jb++)
	2114	+ for (ie = 0; ie < ne; ie++)
	2115	+ ierr = PyList_SetItem(wght_po, ik * nb * nb * ne + ib * nb * ne + jb * ne + ie,
	2116	+ PyFloat_FromDouble(wght[ik][ib][jb][ie]));
	2117	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	2118	+
	2119	+ free(eig1[0]);
	2120	+ free(eig1);
	2121	+ free(eig2[0]);
	2122	+ free(eig2);
	2123	+ free(wght[0][0][0]);
	2124	+ free(wght[0][0]);
	2125	+ free(wght[0]);
	2126	+ free(wght);
	2127	+ free(e0);
	2128	+
	2129	+ return wght_po;
	2130	+}
	2131	+/*
	2132	+Main SUBROUTINE for occupation : Theta(EF - E1)
	2133	+*/
	2134	+void occ_main(
	2135	+ int nk,
	2136	+ int nb,
	2137	+ double** eig,
	2138	+ int **ikv,
	2139	+ double wlsm[20][4],
	2140	+ double **ei1,
	2141	+ double **w1,
	2142	+ double **wght
	2143	+) {
	2144	+ int it, ik, ib, ii, jj, indx[4];
	2145	+ double e[4], v, tsmall[4][4];
	2146	+
	2147	+ for (ik = 0; ik < nk; ik++)
	2148	+ for (ib = 0; ib < nb; ib++)
	2149	+ wght[ik][ib] = 0.0;
	2150	+
	2151	+ for (it = 0; it < 6 * nk; it++) {
	2152	+
	2153	+ for (ii = 0; ii < 4; ii++)
	2154	+ for (ib = 0; ib < nb; ib++)
	2155	+ ei1[ii][ib] = 0.0;
	2156	+ for (jj = 0; jj < 20; jj++)
	2157	+ for (ii = 0; ii < 4; ii++)
	2158	+ for (ib = 0; ib < nb; ib++)
	2159	+ ei1[ii][ib] += eig[ikv[it][jj]][ib] * wlsm[jj][ii];
	2160	+
	2161	+ for (ib = 0; ib < nb; ib++) {
	2162	+
	2163	+ for (ii = 0; ii < 4; ii++)
	2164	+ w1[ib][ii] = 0.0;
	2165	+
	2166	+ for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
	2167	+ eig_sort(4, e, indx);
	2168	+
	2169	+ if (e[0] <= 0.0 && 0.0 < e[1]) {
	2170	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	2171	+ for (ii = 0; ii < 4; ii++)
	2172	+ for (jj = 0; jj < 4; jj++)
	2173	+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	2174	+ }
	2175	+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
	2176	+
	2177	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	2178	+ for (ii = 0; ii < 4; ii++)
	2179	+ for (jj = 0; jj < 4; jj++)
	2180	+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	2181	+
	2182	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	2183	+ for (ii = 0; ii < 4; ii++)
	2184	+ for (jj = 0; jj < 4; jj++)
	2185	+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	2186	+
	2187	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	2188	+ for (ii = 0; ii < 4; ii++)
	2189	+ for (jj = 0; jj < 4; jj++)
	2190	+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	2191	+ }
	2192	+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
	2193	+
	2194	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	2195	+ for (ii = 0; ii < 4; ii++)
	2196	+ for (jj = 0; jj < 4; jj++)
	2197	+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	2198	+
	2199	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	2200	+ for (ii = 0; ii < 4; ii++)
	2201	+ for (jj = 0; jj < 4; jj++)
	2202	+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	2203	+
	2204	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	2205	+ for (ii = 0; ii < 4; ii++)
	2206	+ for (jj = 0; jj < 4; jj++)
	2207	+ w1[ib][indx[ii]] += v * tsmall[ii][jj] * 0.25;
	2208	+ }
	2209	+ else if (e[3] <= 0.0) {
	2210	+ for (ii = 0; ii < 4; ii++)
	2211	+ w1[ib][ii] += 0.25;
	2212	+ }
	2213	+ else {
	2214	+ continue;
	2215	+ }
	2216	+ }
	2217	+ for (ii = 0; ii < 20; ii++)
	2218	+ for (ib = 0; ib < nb; ib++)
	2219	+ for (jj = 0; jj < 4; jj++) {
	2220	+ wght[ikv[it][ii]][ib] += wlsm[ii][jj] * w1[ib][jj];
	2221	+ }
	2222	+ }
	2223	+ for (ik = 0; ik < nk; ik++)
	2224	+ for (ib = 0; ib < nb; ib++)
	2225	+ wght[ik][ib] /= (6.0 * (double)nk);
	2226	+}
	2227	+/*
	2228	+
	2229	+*/
	2230	+static PyObject* occ_c(PyObject* self, PyObject* args)
	2231	+{
	2232	+ int ik, ib, ii, nk, nb, ng[3], ierr, ** ikv;
	2233	+ double eig, wght, bvec[3][3], wlsm[20][4], ei1, w1;
	2234	+ PyObject* eig_po, * wght_po;
	2235	+ /*
	2236	+ Read input from python object
	2237	+ */
	2238	+ if (!PyArg_ParseTuple(args, "iiiiidddddddddO",
	2239	+ &ng[0], &ng[1], &ng[2], &nk, &nb,
	2240	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	2241	+ &eig_po))
	2242	+ return NULL;
	2243	+ /*
	2244	+ convert python list object to array
	2245	+ */
	2246	+ eig = (double*)malloc(nk sizeof(double*));
	2247	+ eig[0] = (double)malloc(nk nb * sizeof(double));
	2248	+ wght = (double*)malloc(nk sizeof(double*));
	2249	+ wght[0] = (double)malloc(nk nb * sizeof(double));
	2250	+ for (ik = 0; ik < nk; ik++) {
	2251	+ eig[ik] = eig[0] + ik * nb;
	2252	+ wght[ik] = wght[0] + ik * nb;
	2253	+ }
	2254	+
	2255	+ for (ik = 0; ik < nk; ik++)
	2256	+ for (ib = 0; ib < nb; ib++)
	2257	+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
	2258	+
	2259	+ /*
	2260	+ Start main calculation
	2261	+ */
	2262	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	2263	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	2264	+ for (ik = 0; ik < 6 * nk; ik++) {
	2265	+ ikv[ik] = ikv[0] + ik * 20;
	2266	+ }
	2267	+ ei1 = (double*)malloc(4 sizeof(double*));
	2268	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	2269	+ for (ii = 0; ii < 4; ii++) {
	2270	+ ei1[ii] = ei1[0] + ii * nb;
	2271	+ }
	2272	+
	2273	+ w1 = (double*)malloc(nb sizeof(double*));
	2274	+ w1[0] = (double)malloc(nb 4 * sizeof(double));
	2275	+ for (ib = 0; ib < nb; ib++) {
	2276	+ w1[ib] = w1[0] + ib * 4;
	2277	+ }
	2278	+
	2279	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	2280	+
	2281	+ occ_main(nk, nb, eig, ikv, wlsm, ei1, w1, wght);
	2282	+
	2283	+ free(ikv[0]);
	2284	+ free(ikv);
	2285	+ free(ei1[0]);
	2286	+ free(ei1);
	2287	+ free(w1[0]);
	2288	+ free(w1);
	2289	+ /*
	2290	+ Convert weight to python list object
	2291	+ */
	2292	+ wght_po = PyList_New(nk * nb);
	2293	+ ierr = 0;
	2294	+ for (ik = 0; ik < nk; ik++)
	2295	+ for (ib = 0; ib < nb; ib++)
	2296	+ ierr = PyList_SetItem(wght_po, ik * nb + ib, PyFloat_FromDouble(wght[ik][ib]));
	2297	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	2298	+
	2299	+ free(eig[0]);
	2300	+ free(eig);
	2301	+ free(wght[0]);
	2302	+ free(wght);
	2303	+ return wght_po;
	2304	+}
	2305	+/*
	2306	+Calculate Fermi energy
	2307	+*/
	2308	+static PyObject* fermieng_c(PyObject* self, PyObject* args)
	2309	+{
	2310	+ int maxiter, ik, ib, iteration, nk, nb, ng[3],
	2311	+ **ikv, ii, ierr;
	2312	+ double eps, elw, eup, **eig2, sumkmid, bvec[3][3],
	2313	+ eig, wght, wlsm[20][4], ei1, w1, nelec, ef;
	2314	+ PyObject* eig_po, * wght_po;
	2315	+ /*
	2316	+ Read input from python object
	2317	+ */
	2318	+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOd",
	2319	+ &ng[0], &ng[1], &ng[2], &nk, &nb,
	2320	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	2321	+ &eig_po, &nelec))
	2322	+ return NULL;
	2323	+ /*
	2324	+ convert python list object to array
	2325	+ */
	2326	+ eig = (double*)malloc(nk sizeof(double*));
	2327	+ eig[0] = (double)malloc(nk nb * sizeof(double));
	2328	+ wght = (double*)malloc(nk sizeof(double*));
	2329	+ wght[0] = (double)malloc(nk nb * sizeof(double));
	2330	+ for (ik = 0; ik < nk; ik++) {
	2331	+ eig[ik] = eig[0] + ik * nb;
	2332	+ wght[ik] = wght[0] + ik * nb;
	2333	+ }
	2334	+
	2335	+ for (ik = 0; ik < nk; ik++)
	2336	+ for (ib = 0; ib < nb; ib++)
	2337	+ eig[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig_po, ik * nb + ib));
	2338	+
	2339	+ /*
	2340	+ Start main calculation
	2341	+ */
	2342	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	2343	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	2344	+ for (ik = 0; ik < 6 * nk; ik++) {
	2345	+ ikv[ik] = ikv[0] + ik * 20;
	2346	+ }
	2347	+ ei1 = (double*)malloc(4 sizeof(double*));
	2348	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	2349	+ for (ii = 0; ii < 4; ii++) {
	2350	+ ei1[ii] = ei1[0] + ii * nb;
	2351	+ }
	2352	+
	2353	+ w1 = (double*)malloc(nb sizeof(double*));
	2354	+ w1[0] = (double)malloc(nb 4 * sizeof(double));
	2355	+ for (ib = 0; ib < nb; ib++) {
	2356	+ w1[ib] = w1[0] + ib * 4;
	2357	+ }
	2358	+
	2359	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	2360	+
	2361	+ eig2 = (double*)malloc(nk sizeof(double*));
	2362	+ eig2[0] = (double)malloc(nk nb * sizeof(double));
	2363	+ for (ik = 0; ik < nk; ik++) {
	2364	+ eig2[ik] = eig2[0] + ik * nb;
	2365	+ }
	2366	+
	2367	+ maxiter = 300;
	2368	+ eps = 1.0e-10;
	2369	+
	2370	+ elw = eig[0][0];
	2371	+ eup = eig[0][0];
	2372	+ for (ik = 0; ik < nk; ik++) {
	2373	+ for (ib = 0; ib < nb; ib++) {
	2374	+ if (elw > eig[ik][ib]) elw = eig[ik][ib];
	2375	+ if (eup < eig[ik][ib]) eup = eig[ik][ib];
	2376	+ }
	2377	+ }
	2378	+ /*
	2379	+ Bisection method
	2380	+ */
	2381	+ for (iteration = 0; iteration < maxiter; iteration++) {
	2382	+
	2383	+ ef = (eup + elw) * 0.5;
	2384	+ /*
	2385	+ Calc. # of electrons
	2386	+ */
	2387	+ for (ik = 0; ik < nk; ik++)
	2388	+ for (ib = 0; ib < nb; ib++)
	2389	+ eig2[ik][ib] = eig[ik][ib] - ef;
	2390	+ occ_main(nk, nb, eig2, ikv, wlsm, ei1, w1, wght);
	2391	+
	2392	+ sumkmid = 0.0;
	2393	+ for (ik = 0; ik < nk; ik++) {
	2394	+ for (ib = 0; ib < nb; ib++) {
	2395	+ sumkmid += wght[ik][ib];
	2396	+ }
	2397	+ }
	2398	+ /*
	2399	+ convergence check
	2400	+ */
	2401	+ if (fabs(sumkmid - nelec) < eps) {
	2402	+ break;
	2403	+ }
	2404	+ else if (sumkmid < nelec)
	2405	+ elw = ef;
	2406	+ else
	2407	+ eup = ef;
	2408	+ }
	2409	+ free(eig2);
	2410	+ free(ikv[0]);
	2411	+ free(ikv);
	2412	+ free(ei1[0]);
	2413	+ free(ei1);
	2414	+ free(w1[0]);
	2415	+ free(w1);
	2416	+ /*
	2417	+ Convert weight to python list object
	2418	+ */
	2419	+ wght_po = PyList_New(nk * nb+2);
	2420	+ ierr = 0;
	2421	+ for (ik = 0; ik < nk; ik++)
	2422	+ for (ib = 0; ib < nb; ib++)
	2423	+ ierr = PyList_SetItem(wght_po, ik * nb + ib, PyFloat_FromDouble(wght[ik][ib]));
	2424	+ ierr = PyList_SetItem(wght_po, nk*nb, PyFloat_FromDouble(ef));
	2425	+ ierr = PyList_SetItem(wght_po, nk * nb+1, PyLong_FromLong((long)iteration));
	2426	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	2427	+
	2428	+ free(eig[0]);
	2429	+ free(eig);
	2430	+ free(wght[0]);
	2431	+ free(wght);
	2432	+ return wght_po;
	2433	+}
	2434	+/*
	2435	+Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
	2436	+for 0<x<1, 0<y<1-x, 0<z<1-x-y
	2437	+*/
	2438	+/*
	2439	+1, Different each other
	2440	+*/
	2441	+static void libtetrabz_polcmplx_1234(
	2442	+ double g1,
	2443	+ double g2,
	2444	+ double g3,
	2445	+ double g4,
	2446	+ double* w
	2447	+) {
	2448	+ double w2, w3, w4;
	2449	+ /*
	2450	+ Real
	2451	+ */
	2452	+ w2 = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1)) +
	2453	+ (g2 * g2 - 3.0) * g2 * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2454	+ w2 = -2.0 * (g2 * g2 - 1.0) + w2 / (g2 - g1);
	2455	+ w2 = w2 / (g2 - g1);
	2456	+ w3 = 2.0 * (3.0 * g3 * g3 - 1.0) * (atan(g3) - atan(g1)) +
	2457	+ (g3 * g3 - 3.0) * g3 * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2458	+ w3 = -2.0 * (g3 * g3 - 1.0) + w3 / (g3 - g1);
	2459	+ w3 = w3 / (g3 - g1);
	2460	+ w4 = 2.0 * (3.0 * g4 * g4 - 1.0) * (atan(g4) - atan(g1)) +
	2461	+ (g4 * g4 - 3.0) * g4 * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
	2462	+ w4 = -2.0 * (g4 * g4 - 1.0) + w4 / (g4 - g1);
	2463	+ w4 = w4 / (g4 - g1);
	2464	+ w2 = (w2 - w3) / (g2 - g3);
	2465	+ w4 = (w4 - w3) / (g4 - g3);
	2466	+ w[0] = (w4 - w2) / (2.0 * (g4 - g2));
	2467	+ /*
	2468	+ Imaginal
	2469	+ */
	2470	+ w2 = 2.0 * (3.0 - g2 * g2) * g2 * (atan(g2) - atan(g1)) +
	2471	+ (3.0 * g2 * g2 - 1.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2472	+ w2 = 4.0 * g2 - w2 / (g2 - g1);
	2473	+ w2 = w2 / (g2 - g1);
	2474	+ w3 = 2.0 * (3.0 - g3 * g3) * g3 * (atan(g3) - atan(g1)) +
	2475	+ (3.0 * g3 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2476	+ w3 = 4.0 * g3 - w3 / (g3 - g1);
	2477	+ w3 = w3 / (g3 - g1);
	2478	+ w4 = 2.0 * (3.0 - g4 * g4) * g4 * (atan(g4) - atan(g1)) +
	2479	+ (3.0 * g4 * g4 - 1.0) * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
	2480	+ w4 = 4.0 * g4 - w4 / (g4 - g1);
	2481	+ w4 = w4 / (g4 - g1);
	2482	+ w2 = (w2 - w3) / (g2 - g3);
	2483	+ w4 = (w4 - w3) / (g4 - g3);
	2484	+ w[1] = (w4 - w2) / (2.0 * (g4 - g2));
	2485	+}
	2486	+/*
	2487	+2, g4 = g1
	2488	+*/
	2489	+static void libtetrabz_polcmplx_1231(
	2490	+ double g1,
	2491	+ double g2,
	2492	+ double g3,
	2493	+ double* w
	2494	+) {
	2495	+ double w2, w3;
	2496	+ /*
	2497	+ Real
	2498	+ */
	2499	+ w2 = 2.0 * (-1.0 + 3.0 * g2 * g2) * (atan(g2) - atan(g1))
	2500	+ + g2 * (-3.0 + g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2501	+ w2 = 2.0 * (1.0 - g2 * g2) + w2 / (g2 - g1);
	2502	+ w2 = -g1 + w2 / (g2 - g1);
	2503	+ w2 = w2 / (g2 - g1);
	2504	+ w3 = 2.0 * (-1.0 + 3.0 * g3 * g3) * (atan(g3) - atan(g1))
	2505	+ + g3 * (-3.0 + g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2506	+ w3 = 2.0 * (1 - g3 * g3) + w3 / (g3 - g1);
	2507	+ w3 = -g1 + w3 / (g3 - g1);
	2508	+ w3 = w3 / (g3 - g1);
	2509	+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
	2510	+ /*
	2511	+ Imaginal
	2512	+ */
	2513	+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
	2514	+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2515	+ w2 = 4.0 * g2 - w2 / (g2 - g1);
	2516	+ w2 = 1 + w2 / (g2 - g1);
	2517	+ w2 = w2 / (g2 - g1);
	2518	+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
	2519	+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2520	+ w3 = 4.0 * g3 - w3 / (g3 - g1);
	2521	+ w3 = 1 + w3 / (g3 - g1);
	2522	+ w3 = w3 / (g3 - g1);
	2523	+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
	2524	+}
	2525	+/*
	2526	+3, g4 = g3
	2527	+*/
	2528	+static void libtetrabz_polcmplx_1233(
	2529	+ double g1,
	2530	+ double g2,
	2531	+ double g3,
	2532	+ double* w
	2533	+) {
	2534	+ double w2, w3;
	2535	+ /*
	2536	+ Real
	2537	+ */
	2538	+ w2 = 2.0 * (1.0 - 3.0 * g2 * g2) * (atan(g2) - atan(g1)) +
	2539	+ g2 * (3.0 - g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2540	+ w2 = 2.0 * (1 - g2 * g2) - w2 / (g2 - g1);
	2541	+ w2 = w2 / (g2 - g1);
	2542	+ w3 = 2.0 * (1.0 - 3.0 * g3 * g3) * (atan(g3) - atan(g1)) +
	2543	+ g3 * (3.0 - g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2544	+ w3 = 2.0 * (1 - g3 * g3) - w3 / (g3 - g1);
	2545	+ w3 = w3 / (g3 - g1);
	2546	+ w2 = (w3 - w2) / (g3 - g2);
	2547	+ w3 = 4.0 * (1.0 - 3.0 * g1 * g3) * (atan(g3) - atan(g1))
	2548	+ + (3.0 * g1 + 3.0 * g3 - 3.0 * g1 * g3 * g3 + g3 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2549	+ w3 = -4.0 * (1.0 - g1 * g1) + w3 / (g3 - g1);
	2550	+ w3 = 4.0 * g1 + w3 / (g3 - g1);
	2551	+ w3 = w3 / (g3 - g1);
	2552	+ w[0] = (w3 - w2) / (2.0 * (g3 - g2));
	2553	+ /*
	2554	+ Imaginal
	2555	+ */
	2556	+ w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
	2557	+ (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2558	+ w2 = 4.0 * g2 - w2 / (g2 - g1);
	2559	+ w2 = w2 / (g2 - g1);
	2560	+ w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
	2561	+ (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2562	+ w3 = 4.0 * g3 - w3 / (g3 - g1);
	2563	+ w3 = w3 / (g3 - g1);
	2564	+ w2 = (w3 - w2) / (g3 - g2);
	2565	+ w3 = (3.0 * g1 - 3.0 * g1 * g3 * g3 + 3.0 * g3 + g3 * g3 * g3) * (atan(g3) - atan(g1))
	2566	+ + (3.0 * g1 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
	2567	+ w3 = w3 / (g3 - g1) - 4.0 * g1;
	2568	+ w3 = w3 / (g3 - g1) - 2.0;
	2569	+ w3 = (2.0 * w3) / (g3 - g1);
	2570	+ w[1] = (w3 - w2) / (2.0 * (g3 - g2));
	2571	+}
	2572	+/*
	2573	+4, g4 = g1 and g3 = g2
	2574	+*/
	2575	+static void libtetrabz_polcmplx_1221(
	2576	+ double g1,
	2577	+ double g2,
	2578	+ double* w
	2579	+) {
	2580	+ /*
	2581	+ Real
	2582	+ */
	2583	+ w[0] = -2.0 * (-1.0 + 2.0 * g1 * g2 + g2 * g2) * (atan(g2) - atan(g1))
	2584	+ + (g1 + 2.0 * g2 - g1 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2585	+ w[0] = 2.0 * (-1.0 + g1 * g1) + w[0] / (g2 - g1);
	2586	+ w[0] = 3.0 * g1 + w[0] / (g2 - g1);
	2587	+ w[0] = 2.0 + (3.0 * w[0]) / (g2 - g1);
	2588	+ w[0] = w[0] / (2.0 * (g2 - g1));
	2589	+ /*
	2590	+ Imaginal
	2591	+ */
	2592	+ w[1] = 2.0 * (g1 + 2.0 * g2 - g1 * g2 * g2) * (atan(g2) - atan(g1))
	2593	+ + (-1.0 + 2.0 * g1 * g2 + g2 * g2) * log((1 + g2 * g2) / (1 + g1 * g1));
	2594	+ w[1] = -4.0 * g1 + w[1] / (g2 - g1);
	2595	+ w[1] = -3.0 + w[1] / (g2 - g1);
	2596	+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
	2597	+}
	2598	+/*
	2599	+5, g4 = g3 = g2
	2600	+*/
	2601	+static void libtetrabz_polcmplx_1222(
	2602	+ double g1,
	2603	+ double g2,
	2604	+ double* w
	2605	+) {
	2606	+ /*
	2607	+ Real
	2608	+ */
	2609	+ w[0] = 2.0 * (-1.0 + g1 * g1 + 2.0 * g1 * g2) * (atan(g2) - atan(g1))
	2610	+ + (-2.0 * g1 - g2 + g1 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2611	+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
	2612	+ w[0] = g1 - w[0] / (g2 - g1);
	2613	+ w[0] = 1.0 - (3.0 * w[0]) / (g2 - g1);
	2614	+ w[0] = w[0] / (2.0 * (g2 - g1));
	2615	+ /*
	2616	+ Imaginal
	2617	+ */
	2618	+ w[1] = 2.0 * (-2.0 * g1 - g2 + g1 * g1 * g2) * (atan(g2) - atan(g1))
	2619	+ + (1.0 - g1 * g1 - 2.0 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2620	+ w[1] = 4.0 * g1 + w[1] / (g2 - g1);
	2621	+ w[1] = 1.0 + w[1] / (g2 - g1);
	2622	+ w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
	2623	+}
	2624	+/*
	2625	+6, g4 = g3 = g1
	2626	+*/
	2627	+static void libtetrabz_polcmplx_1211(
	2628	+ double g1,
	2629	+ double g2,
	2630	+ double* w
	2631	+) {
	2632	+ /*
	2633	+ Real
	2634	+ */
	2635	+ w[0] = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1))
	2636	+ + g2 * (g2 * g2 - 3.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2637	+ w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
	2638	+ w[0] = -5.0 * g1 + w[0] / (g2 - g1);
	2639	+ w[0] = -11.0 + (3.0 * w[0]) / (g2 - g1);
	2640	+ w[0] = w[0] / (6.0 * (g2 - g1));
	2641	+ /*
	2642	+ Imaginal
	2643	+ */
	2644	+ w[1] = 2.0 * g2 * (-3.0 + g2 * g2) * (atan(g2) - atan(g1))
	2645	+ + (1.0 - 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
	2646	+ w[1] = 4.0 * g2 + w[1] / (g2 - g1);
	2647	+ w[1] = 1.0 + w[1] / (g2 - g1);
	2648	+ w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
	2649	+}
	2650	+/*
	2651	+Tetrahedron method for delta(om - ep + e)
	2652	+*/
	2653	+static void libtetrabz_polcmplx3(
	2654	+ int ne,
	2655	+ double** e0,
	2656	+ double de[4],
	2657	+ double*** w1
	2658	+) {
	2659	+ int i4, ir, indx[4], ie;
	2660	+ double e[4], x[4], thr, w2[4][2], denom;
	2661	+
	2662	+ for (i4 = 0; i4 < 3; i4++) e[i4] = de[i4];
	2663	+ eig_sort(4, e, indx);
	2664	+
	2665	+ for (ie = 0; ie < ne; ie++) {
	2666	+ /*
	2667	+ I am not sure which one is better.
	2668	+ The former is more stable.
	2669	+ The latter is more accurate ?
	2670	+ */
	2671	+ for (i4 = 0; i4 < 4; i4++) {
	2672	+ denom = (de[i4] + e0[ie][0]) * (de[i4] + e0[ie][0]) + e0[ie][1] * e0[ie][1];
	2673	+ w1[i4][ie][0] = 0.25 * (de[i4] + e0[ie][0]) / denom;
	2674	+ w1[i4][ie][1] = 0.25 * (-e0[ie][1]) / denom;
	2675	+ }
	2676	+ continue;
	2677	+
	2678	+ for (i4 = 0; i4 < 4; i4++)
	2679	+ x[i4] = (e[i4] + e0[ie][0]) / e0[ie][1];
	2680	+ /* thr = maxval(de(1:4)) * 1d-3*/
	2681	+ thr = fabs(x[0]);
	2682	+ for (i4 = 0; i4 < 4; i4++)
	2683	+ if (thr < fabs(x[i4])) thr = fabs(x[i4]);
	2684	+ thr = fmax(1.0e-3, thr * 1.0e-2);
	2685	+
	2686	+ if (fabs(x[3] - x[2]) < thr) {
	2687	+ if (fabs(x[3] - x[1]) < thr) {
	2688	+ if (fabs(x[3] - x[0]) < thr) {
	2689	+ /*
	2690	+ e[3] = e[2] = e[1] = e[0]
	2691	+ */
	2692	+ w2[3][0] = 0.25 * x[3] / (1.0 + x[3] * x[3]);
	2693	+ w2[3][1] = 0.25 / (1.0 + x[3] * x[3]);
	2694	+ for (ir = 0; ir < 2; ir++) {
	2695	+ w2[2][ir] = w2[3][ir];
	2696	+ w2[1][ir] = w2[3][ir];
	2697	+ w2[0][ir] = w2[3][ir];
	2698	+ }
	2699	+ }
	2700	+ else {
	2701	+ /*
	2702	+ e[3] = e[2] = e[1]
	2703	+ */
	2704	+ libtetrabz_polcmplx_1211(x[3], x[0], w2[3]);
	2705	+ for (ir = 0; ir < 2; ir++) {
	2706	+ w2[2][ir] = w2[3][ir];
	2707	+ w2[1][ir] = w2[3][ir];
	2708	+ }
	2709	+ libtetrabz_polcmplx_1222(x[0], x[3], w2[0]);
	2710	+ /*
	2711	+ # IF(ANY(w2(1:2,1:4) < 0.0)):
	2712	+ # WRITE(,) ie
	2713	+ # WRITE(*,'(100e15.5)') x[0:4]
	2714	+ # WRITE(*,'(2e15.5)') w2(1:2,1:4)
	2715	+ # STOP "weighting 4=3=2"
	2716	+ */
	2717	+ }
	2718	+ }
	2719	+ else if (fabs(x[1] - x[0]) < thr) {
	2720	+ /*
	2721	+ e[3] = e[2], e[1] = e[0]
	2722	+ */
	2723	+ libtetrabz_polcmplx_1221(x[3], x[1], w2[3]);
	2724	+ for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
	2725	+ libtetrabz_polcmplx_1221(x[1], x[3], w2[1]);
	2726	+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
	2727	+ /*
	2728	+ IF(ANY(w2(1:2,1:4) < 0.0)):
	2729	+ WRITE(,) ie
	2730	+ WRITE(*,'(100e15.5)') x[0:4]
	2731	+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
	2732	+ STOP "weighting 4=3 2=1"
	2733	+ */
	2734	+ }
	2735	+ else {
	2736	+ /*
	2737	+ e[3] = e[2]
	2738	+ */
	2739	+ libtetrabz_polcmplx_1231(x[3], x[0], x[1], w2[3]);
	2740	+ for (ir = 0; ir < 2; ir++) w2[2][ir] = w2[3][ir];
	2741	+ libtetrabz_polcmplx_1233(x[1], x[0], x[3], w2[1]);
	2742	+ libtetrabz_polcmplx_1233(x[0], x[1], x[3], w2[0]);
	2743	+ /*
	2744	+ IF(ANY(w2(1:2,1:4) < 0.0)):
	2745	+ WRITE(,) ie
	2746	+ WRITE(*,'(100e15.5)') x[0:4]
	2747	+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
	2748	+ STOP "weighting 4=3"
	2749	+ */
	2750	+ }
	2751	+ }
	2752	+ else if (fabs(x[2] - x[1]) < thr) {
	2753	+ if (fabs(x[2] - x[0]) < thr) {
	2754	+ /*
	2755	+ e[2] = e[1] = e[0]
	2756	+ */
	2757	+ libtetrabz_polcmplx_1222(x[3], x[2], w2[3]);
	2758	+ libtetrabz_polcmplx_1211(x[2], x[3], w2[2]);
	2759	+ for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
	2760	+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[2][ir];
	2761	+ /*
	2762	+ IF(ANY(w2(1:2,1:4) < 0.0)):
	2763	+ WRITE(,) ie
	2764	+ WRITE(*,'(100e15.5)') x[0:4]
	2765	+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
	2766	+ STOP "weighting 3=2=1"
	2767	+ */
	2768	+ }
	2769	+ else {
	2770	+ /*
	2771	+ e[2] = e[1]
	2772	+ */
	2773	+ libtetrabz_polcmplx_1233(x[3], x[0], x[2], w2[3]);
	2774	+ libtetrabz_polcmplx_1231(x[2], x[0], x[3], w2[2]);
	2775	+ for (ir = 0; ir < 2; ir++) w2[1][ir] = w2[2][ir];
	2776	+ libtetrabz_polcmplx_1233(x[0], x[3], x[2], w2[0]);
	2777	+ /*
	2778	+ IF(ANY(w2(1:2,1:4) < 0.0)):
	2779	+ WRITE(,) ie
	2780	+ WRITE(*,'(100e15.5)') x[0:4]
	2781	+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
	2782	+ STOP "weighting 3=2"
	2783	+ */
	2784	+ }
	2785	+ }
	2786	+ else if (fabs(x[1] - x[0]) < thr) {
	2787	+ /*
	2788	+ e[1] = e[0]
	2789	+ */
	2790	+ libtetrabz_polcmplx_1233(x[3], x[2], x[1], w2[3]);
	2791	+ libtetrabz_polcmplx_1233(x[2], x[3], x[1], w2[2]);
	2792	+ libtetrabz_polcmplx_1231(x[1], x[2], x[3], w2[1]);
	2793	+ for (ir = 0; ir < 2; ir++) w2[0][ir] = w2[1][ir];
	2794	+ /*
	2795	+ IF(ANY(w2(1:2,1:4) < 0.0)):
	2796	+ WRITE(,) ie
	2797	+ WRITE(*,'(100e15.5)') x[0:4]
	2798	+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
	2799	+ STOP "weighting 2=1"
	2800	+ */
	2801	+ }
	2802	+ else {
	2803	+ /*
	2804	+ Different each other.
	2805	+ */
	2806	+ libtetrabz_polcmplx_1234(x[3], x[0], x[1], x[2], w2[3]);
	2807	+ libtetrabz_polcmplx_1234(x[2], x[0], x[1], x[3], w2[2]);
	2808	+ libtetrabz_polcmplx_1234(x[1], x[0], x[2], x[3], w2[1]);
	2809	+ libtetrabz_polcmplx_1234(x[0], x[1], x[2], x[3], w2[0]);
	2810	+ /*
	2811	+ IF(ANY(w2(1:2,1:4) < 0.0)):
	2812	+ WRITE(,) ie
	2813	+ WRITE(*,'(100e15.5)') x[0:4]
	2814	+ WRITE(*,'(2e15.5)') w2(1:2,1:4)
	2815	+ STOP "weighting"
	2816	+ */
	2817	+ }
	2818	+ for (i4 = 0; i4 < 4; i4++) {
	2819	+ w1[indx[i4]][ie][0] = w2[i4][0] / e0[ie][1];
	2820	+ w1[indx[i4]][ie][1] = w2[i4][1] / (-e0[ie][1]);
	2821	+ }
	2822	+ }
	2823	+}
	2824	+/*
	2825	+Tetrahedron method for theta( - E2)
	2826	+*/
	2827	+static void libtetrabz_polcmplx2(
	2828	+ int nb,
	2829	+ int ne,
	2830	+ double* *e0,
	2831	+ double* ei1,
	2832	+ double** ej1,
	2833	+ double**** w1
	2834	+) {
	2835	+ int ib, i4, j4, ir, ie, indx[4];
	2836	+ double e[4], tsmall[4][4], v, de[4], thr, *** w2;
	2837	+
	2838	+ w2 = (double**)malloc(4 sizeof(double**));
	2839	+ w2[0] = (double*)malloc(4 ne * sizeof(double*));
	2840	+ w2[0][0] = (double)malloc(4 ne * 2 * sizeof(double));
	2841	+ for (i4 = 0; i4 < 4; i4++) {
	2842	+ w2[i4] = w2[0] + i4 * ne;
	2843	+ for (ie = 0; ie < ne; ie++) {
	2844	+ w2[i4][ie] = w2[0][0] + i4 * ne * 2 + ie * 2;
	2845	+ }
	2846	+ }
	2847	+
	2848	+ thr = 1.0e-8;
	2849	+
	2850	+ for (ib = 0; ib < nb; ib++) {
	2851	+
	2852	+ for (i4 = 0; i4 < 4; i4++)
	2853	+ for (ie = 0; ie < ne; ie++)
	2854	+ for (ir = 0; ir < 2; ir++)
	2855	+ w1[ib][i4][ie][ir] = 0.0;
	2856	+
	2857	+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
	2858	+ eig_sort(4, e, indx);
	2859	+
	2860	+ if ((e[0] <= 0.0 && 0.0 < e[1]) \|\| (e[0] < 0.0 && 0.0 <= e[1])) {
	2861	+
	2862	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	2863	+
	2864	+ if (v > thr) {
	2865	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	2866	+ for (i4 = 0; i4 < 4; i4++)
	2867	+ for (j4 = 0; j4 < 4; j4++)
	2868	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	2869	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2870	+ for (i4 = 0; i4 < 4; i4++)
	2871	+ for (j4 = 0; j4 < 4; j4++)
	2872	+ for (ie = 0; ie < ne; ie++)
	2873	+ for (ir = 0; ir < 2; ir++)
	2874	+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
	2875	+ }
	2876	+ }
	2877	+ else if ((e[1] <= 0.0 && 0.0 < e[2]) \|\| (e[1] < 0.0 && 0.0 <= e[2])) {
	2878	+
	2879	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	2880	+
	2881	+ if (v > thr) {
	2882	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	2883	+ for (i4 = 0; i4 < 4; i4++)
	2884	+ for (j4 = 0; j4 < 4; j4++)
	2885	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	2886	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2887	+ for (i4 = 0; i4 < 4; i4++)
	2888	+ for (j4 = 0; j4 < 4; j4++)
	2889	+ for (ie = 0; ie < ne; ie++)
	2890	+ for (ir = 0; ir < 2; ir++)
	2891	+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
	2892	+ }
	2893	+
	2894	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	2895	+
	2896	+ if (v > thr) {
	2897	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	2898	+ for (i4 = 0; i4 < 4; i4++)
	2899	+ for (j4 = 0; j4 < 4; j4++)
	2900	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	2901	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2902	+ for (i4 = 0; i4 < 4; i4++)
	2903	+ for (j4 = 0; j4 < 4; j4++)
	2904	+ for (ie = 0; ie < ne; ie++)
	2905	+ for (ir = 0; ir < 2; ir++)
	2906	+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
	2907	+ }
	2908	+
	2909	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	2910	+
	2911	+ if (v > thr) {
	2912	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	2913	+ for (i4 = 0; i4 < 4; i4++)
	2914	+ for (j4 = 0; j4 < 4; j4++)
	2915	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	2916	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2917	+ for (i4 = 0; i4 < 4; i4++)
	2918	+ for (j4 = 0; j4 < 4; j4++)
	2919	+ for (ie = 0; ie < ne; ie++)
	2920	+ for (ir = 0; ir < 2; ir++)
	2921	+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
	2922	+ }
	2923	+ }
	2924	+ else if ((e[2] <= 0.0 && 0.0 < e[3]) \|\| (e[2] < 0.0 && 0.0 <= e[3])) {
	2925	+
	2926	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	2927	+
	2928	+ if (v > thr) {
	2929	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	2930	+ for (i4 = 0; i4 < 4; i4++)
	2931	+ for (j4 = 0; j4 < 4; j4++)
	2932	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	2933	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2934	+ for (i4 = 0; i4 < 4; i4++)
	2935	+ for (j4 = 0; j4 < 4; j4++)
	2936	+ for (ie = 0; ie < ne; ie++)
	2937	+ for (ir = 0; ir < 2; ir++)
	2938	+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
	2939	+ }
	2940	+
	2941	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	2942	+
	2943	+ if (v > thr) {
	2944	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	2945	+ for (i4 = 0; i4 < 4; i4++)
	2946	+ for (j4 = 0; j4 < 4; j4++)
	2947	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	2948	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2949	+ for (i4 = 0; i4 < 4; i4++)
	2950	+ for (j4 = 0; j4 < 4; j4++)
	2951	+ for (ie = 0; ie < ne; ie++)
	2952	+ for (ir = 0; ir < 2; ir++)
	2953	+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
	2954	+ }
	2955	+
	2956	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	2957	+
	2958	+ if (v > thr) {
	2959	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	2960	+ for (i4 = 0; i4 < 4; i4++)
	2961	+ for (j4 = 0; j4 < 4; j4++)
	2962	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	2963	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2964	+ for (i4 = 0; i4 < 4; i4++)
	2965	+ for (j4 = 0; j4 < 4; j4++)
	2966	+ for (ie = 0; ie < ne; ie++)
	2967	+ for (ir = 0; ir < 2; ir++)
	2968	+ w1[ib][indx[i4]][ie][ir] += v * tsmall[i4][j4] * w2[j4][ie][ir];
	2969	+ }
	2970	+ }
	2971	+ else if (e[3] <= 0.0) {
	2972	+ for (i4 = 0; i4 < 4; i4++)
	2973	+ de[i4] = ej1[ib][i4] - ei1[i4];
	2974	+ libtetrabz_polcmplx3(ne, e0, de, w2);
	2975	+ for (i4 = 0; i4 < 4; i4++)
	2976	+ for (ie = 0; ie < ne; ie++)
	2977	+ for (ir = 0; ir < 2; ir++)
	2978	+ w1[ib][i4][ie][ir] += w2[i4][ie][ir];
	2979	+ }
	2980	+ }
	2981	+ free(w2[0][0]);
	2982	+ free(w2[0]);
	2983	+ free(w2);
	2984	+}
	2985	+/*
	2986	+Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
	2987	+*/
	2988	+static PyObject* polcmplx_c(PyObject* self, PyObject* args) {
	2989	+ int it, ik, ie, ib, i4, j4, ir, jb, **ikv, indx[4], i20, ierr, ng[3], nk, nb, ne;
	2990	+ double wlsm[20][4], ei1, ej1, ei2[4], ej2, e[4], * w1, ** w2, v, tsmall[4][4], thr,
	2991	+ bvec[3][3], eig1, eig2, e0, *** wght;
	2992	+ PyObject* eig1_po, * eig2_po, * e0_po, * wght_po;
	2993	+ /*
	2994	+ Read input from python object
	2995	+ */
	2996	+ if (!PyArg_ParseTuple(args, "iiiiiidddddddddOOO",
	2997	+ &ng[0], &ng[1], &ng[2], &nk, &nb, &ne,
	2998	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	2999	+ &eig1_po, &eig2_po, &e0_po))
	3000	+ return NULL;
	3001	+ /*
	3002	+ convert python list object to array
	3003	+ */
	3004	+ eig1 = (double*)malloc(nk sizeof(double*));
	3005	+ eig1[0] = (double)malloc(nk nb * sizeof(double));
	3006	+ eig2 = (double*)malloc(nk sizeof(double*));
	3007	+ eig2[0] = (double)malloc(nk nb * sizeof(double));
	3008	+ wght = (double****)malloc(nk sizeof(double****));
	3009	+ wght[0] = (double***)malloc(nk nb * sizeof(double***));
	3010	+ wght[0][0] = (double**)malloc(nk nb * nb * sizeof(double**));
	3011	+ wght[0][0][0] = (double*)malloc(nk nb * nb * ne * sizeof(double*));
	3012	+ wght[0][0][0][0] = (double)malloc(nk nb * nb * ne * 2 * sizeof(double));
	3013	+ for (ik = 0; ik < nk; ik++) {
	3014	+ eig1[ik] = eig1[0] + ik * nb;
	3015	+ eig2[ik] = eig2[0] + ik * nb;
	3016	+ wght[ik] = wght[0] + ik * nb;
	3017	+ for (ib = 0; ib < nb; ib++) {
	3018	+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
	3019	+ for (jb = 0; jb < nb; jb++) {
	3020	+ wght[ik][ib][jb] = wght[0][0][0]
	3021	+ + ik * nb * nb * ne + ib * nb * ne + jb * ne;
	3022	+ for (ie = 0; ie < ne; ie++) {
	3023	+ wght[ik][ib][jb][ie] = wght[0][0][0][0]
	3024	+ + ik * nb * nb * ne * 2 + ib * nb * ne * 2 + jb * ne * 2 + ie * 2;
	3025	+ }
	3026	+ }
	3027	+ }
	3028	+ }
	3029	+ e0 = (double*)malloc(ne sizeof(double*));
	3030	+ e0[0] = (double)malloc(ne 2 * sizeof(double));
	3031	+ for (ie = 0; ie < ne; ie++)
	3032	+ e0[ie] = e0[0] + ie * 2;
	3033	+
	3034	+ for (ik = 0; ik < nk; ik++)
	3035	+ for (ib = 0; ib < nb; ib++) {
	3036	+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
	3037	+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
	3038	+ }
	3039	+ for (ie = 0; ie < ne; ie++)
	3040	+ for (ir = 0; ir < 2; ir++)
	3041	+ e0[ie][ir] = PyFloat_AsDouble(PyList_GetItem(e0_po, ie * 2 + ir));
	3042	+ /*
	3043	+ Start main calculation
	3044	+ */
	3045	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	3046	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	3047	+ for (it = 0; it < 6 * nk; it++) {
	3048	+ ikv[it] = ikv[0] + it * 20;
	3049	+ }
	3050	+
	3051	+ ei1 = (double*)malloc(4 sizeof(double*));
	3052	+ ej1 = (double*)malloc(4 sizeof(double*));
	3053	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	3054	+ ej1[0] = (double)malloc(4 nb * sizeof(double));
	3055	+ for (i4 = 0; i4 < 4; i4++) {
	3056	+ ei1[i4] = ei1[0] + i4 * nb;
	3057	+ ej1[i4] = ej1[0] + i4 * nb;
	3058	+ }
	3059	+
	3060	+ ej2 = (double*)malloc(nb sizeof(double*));
	3061	+ ej2[0] = (double)malloc(nb 4 * sizeof(double));
	3062	+ for (ib = 0; ib < nb; ib++)
	3063	+ ej2[ib] = ej2[0] + ib * 4;
	3064	+
	3065	+ w1 = (double****)malloc(nb sizeof(double****));
	3066	+ w1[0] = (double***)malloc(nb 4 * sizeof(double***));
	3067	+ w1[0][0] = (double**)malloc(nb 4 * nb * sizeof(double**));
	3068	+ w1[0][0][0] = (double*)malloc(nb 4 * nb * ne * sizeof(double*));
	3069	+ w1[0][0][0][0] = (double)malloc(nb 4 * nb * ne * 2 * sizeof(double));
	3070	+ for (ib = 0; ib < nb; ib++) {
	3071	+ w1[ib] = w1[0] + ib * 4;
	3072	+ for (i4 = 0; i4 < 4; i4++) {
	3073	+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
	3074	+ for (jb = 0; jb < nb; jb++) {
	3075	+ w1[ib][i4][jb] = w1[0][0][0] + ib * 4 * nb * ne + i4 * nb * ne + jb * ne;
	3076	+ for (ie = 0; ie < ne; ie++) {
	3077	+ w1[ib][i4][jb][ie] = w1[0][0][0][0] + ib * 4 * nb * ne * 2 + i4 * nb * ne * 2 + jb * ne * 2 + ie * 2;
	3078	+ }
	3079	+ }
	3080	+ }
	3081	+ }
	3082	+
	3083	+ w2 = (double***)malloc(nb sizeof(double***));
	3084	+ w2[0] = (double**)malloc(nb 4 * sizeof(double**));
	3085	+ w2[0][0] = (double*)malloc(nb 4 * ne * sizeof(double*));
	3086	+ w2[0][0][0] = (double)malloc(nb 4 * ne * 2 * sizeof(double));
	3087	+ for (ib = 0; ib < nb; ib++) {
	3088	+ w2[ib] = w2[0] + ib * 4;
	3089	+ for (i4 = 0; i4 < 4; i4++) {
	3090	+ w2[ib][i4] = w2[0][0] + ib * 4 * ne + i4 * ne;
	3091	+ for (ie = 0; ie < ne; ie++) {
	3092	+ w2[ib][i4][ie] = w2[0][0][0] + ib * 4 * ne * 2 + i4 * ne * 2 + ie * 2;
	3093	+ }
	3094	+ }
	3095	+ }
	3096	+
	3097	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	3098	+
	3099	+ for (ik = 0; ik < nk; ik++)
	3100	+ for (ib = 0; ib < nb; ib++)
	3101	+ for (jb = 0; jb < nb; jb++)
	3102	+ for (ie = 0; ie < ne; ie++)
	3103	+ for (ir = 0; ir < 2; ir++)
	3104	+ wght[ik][ib][jb][ie][ir] = 0.0;
	3105	+
	3106	+ thr = 1.0e-8;
	3107	+
	3108	+ for (it = 0; it < 6 * nk; it++) {
	3109	+
	3110	+ for (i4 = 0; i4 < 4; i4++)
	3111	+ for (ib = 0; ib < nb; ib++) {
	3112	+ ei1[i4][ib] = 0.0;
	3113	+ ej1[i4][ib] = 0.0;
	3114	+ }
	3115	+ for (i20 = 0; i20 < 20; i20++) {
	3116	+ for (i4 = 0; i4 < 4; i4++) {
	3117	+ for (ib = 0; ib < nb; ib++) {
	3118	+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
	3119	+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
	3120	+ }
	3121	+ }
	3122	+ }
	3123	+
	3124	+ for (ib = 0; ib < nb; ib++) {
	3125	+
	3126	+ for (i4 = 0; i4 < 4; i4++)
	3127	+ for (jb = 0; jb < nb; jb++)
	3128	+ for (ie = 0; ie < ne; ie++)
	3129	+ for (ir = 0; ir < 2; ir++)
	3130	+ w1[ib][i4][jb][ie][ir] = 0.0;
	3131	+
	3132	+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
	3133	+ eig_sort(4, e, indx);
	3134	+
	3135	+ if (e[0] <= 0.0 && 0.0 < e[1]) {
	3136	+
	3137	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	3138	+
	3139	+ if (v > thr) {
	3140	+ for (j4 = 0; j4 < 4; j4++) {
	3141	+ ei2[j4] = 0.0;
	3142	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3143	+ }
	3144	+ for (i4 = 0; i4 < 4; i4++)
	3145	+ for (j4 = 0; j4 < 4; j4++) {
	3146	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3147	+ for (jb = 0; jb < nb; jb++)
	3148	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3149	+ }
	3150	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3151	+ for (i4 = 0; i4 < 4; i4++)
	3152	+ for (jb = 0; jb < nb; jb++)
	3153	+ for (j4 = 0; j4 < 4; j4++)
	3154	+ for (ie = 0; ie < ne; ie++)
	3155	+ for (ir = 0; ir < 2; ir++)
	3156	+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
	3157	+ }
	3158	+ }
	3159	+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
	3160	+
	3161	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	3162	+
	3163	+ if (v > thr) {
	3164	+ for (j4 = 0; j4 < 4; j4++) {
	3165	+ ei2[j4] = 0.0;
	3166	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3167	+ }
	3168	+ for (i4 = 0; i4 < 4; i4++)
	3169	+ for (j4 = 0; j4 < 4; j4++) {
	3170	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3171	+ for (jb = 0; jb < nb; jb++)
	3172	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3173	+ }
	3174	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3175	+ for (i4 = 0; i4 < 4; i4++)
	3176	+ for (jb = 0; jb < nb; jb++)
	3177	+ for (j4 = 0; j4 < 4; j4++)
	3178	+ for (ie = 0; ie < ne; ie++)
	3179	+ for (ir = 0; ir < 2; ir++)
	3180	+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
	3181	+ }
	3182	+
	3183	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	3184	+
	3185	+ if (v > thr) {
	3186	+ for (j4 = 0; j4 < 4; j4++) {
	3187	+ ei2[j4] = 0.0;
	3188	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3189	+ }
	3190	+ for (i4 = 0; i4 < 4; i4++)
	3191	+ for (j4 = 0; j4 < 4; j4++) {
	3192	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3193	+ for (jb = 0; jb < nb; jb++)
	3194	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3195	+ }
	3196	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3197	+ for (i4 = 0; i4 < 4; i4++)
	3198	+ for (jb = 0; jb < nb; jb++)
	3199	+ for (j4 = 0; j4 < 4; j4++)
	3200	+ for (ie = 0; ie < ne; ie++)
	3201	+ for (ir = 0; ir < 2; ir++)
	3202	+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
	3203	+ }
	3204	+
	3205	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	3206	+
	3207	+ if (v > thr) {
	3208	+ for (j4 = 0; j4 < 4; j4++) {
	3209	+ ei2[j4] = 0.0;
	3210	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3211	+ }
	3212	+ for (i4 = 0; i4 < 4; i4++)
	3213	+ for (j4 = 0; j4 < 4; j4++) {
	3214	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3215	+ for (jb = 0; jb < nb; jb++)
	3216	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3217	+ }
	3218	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3219	+ for (i4 = 0; i4 < 4; i4++)
	3220	+ for (jb = 0; jb < nb; jb++)
	3221	+ for (j4 = 0; j4 < 4; j4++)
	3222	+ for (ie = 0; ie < ne; ie++)
	3223	+ for (ir = 0; ir < 2; ir++)
	3224	+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
	3225	+ }
	3226	+ }
	3227	+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
	3228	+
	3229	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	3230	+
	3231	+ if (v > thr) {
	3232	+ for (j4 = 0; j4 < 4; j4++) {
	3233	+ ei2[j4] = 0.0;
	3234	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3235	+ }
	3236	+ for (i4 = 0; i4 < 4; i4++)
	3237	+ for (j4 = 0; j4 < 4; j4++) {
	3238	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3239	+ for (jb = 0; jb < nb; jb++)
	3240	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3241	+ }
	3242	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3243	+ for (i4 = 0; i4 < 4; i4++)
	3244	+ for (jb = 0; jb < nb; jb++)
	3245	+ for (j4 = 0; j4 < 4; j4++)
	3246	+ for (ie = 0; ie < ne; ie++)
	3247	+ for (ir = 0; ir < 2; ir++)
	3248	+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
	3249	+ }
	3250	+
	3251	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	3252	+
	3253	+ if (v > thr) {
	3254	+ for (j4 = 0; j4 < 4; j4++) {
	3255	+ ei2[j4] = 0.0;
	3256	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3257	+ }
	3258	+ for (i4 = 0; i4 < 4; i4++)
	3259	+ for (j4 = 0; j4 < 4; j4++) {
	3260	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3261	+ for (jb = 0; jb < nb; jb++)
	3262	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3263	+ }
	3264	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3265	+ for (i4 = 0; i4 < 4; i4++)
	3266	+ for (jb = 0; jb < nb; jb++)
	3267	+ for (j4 = 0; j4 < 4; j4++)
	3268	+ for (ie = 0; ie < ne; ie++)
	3269	+ for (ir = 0; ir < 2; ir++)
	3270	+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
	3271	+ }
	3272	+
	3273	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	3274	+
	3275	+ if (v > thr) {
	3276	+ for (j4 = 0; j4 < 4; j4++) {
	3277	+ ei2[j4] = 0.0;
	3278	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3279	+ }
	3280	+ for (i4 = 0; i4 < 4; i4++)
	3281	+ for (j4 = 0; j4 < 4; j4++) {
	3282	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3283	+ for (jb = 0; jb < nb; jb++)
	3284	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3285	+ }
	3286	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3287	+ for (i4 = 0; i4 < 4; i4++)
	3288	+ for (jb = 0; jb < nb; jb++)
	3289	+ for (j4 = 0; j4 < 4; j4++)
	3290	+ for (ie = 0; ie < ne; ie++)
	3291	+ for (ir = 0; ir < 2; ir++)
	3292	+ w1[ib][indx[i4]][jb][ie][ir] += v * tsmall[i4][j4] * w2[jb][j4][ie][ir];
	3293	+ }
	3294	+ }
	3295	+ else if (e[3] <= 0.0) {
	3296	+ for (i4 = 0; i4 < 4; i4++) {
	3297	+ ei2[i4] = ei1[i4][ib];
	3298	+ for (jb = 0; jb < nb; jb++)
	3299	+ ej2[jb][i4] = ej1[i4][jb];
	3300	+ }
	3301	+ libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
	3302	+ for (i4 = 0; i4 < 4; i4++)
	3303	+ for (jb = 0; jb < nb; jb++)
	3304	+ for (ie = 0; ie < ne; ie++)
	3305	+ for (ir = 0; ir < 2; ir++)
	3306	+ w1[ib][i4][jb][ie][ir] += w2[jb][i4][ie][ir];
	3307	+ }
	3308	+ else {
	3309	+ continue;
	3310	+ }
	3311	+ }
	3312	+ for (i20 = 0; i20 < 20; i20++)
	3313	+ for (ib = 0; ib < nb; ib++)
	3314	+ for (i4 = 0; i4 < 4; i4++)
	3315	+ for (jb = 0; jb < nb; jb++)
	3316	+ for (ie = 0; ie < ne; ie++)
	3317	+ for (ir = 0; ir < 2; ir++)
	3318	+ wght[ikv[it][i20]][ib][jb][ie][ir] += wlsm[i20][i4] * w1[ib][i4][jb][ie][ir];
	3319	+ }
	3320	+ for (ik = 0; ik < nk; ik++)
	3321	+ for (ib = 0; ib < nb; ib++)
	3322	+ for (jb = 0; jb < nb; jb++)
	3323	+ for (ie = 0; ie < ne; ie++)
	3324	+ for (i4 = 0; i4 < 2; i4++)
	3325	+ wght[ik][ib][jb][ie][i4] /= (6.0 * (double) nk);
	3326	+
	3327	+ free(ikv[0]);
	3328	+ free(ikv);
	3329	+ free(ei1[0]);
	3330	+ free(ei1);
	3331	+ free(ej1[0]);
	3332	+ free(ej1);
	3333	+ free(ej2[0]);
	3334	+ free(ej2);
	3335	+ free(w1[0][0][0][0]);
	3336	+ free(w1[0][0][0]);
	3337	+ free(w1[0][0]);
	3338	+ free(w1[0]);
	3339	+ free(w1);
	3340	+ free(w2[0][0][0]);
	3341	+ free(w2[0][0]);
	3342	+ free(w2[0]);
	3343	+ free(w2);
	3344	+ /*
	3345	+ Convert weight to python list object
	3346	+ */
	3347	+ wght_po = PyList_New(nk * nb * nb * ne * 2);
	3348	+ ierr = 0;
	3349	+
	3350	+ for (ik = 0; ik < nk; ik++)
	3351	+ for (ib = 0; ib < nb; ib++)
	3352	+ for (jb = 0; jb < nb; jb++)
	3353	+ for (ie = 0; ie < ne; ie++)
	3354	+ for (ir = 0; ir < 2; ir++)
	3355	+ ierr = PyList_SetItem(wght_po,
	3356	+ ik * nb * nb * ne * 2 + ib * nb * ne * 2 + jb * ne * 2 + ie * 2 + ir,
	3357	+ PyFloat_FromDouble(wght[ik][ib][jb][ie][ir]));
	3358	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	3359	+
	3360	+ free(eig1[0]);
	3361	+ free(eig1);
	3362	+ free(eig2[0]);
	3363	+ free(eig2);
	3364	+ free(wght[0][0][0][0]);
	3365	+ free(wght[0][0][0]);
	3366	+ free(wght[0][0]);
	3367	+ free(wght[0]);
	3368	+ free(wght);
	3369	+ free(e0[0]);
	3370	+ free(e0);
	3371	+
	3372	+ return wght_po;
	3373	+}
	3374	+/*
	3375	+ Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
	3376	+ for 0<x<1, 0<y<1-x, 0<z<1-x-y
	3377	+*/
	3378	+/*
	3379	+1, Different each other
	3380	+*/
	3381	+static double libtetrabz_polstat_1234(
	3382	+ double g1,
	3383	+ double g2,
	3384	+ double g3,
	3385	+ double g4,
	3386	+ double lng1,
	3387	+ double lng2,
	3388	+ double lng3,
	3389	+ double lng4)
	3390	+{
	3391	+ double w2, w3, w4, w;
	3392	+ w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 / (g2 - g1);
	3393	+ w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 / (g3 - g1);
	3394	+ w4 = ((lng4 - lng1) / (g4 - g1) * g4 - 1.0) * g4 / (g4 - g1);
	3395	+ w2 = ((w2 - w3) * g2) / (g2 - g3);
	3396	+ w4 = ((w4 - w3) * g4) / (g4 - g3);
	3397	+ w = (w4 - w2) / (g4 - g2);
	3398	+ return w;
	3399	+}
	3400	+/*
	3401	+2, g4 = g1
	3402	+*/
	3403	+static double libtetrabz_polstat_1231(
	3404	+ double g1,
	3405	+ double g2,
	3406	+ double g3,
	3407	+ double lng1,
	3408	+ double lng2,
	3409	+ double lng3
	3410	+) {
	3411	+ double w2, w3, w;
	3412	+ w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 * g2 / (g2 - g1) - g1 / 2.0;
	3413	+ w2 = w2 / (g2 - g1);
	3414	+ w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 * g3 / (g3 - g1) - g1 / 2.0;
	3415	+ w3 = w3 / (g3 - g1);
	3416	+ w = (w3 - w2) / (g3 - g2);
	3417	+ return w;
	3418	+}
	3419	+/*
	3420	+# 3, g4 = g3
	3421	+*/
	3422	+static double libtetrabz_polstat_1233(
	3423	+ double g1,
	3424	+ double g2,
	3425	+ double g3,
	3426	+ double lng1,
	3427	+ double lng2,
	3428	+ double lng3
	3429	+) {
	3430	+ double w2, w3, w;
	3431	+ w2 = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
	3432	+ w2 = (g2 * w2) / (g2 - g1);
	3433	+ w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
	3434	+ w3 = (g3 * w3) / (g3 - g1);
	3435	+ w2 = (w3 - w2) / (g3 - g2);
	3436	+ w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
	3437	+ w3 = 1.0 - (2.0 * w3 * g1) / (g3 - g1);
	3438	+ w3 = w3 / (g3 - g1);
	3439	+ w = (g3 * w3 - g2 * w2) / (g3 - g2);
	3440	+ return w;
	3441	+}
	3442	+/*
	3443	+4, g4 = g1 and g3 = g2
	3444	+*/
	3445	+static double libtetrabz_polstat_1221(
	3446	+ double g1,
	3447	+ double g2,
	3448	+ double lng1,
	3449	+ double lng2
	3450	+) {
	3451	+ double w;
	3452	+ w = 1.0 - (lng2 - lng1) / (g2 - g1) * g1;
	3453	+ w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
	3454	+ w = -1.0 + (3.0 * g2 * w) / (g2 - g1);
	3455	+ w = w / (2.0 * (g2 - g1));
	3456	+ return w;
	3457	+}
	3458	+/*
	3459	+5, g4 = g3 = g2
	3460	+*/
	3461	+static double libtetrabz_polstat_1222(
	3462	+ double g1,
	3463	+ double g2,
	3464	+ double lng1,
	3465	+ double lng2
	3466	+) {
	3467	+ double w;
	3468	+ w = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
	3469	+ w = (2.0 * g1 * w) / (g2 - g1) - 1.0;
	3470	+ w = (3.0 * g1 * w) / (g2 - g1) + 1.0;
	3471	+ w = w / (2.0 * (g2 - g1));
	3472	+ return w;
	3473	+}
	3474	+/*
	3475	+6, g4 = g3 = g1
	3476	+*/
	3477	+static double libtetrabz_polstat_1211(
	3478	+ double g1,
	3479	+ double g2,
	3480	+ double lng1,
	3481	+ double lng2
	3482	+) {
	3483	+ double w;
	3484	+ w = -1.0 + (lng2 - lng1) / (g2 - g1) * g2;
	3485	+ w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
	3486	+ w = -1.0 + (3.0 * g2 * w) / (2.0 * (g2 - g1));
	3487	+ w = w / (3.0 * (g2 - g1));
	3488	+ return w;
	3489	+}
	3490	+/*
	3491	+Tetrahedron method for delta(om - ep + e)
	3492	+*/
	3493	+static void libtetrabz_polstat3(
	3494	+ double de[4],
	3495	+ double w1[4]
	3496	+)
	3497	+{
	3498	+ int i4, indx[4];
	3499	+ double thr, thr2, e[4], ln[4];
	3500	+
	3501	+ for (i4 = 0; i4 < 4; i4++) e[i4] = de[i4];
	3502	+ eig_sort(4, e, indx);
	3503	+
	3504	+ thr = e[3] * 1.0e-3;
	3505	+ thr2 = 1.0e-8;
	3506	+
	3507	+ for(i4 =0; i4 <4; i4++){
	3508	+ if (e[i4] < thr2) {
	3509	+ if (i4 == 3) {
	3510	+ printf(" Nesting # \n");
	3511	+ }
	3512	+ ln[i4] = 0.0;
	3513	+ e[i4] = 0.0;
	3514	+ }
	3515	+ else{
	3516	+ ln[i4] = log(e[i4]);
	3517	+ }
	3518	+ }
	3519	+
	3520	+ if (fabs(e[3] - e[2]) < thr) {
	3521	+ if (fabs(e[3] - e[1]) < thr) {
	3522	+ if (fabs(e[3] - e[0]) < thr) {
	3523	+ /*
	3524	+ e[3] = e[2] = e[1] = e[0]
	3525	+ */
	3526	+ w1[indx[3]] = 0.25 / e[3];
	3527	+ w1[indx[2]] = w1[indx[3]];
	3528	+ w1[indx[1]] = w1[indx[3]];
	3529	+ w1[indx[0]] = w1[indx[3]];
	3530	+ }
	3531	+ else {
	3532	+ /*
	3533	+ e[3] = e[2] = e[1]
	3534	+ */
	3535	+ w1[indx[3]] = libtetrabz_polstat_1211(e[3], e[0], ln[3], ln[0]);
	3536	+ w1[indx[2]] = w1[indx[3]];
	3537	+ w1[indx[1]] = w1[indx[3]];
	3538	+ w1[indx[0]] = libtetrabz_polstat_1222(e[0], e[3], ln[0], ln[3]);
	3539	+ if (w1[0] < 0.0 \|\| w1[1] < 0.0 \|\| w1[2] < 0.0 \|\| w1[3] < 0.0) {
	3540	+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
	3541	+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
	3542	+ printf("weighting 4=3=2\n");
	3543	+ }
	3544	+ }
	3545	+ }
	3546	+ else if (fabs(e[1] - e[0]) < thr) {
	3547	+ /*
	3548	+ e[3] = e[2], e[1] = e[0]
	3549	+ */
	3550	+ w1[indx[3]] = libtetrabz_polstat_1221(e[3], e[1], ln[3], ln[1]);
	3551	+ w1[indx[2]] = w1[indx[3]];
	3552	+ w1[indx[1]] = libtetrabz_polstat_1221(e[1], e[3], ln[1], ln[3]);
	3553	+ w1[indx[0]] = w1[indx[1]];
	3554	+
	3555	+ if (w1[0] < 0.0 \|\| w1[1] < 0.0 \|\| w1[2] < 0.0 \|\| w1[3] < 0.0) {
	3556	+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
	3557	+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
	3558	+ printf("weighting 4=3 2=1\n");
	3559	+ }
	3560	+ }
	3561	+ else {
	3562	+ /*
	3563	+ e[3] = e[2]
	3564	+ */
	3565	+ w1[indx[3]] = libtetrabz_polstat_1231(e[3], e[0], e[1], ln[3], ln[0], ln[1]);
	3566	+ w1[indx[2]] = w1[indx[3]];
	3567	+ w1[indx[1]] = libtetrabz_polstat_1233(e[1], e[0], e[3], ln[1], ln[0], ln[3]);
	3568	+ w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[1], e[3], ln[0], ln[1], ln[3]);
	3569	+
	3570	+ if (w1[0] < 0.0 \|\| w1[1] < 0.0 \|\| w1[2] < 0.0 \|\| w1[3] < 0.0) {
	3571	+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
	3572	+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
	3573	+ printf("weighting 4=3\n");
	3574	+ }
	3575	+ }
	3576	+ }
	3577	+ else if (fabs(e[2] - e[1]) < thr) {
	3578	+ if (fabs(e[2] - e[0]) < thr) {
	3579	+ /*
	3580	+ e[2] = e[1] = e[0]
	3581	+ */
	3582	+ w1[indx[3]] = libtetrabz_polstat_1222(e[3], e[2], ln[3], ln[2]);
	3583	+ w1[indx[2]] = libtetrabz_polstat_1211(e[2], e[3], ln[2], ln[3]);
	3584	+ w1[indx[1]] = w1[indx[2]];
	3585	+ w1[indx[0]] = w1[indx[2]];
	3586	+
	3587	+ if (w1[0] < 0.0 \|\| w1[1] < 0.0 \|\| w1[2] < 0.0 \|\| w1[3] < 0.0) {
	3588	+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
	3589	+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
	3590	+ printf("weighting 3=2=1\n");
	3591	+ }
	3592	+ }
	3593	+ else {
	3594	+ /*
	3595	+ e[2] = e[1]
	3596	+ */
	3597	+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[0], e[2], ln[3], ln[0], ln[2]);
	3598	+ w1[indx[2]] = libtetrabz_polstat_1231(e[2], e[0], e[3], ln[2], ln[0], ln[3]);
	3599	+ w1[indx[1]] = w1[indx[2]];
	3600	+ w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[3], e[2], ln[0], ln[3], ln[2]);
	3601	+
	3602	+ if (w1[0] < 0.0 \|\| w1[1] < 0.0 \|\| w1[2] < 0.0 \|\| w1[3] < 0.0) {
	3603	+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
	3604	+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
	3605	+ printf("weighting 3=2\n");
	3606	+ }
	3607	+ }
	3608	+ }
	3609	+ else if (fabs(e[1] - e[0]) < thr) {
	3610	+ /*
	3611	+ e[1] = e[0]
	3612	+ */
	3613	+ w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[2], e[1], ln[3], ln[2], ln[1]);
	3614	+ w1[indx[2]] = libtetrabz_polstat_1233(e[2], e[3], e[1], ln[2], ln[3], ln[1]);
	3615	+ w1[indx[1]] = libtetrabz_polstat_1231(e[1], e[2], e[3], ln[1], ln[2], ln[3]);
	3616	+ w1[indx[0]] = w1[indx[1]];
	3617	+
	3618	+ if (w1[0] < 0.0 \|\| w1[1] < 0.0 \|\| w1[2] < 0.0 \|\| w1[3] < 0.0) {
	3619	+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
	3620	+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
	3621	+ printf("weighting 2=1\n");
	3622	+ }
	3623	+ }
	3624	+ else {
	3625	+ /*
	3626	+ Different each other.
	3627	+ */
	3628	+ w1[indx[3]] = libtetrabz_polstat_1234(e[3], e[0], e[1], e[2], ln[3], ln[0], ln[1], ln[2]);
	3629	+ w1[indx[2]] = libtetrabz_polstat_1234(e[2], e[0], e[1], e[3], ln[2], ln[0], ln[1], ln[3]);
	3630	+ w1[indx[1]] = libtetrabz_polstat_1234(e[1], e[0], e[2], e[3], ln[1], ln[0], ln[2], ln[3]);
	3631	+ w1[indx[0]] = libtetrabz_polstat_1234(e[0], e[1], e[2], e[3], ln[0], ln[1], ln[2], ln[3]);
	3632	+
	3633	+ if (w1[0] < 0.0 \|\| w1[1] < 0.0 \|\| w1[2] < 0.0 \|\| w1[3] < 0.0) {
	3634	+ printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
	3635	+ printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
	3636	+ printf("weighting\n");
	3637	+ }
	3638	+ }
	3639	+}
	3640	+/*
	3641	+Tetrahedron method for theta( - E2)
	3642	+*/
	3643	+static void libtetrabz_polstat2(
	3644	+ int nb,
	3645	+ double *ei1,
	3646	+ double **ej1,
	3647	+ double **w1
	3648	+) {
	3649	+ int i4, j4, ib, indx[4];
	3650	+ double v, thr, e[4], de[4], w2[4], tsmall[4][4];
	3651	+
	3652	+ thr = 1.0e-8;
	3653	+
	3654	+ for (ib = 0; ib < nb; ib++) {
	3655	+
	3656	+ for (i4 = 0; i4 < 4; i4++)
	3657	+ w1[ib][i4] = 0.0;
	3658	+
	3659	+ for (i4 = 0; i4 < 4; i4++) e[i4] = -ej1[ib][i4];
	3660	+ eig_sort(4, e, indx);
	3661	+
	3662	+ if ((e[0] <= 0.0 && 0.0 < e[1]) \|\| (e[0] < 0.0 && 0.0 <= e[1])) {
	3663	+
	3664	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	3665	+
	3666	+ if (v > thr) {
	3667	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	3668	+ for (i4 = 0; i4 < 4; i4++)
	3669	+ for (j4 = 0; j4 < 4; j4++)
	3670	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	3671	+ libtetrabz_polstat3(de, w2);
	3672	+ for (i4 = 0; i4 < 4; i4++)
	3673	+ for (j4 = 0; j4 < 4; j4++)
	3674	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
	3675	+ }
	3676	+ }
	3677	+ else if ((e[1] <= 0.0 && 0.0 < e[2]) \|\| (e[1] < 0.0 && 0.0 <= e[2])) {
	3678	+
	3679	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	3680	+
	3681	+ if (v > thr) {
	3682	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	3683	+ for (i4 = 0; i4 < 4; i4++)
	3684	+ for (j4 = 0; j4 < 4; j4++)
	3685	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	3686	+ libtetrabz_polstat3(de, w2);
	3687	+ for (i4 = 0; i4 < 4; i4++)
	3688	+ for (j4 = 0; j4 < 4; j4++)
	3689	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
	3690	+ }
	3691	+
	3692	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	3693	+
	3694	+ if (v > thr) {
	3695	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	3696	+ for (i4 = 0; i4 < 4; i4++)
	3697	+ for (j4 = 0; j4 < 4; j4++)
	3698	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	3699	+ libtetrabz_polstat3(de, w2);
	3700	+ for (i4 = 0; i4 < 4; i4++)
	3701	+ for (j4 = 0; j4 < 4; j4++)
	3702	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
	3703	+ }
	3704	+
	3705	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	3706	+
	3707	+ if (v > thr) {
	3708	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	3709	+ for (i4 = 0; i4 < 4; i4++)
	3710	+ for (j4 = 0; j4 < 4; j4++)
	3711	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	3712	+ libtetrabz_polstat3(de, w2);
	3713	+ for (i4 = 0; i4 < 4; i4++)
	3714	+ for (j4 = 0; j4 < 4; j4++)
	3715	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
	3716	+ }
	3717	+ }
	3718	+ else if ((e[2] <= 0.0 && 0.0 < e[3]) \|\| (e[2] < 0.0 && 0.0 <= e[3])) {
	3719	+
	3720	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	3721	+
	3722	+ if (v > thr) {
	3723	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	3724	+ for (i4 = 0; i4 < 4; i4++)
	3725	+ for (j4 = 0; j4 < 4; j4++)
	3726	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	3727	+ libtetrabz_polstat3(de, w2);
	3728	+ for (i4 = 0; i4 < 4; i4++)
	3729	+ for (j4 = 0; j4 < 4; j4++)
	3730	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
	3731	+ }
	3732	+
	3733	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	3734	+
	3735	+ if (v > thr) {
	3736	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	3737	+ for (i4 = 0; i4 < 4; i4++)
	3738	+ for (j4 = 0; j4 < 4; j4++)
	3739	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	3740	+ libtetrabz_polstat3(de, w2);
	3741	+ for (i4 = 0; i4 < 4; i4++)
	3742	+ for (j4 = 0; j4 < 4; j4++)
	3743	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
	3744	+ }
	3745	+
	3746	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	3747	+
	3748	+ if (v > thr) {
	3749	+ for (j4 = 0; j4 < 4; j4++) de[j4] = 0.0;
	3750	+ for (i4 = 0; i4 < 4; i4++)
	3751	+ for (j4 = 0; j4 < 4; j4++)
	3752	+ de[j4] += (ej1[ib][indx[i4]] - ei1[indx[i4]]) * tsmall[i4][j4];
	3753	+ libtetrabz_polstat3(de, w2);
	3754	+ for (i4 = 0; i4 < 4; i4++)
	3755	+ for (j4 = 0; j4 < 4; j4++)
	3756	+ w1[ib][indx[i4]] += v * tsmall[i4][j4] * w2[j4];
	3757	+ }
	3758	+ }
	3759	+ else if (e[3] <= 0.0) {
	3760	+ for (i4 = 0; i4 < 4; i4++)
	3761	+ de[i4] = ej1[ib][i4] - ei1[i4];
	3762	+ libtetrabz_polstat3(de, w2);
	3763	+ for (i4 = 0; i4 < 4; i4++)
	3764	+ w1[ib][i4] += w2[i4];
	3765	+ }
	3766	+ }
	3767	+}
	3768	+/*
	3769	+Compute Static polarization function
	3770	+*/
	3771	+static PyObject* polstat_c(PyObject* self, PyObject* args)
	3772	+{
	3773	+ int it, ik, ib, jb, i20, i4, j4, **ikv, indx[4], ierr, ng[3], nk, nb;
	3774	+ double wlsm[20][4], ei1, ej1, ei2[4], ej2, e[4], * w1, ** w2, v, tsmall[4][4], thr,
	3775	+ bvec[3][3], eig1, eig2, *** wght;
	3776	+ PyObject* eig1_po, * eig2_po, * wght_po;
	3777	+ /*
	3778	+ Read input from python object
	3779	+ */
	3780	+ if (!PyArg_ParseTuple(args, "iiiiidddddddddOO",
	3781	+ &ng[0], &ng[1], &ng[2], &nk, &nb,
	3782	+ &bvec[0][0], &bvec[0][1], &bvec[0][2], &bvec[1][0], &bvec[1][1], &bvec[1][2], &bvec[2][0], &bvec[2][1], &bvec[2][2],
	3783	+ &eig1_po, &eig2_po))
	3784	+ return NULL;
	3785	+ /*
	3786	+ convert python list object to array
	3787	+ */
	3788	+ eig1 = (double*)malloc(nk sizeof(double*));
	3789	+ eig1[0] = (double)malloc(nk nb * sizeof(double));
	3790	+ eig2 = (double*)malloc(nk sizeof(double*));
	3791	+ eig2[0] = (double)malloc(nk nb * sizeof(double));
	3792	+ wght = (double**)malloc(nk sizeof(double**));
	3793	+ wght[0] = (double*)malloc(nk nb * sizeof(double*));
	3794	+ wght[0][0] = (double)malloc(nk nb * nb * sizeof(double));
	3795	+ for (ik = 0; ik < nk; ik++) {
	3796	+ eig1[ik] = eig1[0] + ik * nb;
	3797	+ eig2[ik] = eig2[0] + ik * nb;
	3798	+ wght[ik] = wght[0] + ik * nb;
	3799	+ for (ib = 0; ib < nb; ib++) {
	3800	+ wght[ik][ib] = wght[0][0] + ik * nb * nb + ib * nb;
	3801	+ }
	3802	+ }
	3803	+
	3804	+ for (ik = 0; ik < nk; ik++)
	3805	+ for (ib = 0; ib < nb; ib++) {
	3806	+ eig1[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig1_po, ik * nb + ib));
	3807	+ eig2[ik][ib] = PyFloat_AsDouble(PyList_GetItem(eig2_po, ik * nb + ib));
	3808	+ }
	3809	+ /*
	3810	+ Start main calculation
	3811	+ */
	3812	+
	3813	+
	3814	+ thr = 1.0e-10;
	3815	+
	3816	+ ikv = (int*)malloc(6 nk * sizeof(int*));
	3817	+ ikv[0] = (int)malloc(6 nk * 20 * sizeof(int));
	3818	+ for (it = 0; it < 6 * nk; it++) {
	3819	+ ikv[it] = ikv[0] + it * 20;
	3820	+ }
	3821	+
	3822	+ ei1 = (double*)malloc(4 sizeof(double*));
	3823	+ ej1 = (double*)malloc(4 sizeof(double*));
	3824	+ ei1[0] = (double)malloc(4 nb * sizeof(double));
	3825	+ ej1[0] = (double)malloc(4 nb * sizeof(double));
	3826	+ for (i4 = 0; i4 < 4; i4++) {
	3827	+ ei1[i4] = ei1[0] + i4 * nb;
	3828	+ ej1[i4] = ej1[0] + i4 * nb;
	3829	+ }
	3830	+
	3831	+ ej2 = (double*)malloc(nb sizeof(double*));
	3832	+ ej2[0] = (double)malloc(nb 4 * sizeof(double));
	3833	+ for (ib = 0; ib < nb; ib++)
	3834	+ ej2[ib] = ej2[0] + ib * 4;
	3835	+
	3836	+ w1 = (double**)malloc(nb sizeof(double**));
	3837	+ w1[0] = (double*)malloc(nb 4 * sizeof(double*));
	3838	+ w1[0][0] = (double)malloc(nb 4 * nb * sizeof(double));
	3839	+ for (ib = 0; ib < nb; ib++) {
	3840	+ w1[ib] = w1[0] + ib * 4;
	3841	+ for (i4 = 0; i4 < 4; i4++) {
	3842	+ w1[ib][i4] = w1[0][0] + ib * 4 * nb + i4 * nb;
	3843	+ }
	3844	+ }
	3845	+
	3846	+ w2 = (double*)malloc(nb sizeof(double*));
	3847	+ w2[0] = (double)malloc(nb 4 * sizeof(double));
	3848	+ for (ib = 0; ib < nb; ib++) {
	3849	+ w2[ib] = w2[0] + ib * 4;
	3850	+ }
	3851	+
	3852	+ libtetrabz_initialize(ng, bvec, wlsm, ikv);
	3853	+
	3854	+ for (ik = 0; ik < nk; ik++)
	3855	+ for (ib = 0; ib < nb; ib++)
	3856	+ for (jb = 0; jb < nb; jb++)
	3857	+ wght[ik][ib][jb] = 0.0;
	3858	+
	3859	+ for(it = 0; it < 6*nk; it++) {
	3860	+
	3861	+ for (i4 = 0; i4 < 4; i4++)
	3862	+ for (ib = 0; ib < nb; ib++) {
	3863	+ ei1[i4][ib] = 0.0;
	3864	+ ej1[i4][ib] = 0.0;
	3865	+ }
	3866	+ for (i20 = 0; i20 < 20; i20++) {
	3867	+ for (i4 = 0; i4 < 4; i4++) {
	3868	+ for (ib = 0; ib < nb; ib++) {
	3869	+ ei1[i4][ib] += eig1[ikv[it][i20]][ib] * wlsm[i20][i4];
	3870	+ ej1[i4][ib] += eig2[ikv[it][i20]][ib] * wlsm[i20][i4];
	3871	+ }
	3872	+ }
	3873	+ }
	3874	+
	3875	+ for (ib = 0; ib < nb; ib++) {
	3876	+
	3877	+ for (i4 = 0; i4 < 4; i4++)
	3878	+ for (jb = 0; jb < nb; jb++)
	3879	+ w1[ib][i4][jb] = 0.0;
	3880	+
	3881	+ for (i4 = 0; i4 < 4; i4++) e[i4] = ei1[i4][ib];
	3882	+ eig_sort(4, e, indx);
	3883	+
	3884	+ if (e[0] <= 0.0 && 0.0 < e[1]) {
	3885	+
	3886	+ libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
	3887	+
	3888	+ if (v > thr) {
	3889	+ for (j4 = 0; j4 < 4; j4++) {
	3890	+ ei2[j4] = 0.0;
	3891	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3892	+ }
	3893	+ for (i4 = 0; i4 < 4; i4++)
	3894	+ for (j4 = 0; j4 < 4; j4++) {
	3895	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3896	+ for (jb = 0; jb < nb; jb++)
	3897	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3898	+ }
	3899	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	3900	+ for (i4 = 0; i4 < 4; i4++)
	3901	+ for (jb = 0; jb < nb; jb++)
	3902	+ for (j4 = 0; j4 < 4; j4++)
	3903	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	3904	+ }
	3905	+ }
	3906	+ else if (e[1] <= 0.0 && 0.0 < e[2]) {
	3907	+
	3908	+ libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
	3909	+
	3910	+ if (v > thr) {
	3911	+ for (j4 = 0; j4 < 4; j4++) {
	3912	+ ei2[j4] = 0.0;
	3913	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3914	+ }
	3915	+ for (i4 = 0; i4 < 4; i4++)
	3916	+ for (j4 = 0; j4 < 4; j4++) {
	3917	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3918	+ for (jb = 0; jb < nb; jb++)
	3919	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3920	+ }
	3921	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	3922	+ for (i4 = 0; i4 < 4; i4++)
	3923	+ for (jb = 0; jb < nb; jb++)
	3924	+ for (j4 = 0; j4 < 4; j4++)
	3925	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	3926	+ }
	3927	+
	3928	+ libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
	3929	+
	3930	+ if (v > thr) {
	3931	+ for (j4 = 0; j4 < 4; j4++) {
	3932	+ ei2[j4] = 0.0;
	3933	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3934	+ }
	3935	+ for (i4 = 0; i4 < 4; i4++)
	3936	+ for (j4 = 0; j4 < 4; j4++) {
	3937	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3938	+ for (jb = 0; jb < nb; jb++)
	3939	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3940	+ }
	3941	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	3942	+ for (i4 = 0; i4 < 4; i4++)
	3943	+ for (jb = 0; jb < nb; jb++)
	3944	+ for (j4 = 0; j4 < 4; j4++)
	3945	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	3946	+ }
	3947	+
	3948	+ libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
	3949	+
	3950	+ if (v > thr) {
	3951	+ for (j4 = 0; j4 < 4; j4++) {
	3952	+ ei2[j4] = 0.0;
	3953	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3954	+ }
	3955	+ for (i4 = 0; i4 < 4; i4++)
	3956	+ for (j4 = 0; j4 < 4; j4++) {
	3957	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3958	+ for (jb = 0; jb < nb; jb++)
	3959	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3960	+ }
	3961	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	3962	+ for (i4 = 0; i4 < 4; i4++)
	3963	+ for (jb = 0; jb < nb; jb++)
	3964	+ for (j4 = 0; j4 < 4; j4++)
	3965	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	3966	+ }
	3967	+ }
	3968	+ else if (e[2] <= 0.0 && 0.0 < e[3]) {
	3969	+
	3970	+ libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
	3971	+
	3972	+ if (v > thr) {
	3973	+ for (j4 = 0; j4 < 4; j4++) {
	3974	+ ei2[j4] = 0.0;
	3975	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3976	+ }
	3977	+ for (i4 = 0; i4 < 4; i4++)
	3978	+ for (j4 = 0; j4 < 4; j4++) {
	3979	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	3980	+ for (jb = 0; jb < nb; jb++)
	3981	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	3982	+ }
	3983	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	3984	+ for (i4 = 0; i4 < 4; i4++)
	3985	+ for (jb = 0; jb < nb; jb++)
	3986	+ for (j4 = 0; j4 < 4; j4++)
	3987	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	3988	+ }
	3989	+
	3990	+ libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
	3991	+
	3992	+ if (v > thr) {
	3993	+ for (j4 = 0; j4 < 4; j4++) {
	3994	+ ei2[j4] = 0.0;
	3995	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	3996	+ }
	3997	+ for (i4 = 0; i4 < 4; i4++)
	3998	+ for (j4 = 0; j4 < 4; j4++) {
	3999	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	4000	+ for (jb = 0; jb < nb; jb++)
	4001	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	4002	+ }
	4003	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	4004	+ for (i4 = 0; i4 < 4; i4++)
	4005	+ for (jb = 0; jb < nb; jb++)
	4006	+ for (j4 = 0; j4 < 4; j4++)
	4007	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	4008	+ }
	4009	+
	4010	+ libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
	4011	+
	4012	+ if (v > thr) {
	4013	+ for (j4 = 0; j4 < 4; j4++) {
	4014	+ ei2[j4] = 0.0;
	4015	+ for (jb = 0; jb < nb; jb++) ej2[jb][j4] = 0.0;
	4016	+ }
	4017	+ for (i4 = 0; i4 < 4; i4++)
	4018	+ for (j4 = 0; j4 < 4; j4++) {
	4019	+ ei2[j4] += ei1[indx[i4]][ib] * tsmall[i4][j4];
	4020	+ for (jb = 0; jb < nb; jb++)
	4021	+ ej2[jb][j4] += ej1[indx[i4]][jb] * tsmall[i4][j4];
	4022	+ }
	4023	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	4024	+ for (i4 = 0; i4 < 4; i4++)
	4025	+ for (jb = 0; jb < nb; jb++)
	4026	+ for (j4 = 0; j4 < 4; j4++)
	4027	+ w1[ib][indx[i4]][jb] += v * tsmall[i4][j4] * w2[jb][j4];
	4028	+ }
	4029	+ }
	4030	+ else if (e[3] <= 0.0) {
	4031	+ for (i4 = 0; i4 < 4; i4++) {
	4032	+ ei2[i4] = ei1[i4][ib];
	4033	+ for (jb = 0; jb < nb; jb++)
	4034	+ ej2[jb][i4] = ej1[i4][jb];
	4035	+ }
	4036	+ libtetrabz_polstat2(nb, ei2, ej2, w2);
	4037	+ for (i4 = 0; i4 < 4; i4++)
	4038	+ for (jb = 0; jb < nb; jb++)
	4039	+ w1[ib][i4][jb] += w2[jb][i4];
	4040	+ }
	4041	+ else {
	4042	+ continue;
	4043	+ }
	4044	+ }
	4045	+ for (i20 = 0; i20 < 20; i20++)
	4046	+ for (ib = 0; ib < nb; ib++)
	4047	+ for (i4 = 0; i4 < 4; i4++)
	4048	+ for (jb = 0; jb < nb; jb++)
	4049	+ wght[ikv[it][i20]][ib][jb] += wlsm[i20][i4] * w1[ib][i4][jb];
	4050	+ }
	4051	+ for (ik = 0; ik < nk; ik++)
	4052	+ for (ib = 0; ib < nb; ib++)
	4053	+ for (jb = 0; jb < nb; jb++)
	4054	+ wght[ik][ib][jb] /= (6.0 * (double) nk);
	4055	+
	4056	+ free(ikv[0]);
	4057	+ free(ikv);
	4058	+ free(ei1[0]);
	4059	+ free(ei1);
	4060	+ free(ej1[0]);
	4061	+ free(ej1);
	4062	+ free(ej2[0]);
	4063	+ free(ej2);
	4064	+ free(w1[0][0]);
	4065	+ free(w1[0]);
	4066	+ free(w1);
	4067	+ free(w2[0]);
	4068	+ free(w2);
	4069	+ /*
	4070	+ Convert weight to python list object
	4071	+ */
	4072	+ wght_po = PyList_New(nk * nb * nb);
	4073	+ ierr = 0;
	4074	+ for (ik = 0; ik < nk; ik++)
	4075	+ for (ib = 0; ib < nb; ib++)
	4076	+ for (jb = 0; jb < nb; jb++)
	4077	+ ierr = PyList_SetItem(wght_po, ik * nb * nb + ib * nb + jb, PyFloat_FromDouble(wght[ik][ib][jb]));
	4078	+ if (ierr != 0)printf("Error in PyList_SetItem\n");
	4079	+
	4080	+ free(eig1[0]);
	4081	+ free(eig1);
	4082	+ free(eig2[0]);
	4083	+ free(eig2);
	4084	+ free(wght[0][0]);
	4085	+ free(wght[0]);
	4086	+ free(wght);
	4087	+
	4088	+ return wght_po;
	4089	+}
	4090	+static PyMethodDef LibTetraBZCMethods[] = {
	4091	+ {"dos_c", dos_c, METH_VARARGS, "Density of States"},
	4092	+ {"intdos_c", intdos_c, METH_VARARGS, "Integrated density of States"},
	4093	+ {"occ_c", occ_c, METH_VARARGS, "Occupations"},
	4094	+ {"fermieng_c", fermieng_c, METH_VARARGS, "Fermi energy"},
	4095	+ {"dbldelta_c", dbldelta_c, METH_VARARGS, "Double delta"},
	4096	+ {"dblstep_c", dblstep_c, METH_VARARGS, "Double step function"},
	4097	+ {"polstat_c", polstat_c, METH_VARARGS, "Static polarization function"},
	4098	+ {"fermigr_c", fermigr_c, METH_VARARGS, "Fermi's golden rule"},
	4099	+ {"polcmplx_c", polcmplx_c, METH_VARARGS, "Dynamical polarization function on complex frequency"},
	4100	+ {NULL, NULL, 0, NULL} /* Sentinel */
	4101	+};
	4102	+
	4103	+static struct PyModuleDef libtetrabzcmodule = {
	4104	+ PyModuleDef_HEAD_INIT,
	4105	+ "libtetrabzc", /* name of module */
	4106	+ NULL, /* module documentation, may be NULL */
	4107	+ -1, /* size of per-interpreter state of the module,
	4108	+ or -1 if the module keeps state in global variables. */
	4109	+ LibTetraBZCMethods
	4110	+};
	4111	+
	4112	+PyMODINIT_FUNC
	4113	+PyInit_libtetrabzc(void)
	4114	+{
	4115	+ return PyModule_Create(&libtetrabzcmodule);
	4116	+}

旧リポジトリブラウザで表示

libtetrabz Forkpythonlibtetrabz

python: コミット

コミットメタ情報

ログメッセージ

変更サマリ

差分

libtetrabz
Fork
python
libtetrabz