OSDN Git Service

Backup
authorMitsuaki Kawamura <kawamitsuaki@gmail.com>
Mon, 14 Mar 2022 17:18:57 +0000 (02:18 +0900)
committerMitsuaki Kawamura <kawamitsuaki@gmail.com>
Mon, 14 Mar 2022 17:18:57 +0000 (02:18 +0900)
src/c/__init__.c [new file with mode: 0644]
src/c/libtetrabz_common.c [new file with mode: 0644]
src/c/libtetrabz_common.h [new file with mode: 0644]
src/c/libtetrabz_dbldelta.c [new file with mode: 0644]
src/c/libtetrabz_dblstep.c [new file with mode: 0644]
src/c/libtetrabz_dos.c [new file with mode: 0644]
src/c/libtetrabz_fermigr.c [new file with mode: 0644]
src/c/libtetrabz_occ.c [new file with mode: 0644]
src/c/libtetrabz_polcmplx.c [new file with mode: 0644]
src/c/libtetrabz_polstat.c [new file with mode: 0644]

diff --git a/src/c/__init__.c b/src/c/__init__.c
new file mode 100644 (file)
index 0000000..35566e0
--- /dev/null
@@ -0,0 +1,31 @@
+//\r
+// Copyright (c) 2014 Mitsuaki Kawamura\r
+//\r
+// Permission is hereby granted, free of charge, to any person obtaining a\r
+// copy of this software and associated documentation files (the\r
+// "Software"), to deal in the Software without restriction, including\r
+// without limitation the rights to use, copy, modify, merge, publish,\r
+// distribute, sublicense, and/or sell copies of the Software, and to\r
+// permit persons to whom the Software is furnished to do so, subject to\r
+// the following conditions:\r
+//\r
+// The above copyright notice and this permission notice shall be included\r
+// in all copies or substantial portions of the Software.\r
+//\r
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS\r
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF\r
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.\r
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY\r
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,\r
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE\r
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
+//\r
+from .libtetrabz_dos import dos\r
+from .libtetrabz_dos import intdos\r
+from .libtetrabz_occ import occ\r
+from .libtetrabz_occ import fermieng\r
+from .libtetrabz_polstat import polstat\r
+from .libtetrabz_fermigr import fermigr\r
+from .libtetrabz_dbldelta import dbldelta\r
+from .libtetrabz_dblstep import dblstep\r
+from .libtetrabz_polcmplx import polcmplx\r
diff --git a/src/c/libtetrabz_common.c b/src/c/libtetrabz_common.c
new file mode 100644 (file)
index 0000000..6e662e8
--- /dev/null
@@ -0,0 +1,581 @@
+/*
+// Copyright (c) 2014 Mitsuaki Kawamura
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+// 
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+
+void libtetrabz_initialize(
+    int *ng,
+    double **bvec,
+    double **wlsm,
+    int **ikv
+    )
+    {
+  /*
+  define shortest diagonal line & define type of tetragonal
+  */
+  int ii, i0, i1, i2, i3, it, itype, ivvec0[4], divvec[4][4], ivvec[6][20][3], ikv0[3], nt;
+  double bvec2[3][3], bvec3[4][3], *bnorm;
+  double wlsm1[4][4] = { {1440.0, 0.0, 30.0, 0.0},
+  {0.0, 1440.0, 0.0, 30.0},
+  {30.0, 0.0, 1440.0, 0.0},
+  {0.0, 30.0, 0.0, 1440.0} },
+  wlsm2[4][4] = {{-38.0, 7.0, 17.0, -28.0},
+  {-28.0, -38.0, 7.0, 17.0},
+  {17.0, -28.0, -38.0, 7.0},
+  {7.0, 17.0, -28.0, -38.0}},
+  wlsm3[4][4] = {{-56.0, 9.0, -46.0, 9.0},
+  {9.0, -56.0, 9.0, -46.0},
+  {-46.0, 9.0, -56.0, 9.0},
+  {9.0, -46.0, 9.0, -56.0}},
+  wlsm4[4][4] = {{-38.0, -28.0, 17.0, 7.0},
+  {7.0, -38.0, -28.0, 17.0},
+  {17.0, 7.0, -38.0, -28.0},
+  {-28.0, 17.0, 7.0, -38.0}},
+  wlsm5[4][4] = {{-18.0, -18.0, 12.0, -18.0},
+  {-18.0, -18.0, -18.0, 12.0},
+  {12.0, -18.0, -18.0, -18.0},
+  {-18.0, 12.0, -18.0, -18.0}};
+
+  for (i1 = 0; i1 < 3; i1++)
+    for (i2 = 0; i2 < 3; i2++)
+      bvec2[i1][i2] = bvec[i1][i2] / (double) ng[i1];
+
+  for (i1 = 0; i1 < 3; i1++) {
+    bvec3[0][i1] = -bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
+    bvec3[1][i1] = bvec2[0][i1] - bvec2[1][i1] + bvec2[2][i1];
+    bvec3[2][i1] = bvec2[0][i1] + bvec2[1][i1] - bvec2[2][i1];
+    bvec3[3][i1] = bvec2[0][i1] + bvec2[1][i1] + bvec2[2][i1];
+  }
+  /*
+  length of delta bvec
+  */
+  for (i1 = 0; i1 < 4; i1++)
+    bnorm[i1] = 0.0;
+  for (i2 = 0; i2 < 3; i2++)
+    bnorm[i1] += bvec3[i1][i2] * bvec3[i1][i2];
+
+  itype = 0;
+  for (i1 = 1; i1 < 4; i1++)
+    if (bnorm[i1] < bnorm[itype]) itype = i1;
+  /*
+  start & last
+  */
+  for (i0 = 0; i0 < 4; i0++) {
+    ivvec0[i0] = 0;
+    for (i1 = 0; i1 < 4; i1++)divvec[i0][i1] = 0;
+    divvec[i0][i0] = 1;
+  }
+  ivvec0[itype] = 1;
+  divvec[itype][itype] = -1;
+  /*
+  Corners of tetrahedron
+  */
+  it = 0;
+  for (i0 = 0; i0 < 3; i0++) {
+    for (i1 = 0; i1 < 3; i1++) {
+      if (i1 == i0) continue;
+      for (i2 = 0; i2 < 3; i2++) {
+        if (i2 == i1 || i2 == i0) continue;
+
+        for (i3 = 0; i3 < 3; i3++) {
+          ivvec[it][0][i3] = ivvec0[i3];
+          ivvec[it][1][i3] = ivvec[it][0][i3] + divvec[i0][i3];
+          ivvec[it][2][i3] = ivvec[it][1][i3] + divvec[i1][i3];
+          ivvec[it][3][i3] = ivvec[it][2][i3] + divvec[i2][i3];
+        }
+
+        it += 1;
+      }
+    }
+  }
+  /*
+  Additional points
+  */
+  for (i1 = 0; i1 < 6; i1++) {
+    for (i2 = 0; i2 < 3; i2++) {
+      ivvec[i1][4][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][1][i2];
+      ivvec[i1][5][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][2][i2];
+      ivvec[i1][6][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][3][i2];
+      ivvec[i1][7][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][0][i2];
+
+      ivvec[i1][8][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][2][i2];
+      ivvec[i1][9][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][3][i2];
+      ivvec[i1][10][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][0][i2];
+      ivvec[i1][11][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][1][i2];
+
+      ivvec[i1][12][i2] = 2 * ivvec[i1][0][i2] - ivvec[i1][3][i2];
+      ivvec[i1][13][i2] = 2 * ivvec[i1][1][i2] - ivvec[i1][0][i2];
+      ivvec[i1][14][i2] = 2 * ivvec[i1][2][i2] - ivvec[i1][1][i2];
+      ivvec[i1][15][i2] = 2 * ivvec[i1][3][i2] - ivvec[i1][2][i2];
+
+      ivvec[i1][16][i2] = ivvec[i1][3][i2] - ivvec[i1][0][i2] + ivvec[i1][1][i2];
+      ivvec[i1][17][i2] = ivvec[i1][0][i2] - ivvec[i1][1][i2] + ivvec[i1][2][i2];
+      ivvec[i1][18][i2] = ivvec[i1][1][i2] - ivvec[i1][2][i2] + ivvec[i1][3][i2];
+      ivvec[i1][19][i2] = ivvec[i1][2][i2] - ivvec[i1][3][i2] + ivvec[i1][0][i2];
+    }
+  }
+
+  for (i1 = 0; i1 < 4; i1++) {
+    for (i2 = 0; i2 < 4; i2++) {
+      wlsm[i1][i2] = wlsm1[i1][i2] /= 1260.0;
+      wlsm[i1][i2 + 4] = wlsm2[i1][i2] /= 1260.0;
+      wlsm[i1][i2 + 8] = wlsm3[i1][i2] /= 1260.0;
+      wlsm[i1][i2 + 12] = wlsm4[i1][i2] /= 1260.0;
+      wlsm[i1][i2 + 16] = wlsm5[i1][i2] /= 1260.0;
+    }
+  }
+  /*
+  k-index for energy
+  */
+  nt = 0;
+  for (i2 = 0; i2 < ng[2]; i2++) {
+    for (i1 = 0; i1 < ng[1]; i1++) {
+      for (i0 = 0; i0 < ng[0]; i0++) {
+
+        for (it = 0; it < 6; it++) {
+
+          for (ii = 0; ii < 20; ii++) {
+            ikv0[0] = (i0 + ivvec[it][ii][0]) % ng[0];
+            ikv0[1] = (i1 + ivvec[it][ii][1]) % ng[1];
+            ikv0[2] = (i2 + ivvec[it][ii][2]) % ng[2];
+            ikv[nt][ii] = ikv0[2] + ng[2] * ikv0[1] + ng[2] * ng[1] * ikv0[0];
+          }
+          nt += 1;
+        }
+      }
+    }
+  }
+}
+
+
+void libtetrabz_tsmall_a1(
+    double *e,
+    double e0,
+    double *v,
+    double **tsmall
+) {
+  /*
+  Cut small tetrahedron A1
+  */
+  double a10, a20, a30;
+  a10 = (e0 - e[0]) / (e[1] - e[0]);
+  a20 = (e0 - e[0]) / (e[2] - e[0]);
+  a30 = (e0 - e[0]) / (e[3] - e[0]);
+
+  *v = a10 * a20 * a30;
+
+  tsmall[0][0] = 1.0;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 1.0 - a10;
+  tsmall[1][1] = a10;
+  tsmall[1][2] = 0.0;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 1.0 - a20;
+  tsmall[2][1] = 0.0;
+  tsmall[2][2] = a20;
+  tsmall[2][3] = 0.0;
+  tsmall[3][0] = 1.0 - a30;
+  tsmall[3][1] = 0.0;
+  tsmall[3][2] = 0.0;
+  tsmall[3][3] = a30;
+}
+
+
+void libtetrabz_tsmall_b1(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+  /*
+  Cut small tetrahedron B1
+  :param e: numpy
+  :return:
+  */
+  double a13, a20, a30;
+  a13 = (e0 - e[3]) / (e[1] - e[3]);
+  a20 = (e0 - e[0]) / (e[2] - e[0]);
+  a30 = (e0 - e[0]) / (e[3] - e[0]);
+
+  *v = a20 * a30 * a13;
+
+  tsmall[0][0] = 1.0;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 1.0 - a20;
+  tsmall[1][1] = 0.0;
+  tsmall[1][2] = a20;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 1.0 - a30;
+  tsmall[2][1] = 0.0;
+  tsmall[2][2] = 0.0;
+  tsmall[2][3] = a30;
+  tsmall[3][0] = 0.0;
+  tsmall[3][1] = a13;
+  tsmall[3][2] = 0.0;
+  tsmall[3][3] = 1.0 - a13;
+}
+
+
+void libtetrabz_tsmall_b2(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+  /*
+  Cut small tetrahedron B2
+  :param e:
+  :return:
+  */
+  double a21, a31;
+  a21 = (e0 - e[1]) / (e[2] - e[1]);
+  a31 = (e0 - e[1]) / (e[3] - e[1]);
+
+  *v = a21 * a31;
+
+  tsmall[0][0] = 1.0;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 0.0;
+  tsmall[1][1] = 1.0;
+  tsmall[1][2] = 0.0;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = 1.0 - a21;
+  tsmall[2][2] = a21;
+  tsmall[2][3] = 0.0;
+  tsmall[3][0] = 0.0;
+  tsmall[3][1] = 1.0 - a31;
+  tsmall[3][2] = 0.0;
+  tsmall[3][3] = a31;
+}
+
+void libtetrabz_tsmall_b3(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+  /*
+  Cut small tetrahedron B3
+  :param e:
+  :return:
+  */
+  double a12, a20, a31;
+  a12 = (e0 - e[2]) / (e[1] - e[2]);
+  a20 = (e0 - e[0]) / (e[2] - e[0]);
+  a31 = (e0 - e[1]) / (e[3] - e[1]);
+
+  *v = a12 * a20 * a31;
+
+  tsmall[0][0] = 1.0;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 1.0 - a20;
+  tsmall[1][1] = 0.0;
+  tsmall[1][2] = a20;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = a12;
+  tsmall[2][2] = 1.0 - a12;
+  tsmall[2][3] = 0.0;
+  tsmall[3][0] = 0.0;
+  tsmall[3][1] = 1.0 - a31;
+  tsmall[3][2] = 0.0;
+  tsmall[3][3] = a31;
+}
+
+
+void libtetrabz_tsmall_c1(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+  /*
+  Cut small tetrahedron C1
+  :param e:
+  :return:
+  */
+  double a32;
+  a32 = (e0 - e[2]) / (e[3] - e[2]);
+
+  *v = a32;
+
+  tsmall[0][0] = 1.0;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 0.0;
+  tsmall[1][1] = 1.0;
+  tsmall[1][2] = 0.0;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = 0.0;
+  tsmall[2][2] = 1.0;
+  tsmall[2][3] = 0.0;
+  tsmall[3][0] = 0.0;
+  tsmall[3][1] = 0.0;
+  tsmall[3][2] = 1.0 - a32;
+  tsmall[3][3] = a32;
+}
+
+
+void libtetrabz_tsmall_c2(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+  /*
+  Cut small tetrahedron C2
+  :param e:
+  :return:
+  */
+  double a23, a31;
+  a23 = (e0 - e[3]) / (e[2] - e[3]);
+  a31 = (e0 - e[1]) / (e[3] - e[1]);
+
+  *v = a23 * a31;
+
+  tsmall[0][0] = 1.0;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 0.0;
+  tsmall[1][1] = 1.0;
+  tsmall[1][2] = 0.0;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = 1.0 - a31;
+  tsmall[2][2] = 0.0;
+  tsmall[2][3] = a31;
+  tsmall[3][0] = 0.0;
+  tsmall[3][1] = 0.0;
+  tsmall[3][2] = a23;
+  tsmall[3][3] = 1.0 - a23;
+}
+
+void libtetrabz_tsmall_c3(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+  /*
+  Cut small tetrahedron C3
+  :param e:
+  :return:
+  */
+  double a23, a13, a30;
+  a23 = (e0 - e[3]) / (e[2] - e[3]);
+  a13 = (e0 - e[3]) / (e[1] - e[3]);
+  a30 = (e0 - e[0]) / (e[3] - e[0]);
+
+  *v = a23 * a13 * a30;
+
+  tsmall[0][0] = 1.0;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 1.0 - a30;
+  tsmall[1][1] = 0.0;
+  tsmall[1][2] = 0.0;
+  tsmall[1][3] = a30;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = a13;
+  tsmall[2][2] = 0.0;
+  tsmall[2][3] = 1.0 - a13;
+  tsmall[3][0] = 0.0;
+  tsmall[3][1] = 0.0;
+  tsmall[3][2] = a23;
+  tsmall[3][3] = 1.0 - a23;
+}
+
+
+void libtetrabz_triangle_a1(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+  /*
+  Cut triangle A1
+  :param e:
+  :return:
+  */
+  double a10, a20, a30;
+  a10 = (e0 - e[0]) / (e[1] - e[0]);
+  a20 = (e0 - e[0]) / (e[2] - e[0]);
+  a30 = (e0 - e[0]) / (e[3] - e[0]);
+
+  *v = 3.0 * a10 * a20 / (e[3] - e[0]);
+
+  tsmall[0][0] = 1.0 - a10;
+  tsmall[0][1] = a10;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 1.0 - a20;
+  tsmall[1][1] = 0.0;
+  tsmall[1][2] = a20;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 1.0 - a30;
+  tsmall[2][1] = 0.0;
+  tsmall[2][2] = 0.0;
+  tsmall[2][3] = a30;
+}
+
+void libtetrabz_triangle_b1(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+) {
+  /*
+  Cut triangle B1
+  :param e:
+  :return:
+  */
+  double a30, a13, a20;
+  a30 = (e0 - e[0]) / (e[3] - e[0]);
+  a13 = (e0 - e[3]) / (e[1] - e[3]);
+  a20 = (e0 - e[0]) / (e[2] - e[0]);
+
+  *v = 3.0 * a30 * a13 / (e[2] - e[0]);
+
+  tsmall[0][0] = 1.0 - a20;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = a20;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 1.0 - a30;
+  tsmall[1][1] = 0.0;
+  tsmall[1][2] = 0.0;
+  tsmall[1][3] = a30;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = a13;
+  tsmall[2][2] = 0.0;
+  tsmall[2][3] = 1.0 - a13;
+}
+
+
+void libtetrabz_triangle_b2(
+    double *e,
+    double e0,
+    double *v,
+    double **tsmall
+)
+{
+  /*
+
+  :param e:
+  :returns:
+      - x - first
+      - y - second
+  */
+  double a12, a31, a20;
+  a12 = (e0 - e[2]) / (e[1] - e[2]);
+  a31 = (e0 - e[1]) / (e[3] - e[1]);
+  a20 = (e0 - e[0]) / (e[2] - e[0]);
+
+  *v = 3.0 * a12 * a31 / (e[2] - e[0]);
+
+  tsmall[0][0] = 1.0 - a20;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = a20;
+  tsmall[0][3] = 0.0;
+  tsmall[1][0] = 0.0;
+  tsmall[1][1] = a12;
+  tsmall[1][2] = 1.0 - a12;
+  tsmall[1][3] = 0.0;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = 1.0 - a31;
+  tsmall[2][2] = 0.0;
+  tsmall[2][3] = a31;
+}
+
+void libtetrabz_triangle_c1(
+  double *e,
+  double e0,
+  double *v,
+  double **tsmall
+)
+{
+/*
+Cut triangle C1
+*/
+  double a03, a13, a23;
+  a03 = (e0 - e[3]) / (e[0] - e[3]);
+  a13 = (e0 - e[3]) / (e[1] - e[3]);
+  a23 = (e0 - e[3]) / (e[2] - e[3]);
+
+  *v = 3.0 * a03 * a13 / (e[3] - e[2]);
+
+  tsmall[0][0] = a03;
+  tsmall[0][1] = 0.0;
+  tsmall[0][2] = 0.0;
+  tsmall[0][3] = 1.0 - a03;
+  tsmall[1][0] = 0.0;
+  tsmall[1][1] = a13;
+  tsmall[1][2] = 0.0;
+  tsmall[1][3] = 1.0 - a13;
+  tsmall[2][0] = 0.0;
+  tsmall[2][1] = 0.0;
+  tsmall[2][2] = a23;
+  tsmall[2][3] = 1.0 - a23;
+}
+
+void eig_sort(
+    int n, //!< [in] the number of components
+    double *key, //!< [in] Variables to be sorted [n].
+    int *swap //!< [out] Order of index (sorted)
+)
+{
+  int i, j, k;
+
+  for (i = 0; i < n; ++i) swap[i] = i;
+
+  for (i = 0; i < n - 1; ++i) {
+    for (j = i + 1; j < n; ++j) {
+      if (key[swap[j]] < key[swap[i]]) {
+        /*
+         Swap
+        */
+        k = swap[j];
+        swap[j] = swap[i];
+        swap[i] = k;
+      }
+    }/*for (j = i + 1; j < n; ++j)*/
+  }/*for (i = 0; i < n - 1; ++i)*/
+}/*eig_sort*/
diff --git a/src/c/libtetrabz_common.h b/src/c/libtetrabz_common.h
new file mode 100644 (file)
index 0000000..e9a8c0f
--- /dev/null
@@ -0,0 +1,36 @@
+/*
+ Copyright (c) 2014 Mitsuaki Kawamura
+
+ Permission is hereby granted, free of charge, to any person obtaining a
+ copy of this software and associated documentation files (the
+ "Software"), to deal in the Software without restriction, including
+ without limitation the rights to use, copy, modify, merge, publish,
+ distribute, sublicense, and/or sell copies of the Software, and to
+ permit persons to whom the Software is furnished to do so, subject to
+ the following conditions:
+
+ The above copyright notice and this permission notice shall be included
+ in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+void libtetrabz_initialize(int *ng, double **bvec, double wlsm[4][20], int **ikv);
+void libtetrabz_tsmall_a1(double *e, double e0, double *v, double tsmall[4][4]);
+void libtetrabz_tsmall_b1(double *e, double e0, double *v, double tsmall[4][4]);
+void libtetrabz_tsmall_b2(double *e, double e0, double *v, double tsmall[4][4]);
+void libtetrabz_tsmall_b3(double *e, double e0, double *v, double tsmall[4][4]);
+void libtetrabz_tsmall_c1(double *e, double e0, double *v, double tsmall[4][4]);
+void libtetrabz_tsmall_c2(double *e, double e0, double *v, double tsmall[4][4]);
+void libtetrabz_tsmall_c3(double *e, double e0, double *v, double tsmall[4][4]);
+void libtetrabz_triangle_a1(double *e, double e0, double *v, double tsmall[3][4]);
+void libtetrabz_triangle_b1(double *e, double e0, double *v, double tsmall[3][4]);
+void libtetrabz_triangle_b2(double *e, double e0, double *v, double tsmall[3][4]);
+void libtetrabz_triangle_c1(double *e, double e0, double *v, double tsmall[3][4]);
+void eig_sort(int n, double *key, int *swap);
diff --git a/src/c/libtetrabz_dbldelta.c b/src/c/libtetrabz_dbldelta.c
new file mode 100644 (file)
index 0000000..fb3d1ba
--- /dev/null
@@ -0,0 +1,249 @@
+//
+// Copyright (c) 2014 Mitsuaki Kawamura
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+// 
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+//
+#include "libtetrabz_common.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+void libtetrabz_dbldelta2(
+    int nb,
+    double **ej,
+    double **w
+) {
+  /*
+  2nd step of tetrahedron method.
+  */
+  int ii, ib, indx[4];
+  double a10, a20, a02, a12, v, e[3];
+
+  for (ib = 0; ib < nb; ib++) {
+
+    for (ii = 0; ii < 3; ii++)
+      w[ib][ii] = 0.0;
+
+    for (ii = 0; ii < 3; ii++) e[ii] = ej[ii][ib];
+    eig_sort(3, e, indx);
+
+    if (e[2] < 1.0e-10) {
+      printf("Nesting ##\n");
+    }
+
+    if ((e[0] < 0.0 && 0.0 <= e[1]) || (e[0] <= 0.0 && 0.0 < e[1])) {
+
+      a10 = (0.0 - e[0]) / (e[1] - e[0]);
+      a20 = (0.0 - e[0]) / (e[2] - e[0]);
+
+      v = a10 / (e[2] - e[0]);
+
+      w[ib][indx[0]] = v * (2.0 - a10 - a20);
+      w[ib][indx[1]] = v * a10;
+      w[ib][indx[2]] = v * a20;
+    }
+    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
+
+      a02 = (0.0 - e[2]) / (e[0] - e[2]);
+      a12 = (0.0 - e[2]) / (e[1] - e[2]);
+
+      v = a12 / (e[2] - e[0]);
+
+      w[ib][indx[0]] = v * a02;
+      w[ib][indx[1]] = v * a12;
+      w[ib][indx[2]] = v * (2.0 - a02 - a12);
+    }
+  }
+}
+
+void dbldelta(
+  int* ng,
+  int nk,
+  int nb,
+  double** bvec,
+  double** eig1,
+  double** eig2,
+  double*** wght
+)
+{
+  /*
+  Main SUBROUTINE for Delta(E1) * Delta(E2)
+  */
+  int it, ik, ib, ii, jj, jb, ** ikv, indx[4];
+  double wlsm[4][20], ** ei1, ** ej1, ** ej2, e[4], *** w1, **w2, v, tsmall[3][4], thr;
+
+  thr = 1.0e-10;
+
+  ikv = (int**)calloc(6 * nk, sizeof(int*));
+  ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
+  for (ik = 0; ik < 6 * nk; ik++) {
+    ikv[ik] = ikv[0] + ik * 20;
+  }
+
+  ei1 = (double**)calloc(4, sizeof(double*));
+  ej1 = (double**)calloc(4, sizeof(double*));
+  ei1[0] = (double*)calloc(4 * nb, sizeof(double));
+  ej1[0] = (double*)calloc(4 * nb, sizeof(double));
+  for (ii = 0; ii < 4; ii++) {
+    ei1[ii] = ei1[0] + ii * nb;
+    ej1[ii] = ej1[0] + ii * nb;
+  }
+
+  w1 = (double***)malloc(nb, sizeof(double**));
+  w1[0] = (double**)malloc(nb * nb, sizeof(double*));
+  w1[0][0] = (double*)malloc(nb * nb * 4, sizeof(double));
+  for (ib = 0; ib < nb; ib++) {
+    w1[ii] = w1[0] + ib * nb;
+    for (jb = 0; jb < nb; jb++) {
+      w1[ib][jb] = w1[0][0] + ib * nb * 4 + jb * 4;
+    }
+  }
+
+  ej2 = (double**)calloc(3, sizeof(double*));
+  ej2[0] = (double*)calloc(3 * nb, sizeof(double));
+  for (ii = 0; ii < 3; ii++) {
+    ej2[ii] = ej2[0] + ii * nb;
+  }
+
+  w2 = (double**)calloc(nb, sizeof(double*));
+  w2[0] = (double*)calloc(nb*3, sizeof(double));
+  for (ib = 0; ib < nb; ib++) {
+    w2[ib] = w2[0] + ib * 3;
+  }
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        wght[ik][ib][jb] = 0.0;
+
+  for (it = 0; it < 6 * nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        ej1[ii][ib] = 0.0;
+        for (jj = 0; jj < 20; jj++) {
+          ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
+          ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
+        }
+      }
+    }
+
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        for (ii = 0; ii < 4; ii++)
+          w1[ib][jb][ii] = 0.0;
+
+    for (ib = 0; ib < nb; ib++) {
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      if (e[0] < 0.0 && e[0] <= e[1]) {
+
+        libtetrabz_triangle_a1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 3; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+          libtetrabz_dbldelta2(nb, ej2, w2);
+            for (ii = 0; ii < 4; ii++)
+              for (jj = 0; jj < 3; jj++)
+                for (jb = 0; jb < nb; jb++)
+                  w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[1] < 0.0 && 0.0 <= e[2]) {
+
+        libtetrabz_triangle_b1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 3; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+          libtetrabz_dbldelta2(nb, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 3; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+        libtetrabz_triangle_b2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 3; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+          libtetrabz_dbldelta2(nb, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 3; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[2] < 0.0 && 0.0 < e[3]) {
+
+        libtetrabz_triangle_c1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 3; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+          libtetrabz_dbldelta2(nb, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 3; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else {
+        continue;
+      }
+    }
+    for (ii = 0; ii < 20; ii++)
+      for (jj = 0; jj < 4; jj++)
+        for (ib = 0; ib < nb; ib++)
+          for (jb = 0; jb < nb; jb++)
+            wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        wght[ik][ib][jb] /= (6.0 * (double)nk);
+
+  free(ikv[0]);
+  free(ikv);
+  free(ei1[0]);
+  free(ei1);
+  free(ej1[0]);
+  free(ej1);
+  free(ej2[0]);
+  free(ej2);
+  free(w1[0][0]);
+  free(w1[0]);
+  free(w1); 
+  free(w2[0]);
+  free(w2);
+}
diff --git a/src/c/libtetrabz_dblstep.c b/src/c/libtetrabz_dblstep.c
new file mode 100644 (file)
index 0000000..f85c0f7
--- /dev/null
@@ -0,0 +1,309 @@
+/*
+ Copyright (c) 2014 Mitsuaki Kawamura
+
+ Permission is hereby granted, free of charge, to any person obtaining a
+ copy of this software and associated documentation files (the
+ "Software"), to deal in the Software without restriction, including
+ without limitation the rights to use, copy, modify, merge, publish,
+ distribute, sublicense, and/or sell copies of the Software, and to
+ permit persons to whom the Software is furnished to do so, subject to
+ the following conditions:
+
+ The above copyright notice and this permission notice shall be included
+ in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+#include "libtetrabz_common.h"
+#include <math.h>
+
+void libtetrabz_dblstep2(
+    int nb,
+    double *ei,
+    double **ej,
+    double **w1
+) {
+  /*
+  Tetrahedron method for theta( - de)
+  */
+  int ii, jj, ib, indx[4];
+  double v, thr, e[4], tsmall[4][4];
+
+  thr = 1.0e-8;
+
+  for (ib = 0; ib < nb; ib++) {
+
+    for (ii = 0; ii < 3; ii++)
+      w1[ib][ii] = 0.0;
+
+    for (ii = 0; ii < 3; ii++) e[ii] = - ei[ii] + ej[ii][ib];
+    eig_sort(4, e, indx);
+
+    if (fabs(e[0]) < thr && fabs(e[3]) < thr) {
+      /*
+      Theta(0) = 0.5
+      */
+      v = 0.5;
+      for (ii = 0; ii < 4; ii++)
+        w1[ib][ii] = v * 0.25;
+    }
+    else if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
+      libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 4; jj++)
+          w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+    }
+    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
+
+      libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 4; jj++)
+          w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+      libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 4; jj++)
+          w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+      libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 4; jj++)
+          w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+    }
+    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
+
+      libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 4; jj++)
+          w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+      libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 4; jj++)
+          w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+      libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 4; jj++)
+          w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+    }
+    else if (e[3] <= 0.0) {
+      for (ii = 0; ii < 4; ii++)
+        w1[ib][ii] = 0.25;
+    }
+  }
+}
+
+void dblstep(
+    int *ng,
+    int nk,
+    int nb,
+    double **bvec,
+    double **eig1,
+    double **eig2,
+    double ***wght
+)
+{
+    /*
+    Main SUBROUTINE for Theta(- E1) * Theta(E1 - E2)
+    */
+  int it, ik, ib, ii, jj, jb, **ikv, indx[4];
+  double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ***w1, **w2, v, tsmall[4][4], thr;
+
+  ikv = (int**)calloc(6 * nk, sizeof(int*));
+  ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
+  for (ik = 0; ik < 6 * nk; ik++) {
+    ikv[ik] = ikv[0] + ik * 20;
+  }
+
+  ei1 = (double**)calloc(4, sizeof(double*));
+  ej1 = (double**)calloc(4, sizeof(double*));
+  ei1[0] = (double*)calloc(4 * nb, sizeof(double));
+  ej1[0] = (double*)calloc(4 * nb, sizeof(double));
+  for (ii = 0; ii < 4; ii++) {
+    ei1[ii] = ei1[0] + ii * nb;
+    ej1[ii] = ej1[0] + ii * nb;
+  }
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        wght[ik][ib][jb] = 0.0;
+
+  thr = 1.0e-8;
+
+  for(it = 0; it < 6*nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        ej1[ii][ib] = 0.0;
+        for (jj = 0; jj < 20; jj++) {
+          ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
+          ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
+        }
+      }
+    }
+
+   
+    for (ib = 0; ib < nb; ib++) {
+
+      for (jb = 0; jb < nb; jb++)
+        for (ii = 0; ii < 4; ii++)
+          w1[ib][jb][ii] = 0.0;
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      if (e[0] <= 0.0 && 0.0 < e[1]) {
+
+        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_dblstep2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[1] <= 0.0 && 0.0 < e[2]) {
+
+        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_dblstep2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_dblstep2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_dblstep2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[2] <= 0.0 && 0.0 < e[3]) {
+
+        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_dblstep2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_dblstep2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_dblstep2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[3] <= 0.0) {
+        for (ii = 0; ii < 4; ii++) {
+          ei2[ii] = ej1[ii][ib];
+          for (jb = 0; jb < nb; jb++)
+            ej2[ii][jb] = ej1[ii][jb];
+        }
+        libtetrabz_dblstep2(nb, ei2, ej2, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jb = 0; jb < nb; jb++)
+            w1[ib][jb][ii] += v * w2[jb][ii];
+      }
+      else {
+      continue;
+      }
+    }
+    for (ii = 0; ii < 20; ii++)
+      for (jj = 0; jj < 4; jj++)
+        for (ib = 0; ib < nb; ib++)
+          for (jb = 0; jb < nb; jb++)
+            wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        wght[ik][ib][jb] /= (6.0 * (double) nk);
+}
diff --git a/src/c/libtetrabz_dos.c b/src/c/libtetrabz_dos.c
new file mode 100644 (file)
index 0000000..ea117f4
--- /dev/null
@@ -0,0 +1,283 @@
+/*
+ Copyright (c) 2022 Mitsuaki Kawamura
+
+ Permission is hereby granted, free of charge, to any person obtaining a
+ copy of this software and associated documentation files (the
+ "Software"), to deal in the Software without restriction, including
+ without limitation the rights to use, copy, modify, merge, publish,
+ distribute, sublicense, and/or sell copies of the Software, and to
+ permit persons to whom the Software is furnished to do so, subject to
+ the following conditions:
+
+ The above copyright notice and this permission notice shall be included
+ in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+#include "libtetrabz_common.h"
+#include <stdlib.h>
+
+void dos(
+  int* ng,
+  int nk,
+  int nb,
+  int ne,
+  double** bvec,
+  double** eig,
+  double* e0,
+  double*** wght) {
+  /*
+  Compute Dos : Delta(E - E1)
+  */
+  int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
+  double wlsm[4][20], ** ei1, e[4], *** w1, v, tsmall[3][4];
+
+  ikv = (int**)calloc(6 * nk, sizeof(int*));
+  ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
+  for (ik = 0; ik < 6 * nk; ik++) {
+    ikv[ik] = ikv[0] + ik * 20;
+  }
+
+  ei1 = (double**)calloc(4, sizeof(double*));
+  ei1[0] = (double*)calloc(4 * nb, sizeof(double));
+  for (ii = 0; ii < 4; ii++) {
+    ei1[ii] = ei1[0] + ii * nb;
+  }
+
+  w1 = (double***)calloc(nb, sizeof(double**));
+  w1[0] = (double**)calloc(nb * ne, sizeof(double*));
+  w1[0][0] = (double*)calloc(nb * ne * 4, sizeof(double));
+  for (ib = 0; ib < nb; ib++) {
+    w1[ii] = w1[0] + ib * ne;
+    for (ie = 0; ie < ne; ie++) {
+      w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
+    }
+  }
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (ie = 0; ie < ne; ie++)
+        wght[ik][ib][ie] = 0.0;
+
+  for (it = 0; it < 6 * nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        for (jj = 0; jj < 20; jj++)
+          ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
+      }
+    }
+
+    for (ib = 0; ib < nb; ib++)
+      for (ie = 0; ie < ne; ie++)
+        for (ii = 0; ii < 4; ii++)
+          w1[ib][ie][ii] = 0.0;
+
+    for (ib = 0; ib < nb; ib++) {
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      for (ie = 0; ie < ne; ie++) {
+
+        if ((e[0] < e0[ie] && e0[ie] <= e[1]) || (e[0] <= e0[ie] && e0[ie] < e[1])) {
+
+          libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 3; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
+
+        }
+        else if ((e[1] < e0[ie] && e0[ie] <= e[2]) || (e[1] <= e0[ie] && e0[ie] < e[2])) {
+
+          libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 3; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
+
+          libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 3; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
+        }
+        else if ((e[2] < e0[ie] && e0[ie] <= e[3]) || (e[2] <= e0[ie] && e0[ie] < e[3])) {
+
+          libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 3; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] / 3.0;
+        }
+        else {
+          continue;
+        }
+        for (ii = 0; ii < 20; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ib = 0; ib < nb; ib++)
+              for (ie = 0; ie < ne; ie++)
+                wght[ikv[it][ii]][ib][ie] += w1[ib][ie][jj] * wlsm[jj][ii];
+
+      }
+    }
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (ie = 0; ie < ne; ie++)
+        wght[ik][ib][ie] /= (6.0 * (double)nk);
+
+  free(ikv[0]);
+  free(ikv);
+  free(ei1[0]);
+  free(ei1);
+  free(w1[0][0]);
+  free(w1[0]);
+  free(w1);
+}
+
+void intdos(
+  int* ng,
+  int nk,
+  int nb,
+  int ne,
+  double** bvec,
+  double** eig,
+  double* e0,
+  double*** wght
+) {
+  /*
+  Compute integrated Dos : theta(E - E1)
+  */
+
+  int it, ik, ib, ii, jj, ie, ** ikv, indx[4];
+  double wlsm[4][20], ** ei1, e[4], *** w1, v, tsmall[4][4];
+
+  ikv = (int**)calloc(6 * nk, sizeof(int*));
+  ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
+  for (ik = 0; ik < 6 * nk; ik++) {
+    ikv[ik] = ikv[0] + ik * 20;
+  }
+
+  ei1 = (double**)calloc(4, sizeof(double*));
+  ei1[0] = (double*)calloc(4 * nb, sizeof(double));
+  for (ii = 0; ii < 4; ii++) {
+    ei1[ii] = ei1[0] + ii * nb;
+  }
+
+  w1 = (double***)calloc(nb, sizeof(double**));
+  w1[0] = (double**)calloc(nb * ne, sizeof(double*));
+  w1[0][0] = (double*)calloc(nb * ne * 4, sizeof(double));
+  for (ib = 0; ib < nb; ib++) {
+    w1[ii] = w1[0] + ib * ne;
+    for (ie = 0; ie < ne; ie++) {
+      w1[ib][ie] = w1[0][0] + ib * ne * 4 + ie * 4;
+    }
+  }
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (ie = 0; ie < ne; ie++)
+        wght[ik][ib][ie] = 0.0;
+
+  for (it = 0; it < 6 * nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        for (jj = 0; jj < 20; jj++)
+          ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
+      }
+    }
+
+    for (ib = 0; ib < nb; ib++)
+      for (ie = 0; ie < ne; ie++)
+        for (ii = 0; ii < 4; ii++)
+          w1[ib][ie][ii] = 0.0;
+
+    for (ib = 0; ib < nb; ib++) {
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      for (ie = 0; ie < ne; ie++) {
+
+        if ((e[0] <= e0[ie] && e0[ie] < e[1]) || (e[0] < e0[ie] && e0[ie] <= e[1])) {
+          libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
+        }
+        else if ((e[1] <= e0[ie] && e0[ie] < e[2]) || (e[1] < e0[ie] && e0[ie] <= e[2])) {
+
+          libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
+
+          libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
+
+          libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
+
+        }
+        else if ((e[2] <= e0[ie] && e0[ie] < e[3]) || (e[2] < e0[ie] && e0[ie] <= e[3])) {
+
+          libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
+
+          libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
+
+          libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              w1[ib][ie][ii] += v * tsmall[jj][ii] * 0.25;
+
+        }
+        else if (e[3] <= e0[ie]) {
+          for (ii = 0; ii < 4; ii++)
+            w1[ib][ie][ii] = 0.25;
+        }
+        else {
+          continue;
+        }
+      }
+    }
+
+    for (ii = 0; ii < 20; ii++)
+      for (jj = 0; jj < 4; jj++)
+        for (ib = 0; ib < nb; ib++)
+          for (ie = 0; ie < ne; ie++)
+            wght[ikv[it][ii]][ib][ie] += w1[ib][ie][jj] * wlsm[jj][ii];
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (ie = 0; ie < ne; ie++)
+        wght[ik][ib][ie] /= (6.0 * (double)nk);
+
+  free(ikv[0]);
+  free(ikv);
+  free(ei1[0]);
+  free(ei1);
+  free(w1[0][0]);
+  free(w1[0]);
+  free(w1);
+}
diff --git a/src/c/libtetrabz_fermigr.c b/src/c/libtetrabz_fermigr.c
new file mode 100644 (file)
index 0000000..8460d5d
--- /dev/null
@@ -0,0 +1,443 @@
+//
+// Copyright (c) 2014 Mitsuaki Kawamura
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+// 
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+//
+#include "libtetrabz_common.h"
+
+void libtetrabz_fermigr3(
+    int ne,
+    double *e0,
+    double de[4],
+    double **w1
+) {
+  int ii, jj, ie, indx[3];
+  double e[4], tsmall[3][4], v;
+
+  for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
+  eig_sort(3, e, indx);
+
+  for (ie = 0; ie < ne; ie++) {
+
+    for (ii = 0; ii < 4; ii++) w1[ie][ii] = 0.0;
+
+    if (e[0] < e0[ie] && e0[ie] <= e[1]) {
+
+      libtetrabz_triangle_a1(e, e0[ie], &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 3; jj++)
+          w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
+    }
+    else if (e[1] < e0[ie] && e0[ie] <= e[2]) {
+
+      libtetrabz_triangle_b1(e, e0[ie], &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 3; jj++)
+          w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
+
+      libtetrabz_triangle_b2(e, e0[ie], &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 3; jj++)
+          w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
+    }
+    else if (e[2] < e0[ie] && e0[ie] < e[3]) {
+
+      libtetrabz_triangle_c1(e, e0[ie], &v, tsmall);
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 3; jj++)
+          w1[ie][ii] += v * tsmall[jj][ii] / 3.0;
+    }
+  }
+}
+
+void libtetrabz_fermigr2(
+    int nb,
+    int ne,
+    double *e0,
+    double *ei1,
+    double **ej1,
+    double ***w1
+) {
+  int ib, ii, jj, ie, indx[4];
+  double e[4], tsmall[4][4], v, de[4], thr, **w2;
+
+  for (ib = 0; ib < nb; ib++) {
+
+    for (ie = 0; ie < ne; ie++)
+      for (ii = 0; ii < 4; ii++)
+        w1[ib][ie][ii] = 0.0;
+
+    for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
+    eig_sort(4, e, indx);
+
+    if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
+
+      libtetrabz_tsmall_a1(e, e0[ie], &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_fermigr3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
+      }
+    }
+    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
+
+      libtetrabz_tsmall_b1(e, e0[ie], &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_fermigr3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_b2(e, e0[ie], &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_fermigr3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_b3(e, e0[ie], &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_fermigr3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
+      }
+    }
+    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
+
+      libtetrabz_tsmall_c1(e, e0[ie], &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_fermigr3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_c2(e, e0[ie], &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_fermigr3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_c3(e, e0[ie], &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_fermigr3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][ie][indx[ii]] += v * w2[ie][jj] * tsmall[jj][ii];
+      }
+    }
+    else if (e[3] <= 0.0) {
+      for (ii = 0; ii < 4; ii++)
+        de[ii] = ej1[ii][ib] - ei1[ii];
+      libtetrabz_fermigr3(ne, e0, de, w2);
+      for (ii = 0; ii < 4; ii++)
+        for (ie = 0; ie < ne; ie++)
+          w1[ib][ie][ii] += w2[ie][ii];
+    }
+  }
+}
+
+void fermigr(
+    int *ng,
+    int nk,
+    int nb,
+    int ne,
+    double **bvec,
+    double **eig1,
+    double **eig2,
+    double *e0,
+    double ****wght
+)
+{
+  /*
+  Main SUBROUTINE for Fermi's Golden rule : Theta(- E1) * Theta(E2) * Delta(E2 - E1 - w)
+  */
+  int it, ik, ib, ii, jj, jb, ie, **ikv, indx[4];
+  double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ****w1, ***w2, v, tsmall[4][4], thr;
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        for (ie = 0; ie < ne; ie++)
+          wght[ik][ib][jb][ie] = 0.0;
+
+  thr = 1.0e-10;
+
+  for(it = 0; it < 6*nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        ej1[ii][ib] = 0.0;
+      }
+    }
+    for (jj = 0; jj < 20; jj++) {
+      for (ib = 0; ib < nb; ib++) {
+        for (ii = 0; ii < 4; ii++) {
+          ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
+          ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
+        }
+      }
+    }
+
+    for (ib = 0; ib < nb; ib++) {
+
+      for (jb = 0; jb < nb; jb++)
+        for (ii = 0; ii < 4; ii++)
+          for (ie = 0; ie < ne; ie++)
+            w1[ib][jb][ie][ii] = 0.0;
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      if (e[0] <= 0.0 && 0.0 < e[1]) {
+
+        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[1] <= 0.0 && 0.0 < e[2]) {
+
+        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[2] <= 0.0 && 0.0 < e[3]) {
+
+        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  w1[ib][jb][ie][indx[ii]] += v * w2[jb][ie][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[3] <= 0.0) {
+
+        for (ii = 0; ii < 4; ii++) {
+          ei2[ii] = ej1[ii][ib];
+          for (jb = 0; jb < nb; jb++)
+            ej2[ii][jb] = ej1[ii][jb];
+        }
+        libtetrabz_fermigr2(nb, ne, e0, ei2, ej2, w2);
+
+        for (ii = 0; ii < 4; ii++)
+          for (jb = 0; jb < nb; jb++)
+            for (ie = 0; ie < ne; ie++)
+              w1[ib][jb][ie][ii] += w2[jb][ie][ii];
+      }
+      else {
+        continue;
+      }
+    }
+    for (ii = 0; ii < 20; ii++)
+      for (jj = 0; jj < 4; jj++)
+        for (ib = 0; ib < nb; ib++)
+          for (jb = 0; jb < nb; jb++)
+            for (ie = 0; ie < ne; ie++)
+              wght[ikv[it][ii]][ib][jb][ie] += w1[ib][jb][ie][jj] * wlsm[jj][ii];
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        for (ie = 0; ie < ne; ie++)
+          wght[ik][ib][jb][ie] /= (6.0 * (double) nk);
+}
diff --git a/src/c/libtetrabz_occ.c b/src/c/libtetrabz_occ.c
new file mode 100644 (file)
index 0000000..7df5fe3
--- /dev/null
@@ -0,0 +1,207 @@
+//
+// Copyright (c) 2014 Mitsuaki Kawamura
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+// 
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+//
+#include <math.h>
+#include <stdlib.h>
+#include "libtetrabz_common.h"
+
+void occ(
+  int* ng,
+  int nk,
+  int nb,
+  double** bvec,
+  double** eig,
+  double** wght
+) {
+  /*
+  Main SUBROUTINE for occupation : Theta(EF - E1)
+  */
+  int it, ik, ib, ii, jj, ** ikv, indx[4];
+  double wlsm[4][20], ** ei1, e[4], ** w1, v, tsmall[4][4];
+
+  ikv = (int**)calloc(6 * nk, sizeof(int*));
+  ikv[0] = (int*)calloc(6 * nk * 20, sizeof(int));
+  for (ik = 0; ik < 6 * nk; ik++) {
+    ikv[ik] = ikv[0] + ik * 20;
+  }
+
+  ei1 = (double**)calloc(4, sizeof(double*));
+  ei1[0] = (double*)calloc(4 * nb, sizeof(double));
+  for (ii = 0; ii < 4; ii++) {
+    ei1[ii] = ei1[0] + ii * nb;
+  }
+
+  w1 = (double**)calloc(nb, sizeof(double*));
+  w1[0] = (double*)calloc(nb * 4, sizeof(double));
+  for (ib = 0; ib < nb; ib++) {
+    w1[ii] = w1[0] + ib * 4;
+  }
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      wght[ik][ib] = 0.0;
+
+  for (it = 0; it < 6 * nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        for (jj = 0; jj < 20; jj++)
+          ei1[ii][ib] += wlsm[ii][jj] * eig[ikv[it][jj]][ib];
+      }
+    }
+
+    for (ib = 0; ib < nb; ib++)
+      for (ii = 0; ii < 4; ii++)
+        w1[ib][ii] = 0.0;
+
+    for (ib = 0; ib < nb; ib++) {
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      if (e[0] <= 0.0 && 0.0 < e[1]) {
+        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+      }
+      else if (e[1] <= 0.0 && 0.0 < e[2]) {
+
+        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+      }
+      else if (e[2] <= 0.0 && 0.0 < e[3]) {
+
+        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+
+        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][ii] += v * tsmall[jj][ii] * 0.25;
+      }
+      else if (e[3] <= 0.0) {
+        for (ii = 0; ii < 4; ii++)
+          w1[ib][ii] = 0.25;
+      }
+      else {
+        continue;
+      }
+    }
+    for (ii = 0; ii < 20; ii++)
+      for (jj = 0; jj < 4; jj++)
+        for (ib = 0; ib < nb; ib++)
+          wght[ikv[it][ii]][ib] += w1[ib][jj] * wlsm[jj][ii];
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      wght[ik][ib] /= (6.0 * (double)nk);
+
+  free(ikv[0]);
+  free(ikv);
+  free(ei1[0]);
+  free(ei1);
+  free(w1[0]);
+  free(w1);
+}
+
+void fermieng(
+    int *ng,
+    int nk,
+    int nb,
+    double **bvec,
+    double **eig,
+    double **wght,
+  double nelec,
+  int *iteration,
+  double *ef
+) {
+  /*
+  Calculate Fermi energy
+  */
+  int maxiter, ik, ib;
+  double eps, elw, eup, **eig2, sumkmid;
+
+  maxiter = 300;
+  eps = 1.0e-10;
+
+  elw = eig[0][0];
+  eup = eig[0][0];
+  for (ik = 0; ik < nk; ik++) {
+    for (ib = 0; ib < nb; ib++) {
+      if (elw > eig[ik][ib]) elw = eig[ik][ib];
+      if (eup < eig[ik][ib]) eup = eig[ik][ib];
+    }
+  }
+  /*
+  Bisection method
+  */
+  for (*iteration = 0; *iteration < maxiter; *iteration++) {
+
+    *ef = (eup + elw) * 0.5;
+    /*
+    Calc. # of electrons
+    */
+    for (ik = 0; ik < nk; ik++)
+      for (ib = 0; ib < nb; ib++)
+        eig2[ik][ib] = eig[ik][ib] - *ef;
+    occ(ng, nk, nb, bvec, eig2, wght);
+
+    sumkmid = 0.0;
+    for (ik = 0; ik < nk; ik++) {
+      for (ib = 0; ib < nb; ib++) {
+        sumkmid += wght[ik][ib];
+      }
+    }
+    /*
+    convergence check
+    */
+    if (fabs(sumkmid - nelec) < eps)
+      return;
+    else if (sumkmid < nelec)
+      elw = *ef;
+    else
+      eup = *ef;
+  }
+}
diff --git a/src/c/libtetrabz_polcmplx.c b/src/c/libtetrabz_polcmplx.c
new file mode 100644 (file)
index 0000000..0ad8053
--- /dev/null
@@ -0,0 +1,823 @@
+/*
+ Copyright (c) 2014 Mitsuaki Kawamura
+
+ Permission is hereby granted, free of charge, to any person obtaining a
+ copy of this software and associated documentation files (the
+ "Software"), to deal in the Software without restriction, including
+ without limitation the rights to use, copy, modify, merge, publish,
+ distribute, sublicense, and/or sell copies of the Software, and to
+ permit persons to whom the Software is furnished to do so, subject to
+ the following conditions:
+
+ The above copyright notice and this permission notice shall be included
+ in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+#include "libtetrabz_common.h"
+#include <math.h>
+
+void polcmplx(
+    int *ng,
+    int nk,
+    int nb,
+    int ne,
+    double **bvec,
+    double **eig1,
+    double **eig2,
+    double **e0,
+    double *****wght) {
+  /*
+  Main SUBROUTINE for Polarization (Imaginary axis) : Theta(- E1) * Theta(E2) / (E2 - E1 - iw)
+  */
+  int it, ik, ie, ib, ii, jj, kk, jb, **ikv, indx[4];
+  double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], *****w1, ****w2, v, tsmall[4][4], thr;
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        for (ie = 0; ie < ne; ie++)
+          for (jj = 0; jj < 2; jj++)
+            wght[ik][ib][jb][ie][jj] = 0.0;
+
+  thr = 1.0e-8;
+
+  for (it = 0; it < 6 * nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        ej1[ii][ib] = 0.0;
+      }
+    }
+    for (jj = 0; jj < 20; jj++) {
+      for (ib = 0; ib < nb; ib++) {
+        for (ii = 0; ii < 4; ii++) {
+          ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
+          ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
+        }
+      }
+    }
+
+    for (ib = 0; ib < nb; ib++) {
+
+      for (jb = 0; jb < nb; jb++)
+        for (ii = 0; ii < 4; ii++)
+          for (ie = 0; ie < ne; ie++)
+            for (jj = 0; jj < 2; jj++)
+              w1[ib][jb][ie][ii][jj] = 0.0;
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      if (e[0] <= 0.0 && 0.0 < e[1]) {
+
+        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  for (kk = 0; kk < 4; kk++)
+                    w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[1] <= 0.0 && 0.0 < e[2]) {
+
+        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  for (kk = 0; kk < 4; kk++)
+                    w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  for (kk = 0; kk < 4; kk++)
+                    w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  for (kk = 0; kk < 4; kk++)
+                    w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[2] <= 0.0 && 0.0 < e[3]) {
+
+        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  for (kk = 0; kk < 4; kk++)
+                    w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  for (kk = 0; kk < 4; kk++)
+                    w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++) {
+            ei2[ii] = 0.0;
+            for (jb = 0; jb < nb; jb++) ej2[ii][jb] = 0.0;
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] += tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] += tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          }
+          libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                for (ie = 0; ie < ne; ie++)
+                  for (kk = 0; kk < 4; kk++)
+                    w1[ib][jb][ie][kk][indx[ii]] += v * w2[jb][ie][kk][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[3] <= 0.0) {
+
+        for (ii = 0; ii < 4; ii++) {
+          ei2[ii] = ej1[ii][ib];
+          for (jb = 0; jb < nb; jb++) ej2[ii][jb] = ej1[ii][jb];
+        }
+        libtetrabz_polcmplx2(nb, ne, e0, ei2, ej2, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jb = 0; jb < nb; jb++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 4; kk++)
+                w1[ib][jb][ie][kk][ii] += w2[jb][ie][kk][ii];
+      }
+      else {
+        continue;
+      }
+    }
+    for (ii = 0; ii < 20; ii++)
+      for (jj = 0; jj < 4; jj++)
+        for (ib = 0; ib < nb; ib++)
+          for (jb = 0; jb < nb; jb++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                wght[ikv[it][ii]][ib][jb][ie][kk] += w1[ib][jb][ie][kk][jj] * wlsm[jj][ii];
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        for (ie = 0; ie < ne; ie++)
+          for (ii = 0; ii < 2; ii++)
+            wght[ik][ib][jb][ie][ii] /= (6.0 * (double) nk);
+}
+
+void libtetrabz_polcmplx2(
+    int nb,
+    int ne,
+    double *e0,
+    double *ei1,
+    double **ej1,
+    double ****w1
+) {
+  /*
+  Tetrahedron method for theta( - E2)
+  */
+  int ib, ii, jj, kk, ie, indx[4];
+  double e[4], tsmall[4][4], v, de[4], thr, ***w2;
+
+  thr = 1.0e-8;
+
+  for (ib = 0; ib < nb; ib++) {
+
+    for (ie = 0; ie < ne; ie++)
+      for (ii = 0; ii < 4; ii++)
+        for (jj = 0; jj < 2; jj++)
+          w1[ib][ie][jj][ii] = 0.0;
+
+    for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
+    eig_sort(4, e, indx);
+
+    if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
+
+      libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polcmplx3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
+      }
+    }
+    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
+
+      libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polcmplx3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polcmplx3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polcmplx3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
+      }
+    }
+    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
+
+      libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polcmplx3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polcmplx3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polcmplx3(ne, e0, de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            for (ie = 0; ie < ne; ie++)
+              for (kk = 0; kk < 2; kk++)
+                w1[ib][ie][kk][indx[ii]] += v * w2[jj][ie][kk] * tsmall[jj][ii];
+      }
+    }
+    else if (e[3] <= 0.0) {
+      for (ii = 0; ii < 4; ii++) {
+        de[ii] = tsmall[ii][jj] * (ej1[ii][ib] - ei1[ii]);
+      }
+      libtetrabz_polcmplx3(ne, e0, de, w2);
+      for (ii = 0; ii < 4; ii++)
+        for (ie = 0; ie < ne; ie++)
+          for (kk = 0; kk < 2; kk++)
+            w1[ib][ie][kk][ii] += w2[ie][kk][ii];
+    }
+  }
+}
+
+
+void libtetrabz_polcmplx3(
+    int ne,
+    double **e0,
+    double de[4],
+    double ***w1
+) {
+  /*
+  Tetrahedron method for delta(om - ep + e)
+  */
+  int ii, jj, indx[4], ie;
+  double e[4], x[4], thr, w2[4][2];
+
+  for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
+  eig_sort(4, e, indx);
+
+  for (ie = 0; ie < ne; ie++) {
+    /*
+    I am not sure which one is better.
+    The former is more stable.
+    The latter is more accurate ?
+    */
+    for (ii = 0; ii < 4; ii++)
+      for (jj = 0; jj < 2; jj++)
+        w1[ie][jj][ii] = 0.25 / (de[ii] + e0[ie][jj]);
+
+    continue;
+
+    for (ii = 0; ii < 4; ii++)
+      x[ii] = (e[ii] + e0[ie][0]) / e0[ie][1];
+    /* thr = maxval(de(1:4)) * 1d-3*/
+    thr = fabs(x[0]);
+    for(ii=0;ii<4;ii++)
+      if(thr < fabs(x[ii])) thr = fabs(x[ii]);
+    thr = max(1.0e-3, thr * 1.0e-2);
+
+    if (fabs(x[3] - x[2]) < thr) {
+      if (fabs(x[3] - x[1]) < thr) {
+        if (fabs(x[3] - x[0]) < thr) {
+          /*
+          e[3] = e[2] = e[1] = e[0]
+          */
+          w2[3][0] = 0.25 * x[3] / (1.0 + x[3] *x[3]);
+          w2[3][1] = 0.25 / (1.0 + x[3] * x[3]);
+          for (ii = 0; ii < 2; ii++) {
+            w2[2][ii] = w2[3][ii];
+            w2[1][ii] = w2[3][ii];
+            w2[0][ii] = w2[3][ii];
+          }
+        }
+        else {
+          /*
+          e[3] = e[2] = e[1]
+          */
+          libtetrabz_polcmplx_1211(x[3], x[0], w2[3]);
+          for (ii = 0; ii < 2; ii++) {
+            w2[2][ii] = w2[3][ii];
+            w2[1][ii] = w2[3][ii];
+          }
+          libtetrabz_polcmplx_1222(x[0], x[3], w2[0]);
+          /*
+          # IF(ANY(w2(1:2,1:4) < 0.0)):
+          #   WRITE(*,*) ie
+          #   WRITE(*,'(100e15.5)') x[0:4]
+          #   WRITE(*,'(2e15.5)') w2(1:2,1:4)
+          #   STOP "weighting 4=3=2"
+          */
+        }
+      }
+      else if (fabs(x[1] - x[0]) < thr) {
+        /*
+        e[3] = e[2], e[1] = e[0]
+        */
+        libtetrabz_polcmplx_1221(x[3], x[1], w2[3]);
+        for (ii = 0; ii < 2; ii++) w2[2][ii] = w2[3][ii];
+        libtetrabz_polcmplx_1221(x[1], x[3], w2[1]);
+        for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[1][ii];
+        /*
+         IF(ANY(w2(1:2,1:4) < 0.0)):
+           WRITE(*,*) ie
+           WRITE(*,'(100e15.5)') x[0:4]
+           WRITE(*,'(2e15.5)') w2(1:2,1:4)
+           STOP "weighting 4=3 2=1"
+        */
+      }
+      else {
+        /*
+        e[3] = e[2]
+        */
+        libtetrabz_polcmplx_1231(x[3], x[0], x[1], w2[3]);
+        for (ii = 0; ii < 2; ii++) w2[2][ii] = w2[3][ii];
+        libtetrabz_polcmplx_1233(x[1], x[0], x[3], w2[1]);
+        libtetrabz_polcmplx_1233(x[0], x[1], x[3], w2[0]);
+        /*
+         IF(ANY(w2(1:2,1:4) < 0.0)):
+           WRITE(*,*) ie
+           WRITE(*,'(100e15.5)') x[0:4]
+           WRITE(*,'(2e15.5)') w2(1:2,1:4)
+           STOP "weighting 4=3"
+        */
+      }
+    }
+    else if (fabs(x[2] - x[1]) < thr) {
+      if (fabs(x[2] - x[0]) < thr) {
+        /*
+        e[2] = e[1] = e[0]
+        */
+        libtetrabz_polcmplx_1222(x[3], x[2], w2[3]);
+        libtetrabz_polcmplx_1211(x[2], x[3], w2[2]);
+        for (ii = 0; ii < 2; ii++) w2[1][ii] = w2[2][ii];
+        for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[2][ii];
+        /*
+        IF(ANY(w2(1:2,1:4) < 0.0)):
+          WRITE(*,*) ie
+          WRITE(*,'(100e15.5)') x[0:4]
+          WRITE(*,'(2e15.5)') w2(1:2,1:4)
+          STOP "weighting 3=2=1"
+        */
+      }
+      else {
+        /*
+        e[2] = e[1]
+        */
+        libtetrabz_polcmplx_1233(x[3], x[0], x[2], w2[3]);
+        libtetrabz_polcmplx_1231(x[2], x[0], x[3], w2[2]);
+        for (ii = 0; ii < 2; ii++) w2[1][ii] = w2[2][ii];
+        libtetrabz_polcmplx_1233(x[0], x[3], x[2], w2[0]);
+        /*
+        IF(ANY(w2(1:2,1:4) < 0.0)):
+          WRITE(*,*) ie
+          WRITE(*,'(100e15.5)') x[0:4]
+          WRITE(*,'(2e15.5)') w2(1:2,1:4)
+          STOP "weighting 3=2"
+        */
+      }
+    }
+    else if (fabs(x[1] - x[0]) < thr) {
+      /*
+      e[1] = e[0]
+      */
+      libtetrabz_polcmplx_1233(x[3], x[2], x[1], w2[3]);
+      libtetrabz_polcmplx_1233(x[2], x[3], x[1], w2[2]);
+      libtetrabz_polcmplx_1231(x[1], x[2], x[3], w2[1]);
+      for (ii = 0; ii < 2; ii++) w2[0][ii] = w2[1][ii];
+      /*
+      IF(ANY(w2(1:2,1:4) < 0.0)):
+        WRITE(*,*) ie
+        WRITE(*,'(100e15.5)') x[0:4]
+        WRITE(*,'(2e15.5)') w2(1:2,1:4)
+        STOP "weighting 2=1"
+      */
+    }
+    else {
+      /*
+      Different each other.
+      */
+      libtetrabz_polcmplx_1234(x[3], x[0], x[1], x[2], w2[3]);
+      libtetrabz_polcmplx_1234(x[2], x[0], x[1], x[3], w2[2]);
+      libtetrabz_polcmplx_1234(x[1], x[0], x[2], x[3], w2[1]);
+      libtetrabz_polcmplx_1234(x[0], x[1], x[2], x[3], w2[0]);
+      /*
+      IF(ANY(w2(1:2,1:4) < 0.0)):
+        WRITE(*,*) ie
+        WRITE(*,'(100e15.5)') x[0:4]
+        WRITE(*,'(2e15.5)') w2(1:2,1:4)
+        STOP "weighting"
+      */
+    }
+    for (ii = 0; ii < 4; ii++) {
+      w1[ie][indx[ii]][0] = w2[ii][0] / e0[ie][1];
+      w1[ie][indx[ii]][1] = w2[ii][1] / (-e0[ie][1]);
+    }
+  }
+}
+//
+// Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
+//  for 0<x<1, 0<y<1-x, 0<z<1-x-y
+
+
+void libtetrabz_polcmplx_1234(
+    double g1,
+    double g2,
+    double g3,
+    double g4,
+    double *w
+    ) {
+  /*
+  1, Different each other
+  */
+  double w2, w3, w4;
+  /*
+  Real
+  */
+  w2 = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1)) +
+       (g2 * g2 - 3.0) * g2 * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w2 = -2.0 * (g2 * g2 - 1.0) + w2 / (g2 - g1);
+  w2 = w2 / (g2 - g1);
+  w3 = 2.0 * (3.0 * g3 * g3 - 1.0) * (atan(g3) - atan(g1)) +
+       (g3 * g3 - 3.0) * g3 * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = -2.0 * (g3 * g3 - 1.0) + w3 / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w4 = 2.0 * (3.0 * g4 * g4 - 1.0) * (atan(g4) - atan(g1)) +
+       (g4 * g4 - 3.0) * g4 * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
+  w4 = -2.0 * (g4 * g4 - 1.0) + w4 / (g4 - g1);
+  w4 = w4 / (g4 - g1);
+  w2 = (w2 - w3) / (g2 - g3);
+  w4 = (w4 - w3) / (g4 - g3);
+  w[0] = (w4 - w2) / (2.0 * (g4 - g2));
+  /*
+  Imaginal
+  */
+  w2 = 2.0 * (3.0 - g2 * g2) * g2 * (atan(g2) - atan(g1)) +
+       (3.0 * g2 * g2 - 1.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w2 = 4.0 * g2 - w2 / (g2 - g1);
+  w2 = w2 / (g2 - g1);
+  w3 = 2.0 * (3.0 - g3 * g3) * g3 * (atan(g3) - atan(g1)) +
+       (3.0 * g3 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = 4.0 * g3 - w3 / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w4 = 2.0 * (3.0 - g4 * g4) * g4 * (atan(g4) - atan(g1)) +
+       (3.0 * g4 * g4 - 1.0) * log((1.0 + g4 * g4) / (1.0 + g1 * g1));
+  w4 = 4.0 * g4 - w4 / (g4 - g1);
+  w4 = w4 / (g4 - g1);
+  w2 = (w2 - w3) / (g2 - g3);
+  w4 = (w4 - w3) / (g4 - g3);
+  w[1] = (w4 - w2) / (2.0 * (g4 - g2));
+}
+
+
+void libtetrabz_polcmplx_1231(
+    double g1,
+    double g2,
+    double g3,
+  double *w
+    ) {
+  /*
+  2, g4 = g1
+  */
+  double w2, w3;
+  /*
+  Real
+  */
+  w2 = 2.0 * (-1.0 + 3.0 * g2 * g2) * (atan(g2) - atan(g1))
+       + g2 * (-3.0 + g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w2 = 2.0 * (1.0 - g2 * g2) + w2 / (g2 - g1);
+  w2 = -g1 + w2 / (g2 - g1);
+  w2 = w2 / (g2 - g1);
+  w3 = 2.0 * (-1.0 + 3.0 * g3 * g3) * (atan(g3) - atan(g1))
+       + g3 * (-3.0 + g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = 2.0 * (1 - g3 * g3) + w3 / (g3 - g1);
+  w3 = -g1 + w3 / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w[0] = (w3 - w2) / (2.0 * (g3 - g2));
+  /*
+  Imaginal
+  */
+  w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
+       (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w2 = 4.0 * g2 - w2 / (g2 - g1);
+  w2 = 1 + w2 / (g2 - g1);
+  w2 = w2 / (g2 - g1);
+  w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
+       (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = 4.0 * g3 - w3 / (g3 - g1);
+  w3 = 1 + w3 / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w[1] = (w3 - w2) / (2.0 * (g3 - g2));
+}
+
+
+void libtetrabz_polcmplx_1233(
+    double g1,
+    double g2,
+    double g3,
+    double *w
+    ) {
+  /*
+  3, g4 = g3
+  */
+  double w2, w3;
+  /*
+  Real
+  */
+  w2 = 2.0 * (1.0 - 3.0 * g2 * g2) * (atan(g2) - atan(g1)) +
+       g2 * (3.0 - g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w2 = 2.0 * (1 - g2 * g2) - w2 / (g2 - g1);
+  w2 = w2 / (g2 - g1);
+  w3 = 2.0 * (1.0 - 3.0 * g3 * g3) * (atan(g3) - atan(g1)) +
+       g3 * (3.0 - g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = 2.0 * (1 - g3 * g3) - w3 / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w2 = (w3 - w2) / (g3 - g2);
+  w3 = 4.0 * (1.0 - 3.0 * g1 * g3) * (atan(g3) - atan(g1))
+       + (3.0 * g1 + 3.0 * g3 - 3.0 * g1 * g3 * g3 + g3 * g3*g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = -4.0 * (1.0 - g1 * g1) + w3 / (g3 - g1);
+  w3 = 4.0 * g1 + w3 / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w[0] = (w3 - w2) / (2.0 * (g3 - g2));
+  /*
+  Imaginal
+  */
+  w2 = 2.0 * g2 * (3.0 - g2 * g2) * (atan(g2) - atan(g1)) +
+       (-1.0 + 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w2 = 4.0 * g2 - w2 / (g2 - g1);
+  w2 = w2 / (g2 - g1);
+  w3 = 2.0 * g3 * (3.0 - g3 * g3) * (atan(g3) - atan(g1)) +
+       (-1.0 + 3.0 * g3 * g3) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = 4.0 * g3 - w3 / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w2 = (w3 - w2) / (g3 - g2);
+  w3 = (3.0 * g1 - 3.0 * g1 * g3 * g3 + 3.0 * g3 + g3 * g3*g3) * (atan(g3) - atan(g1))
+       + (3.0 * g1 * g3 - 1.0) * log((1.0 + g3 * g3) / (1.0 + g1 * g1));
+  w3 = w3 / (g3 - g1) - 4.0 * g1;
+  w3 = w3 / (g3 - g1) - 2.0;
+  w3 = (2.0 * w3) / (g3 - g1);
+  w[1] = (w3 - w2) / (2.0 * (g3 - g2));
+}
+
+void libtetrabz_polcmplx_1221(
+    double g1,
+    double g2,
+    double *w
+    ) {
+  /*
+  4, g4 = g1 and g3 = g2
+  */
+  /*
+  Real
+  */
+  w[0] = -2.0 * (-1.0 + 2.0 * g1 * g2 + g2 * g2) * (atan(g2) - atan(g1))
+         + (g1 + 2.0 * g2 - g1 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w[0] = 2.0 * (-1.0 + g1 * g1) + w[0] / (g2 - g1);
+  w[0] = 3.0 * g1 + w[0] / (g2 - g1);
+  w[0] = 2.0 + (3.0 * w[0]) / (g2 - g1);
+  w[0] = w[0] / (2.0 * (g2 - g1));
+  /*
+  Imaginal
+  */
+  w[1] = 2.0 * (g1 + 2.0 * g2 - g1 * g2 * g2) * (atan(g2) - atan(g1))
+         + (-1.0 + 2.0 * g1 * g2 + g2 * g2) * log((1 + g2 * g2) / (1 + g1 * g1));
+  w[1] = -4.0 * g1 + w[1] / (g2 - g1);
+  w[1] = -3.0 + w[1] / (g2 - g1);
+  w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
+}
+
+
+void libtetrabz_polcmplx_1222(
+    double g1,
+    double g2,
+    double *w
+    ) {
+  /*
+  5, g4 = g3 = g2
+  */
+  /*
+  Real
+  */
+  w[0] = 2.0 * (-1.0 + g1 * g1 + 2.0 * g1 * g2) * (atan(g2) - atan(g1))
+         + (-2.0 * g1 - g2 + g1 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
+  w[0] = g1 - w[0] / (g2 - g1);
+  w[0] = 1.0 - (3.0 * w[0]) / (g2 - g1);
+  w[0] = w[0] / (2.0 * (g2 - g1));
+  /*
+  Imaginal
+  */
+  w[1] = 2.0 * (-2.0 * g1 - g2 + g1 * g1 * g2) * (atan(g2) - atan(g1))
+         + (1.0 - g1 * g1 - 2.0 * g1 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w[1] = 4.0 * g1 + w[1] / (g2 - g1);
+  w[1] = 1.0 + w[1] / (g2 - g1);
+  w[1] = (3.0 * w[1]) / (2.0 * (g2 - g1) * (g2 - g1));
+}
+
+
+void libtetrabz_polcmplx_1211(
+    double g1,
+    double g2,
+    double *w
+) {
+  /*
+  6, g4 = g3 = g1
+  */
+  /*
+  Real
+  */
+  w[0] = 2.0 * (3.0 * g2 * g2 - 1.0) * (atan(g2) - atan(g1))
+         + g2 * (g2 * g2 - 3.0) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w[0] = 2.0 * (1.0 - g1 * g1) + w[0] / (g2 - g1);
+  w[0] = -5.0 * g1 + w[0] / (g2 - g1);
+  w[0] = -11.0 + (3.0 * w[0]) / (g2 - g1);
+  w[0] = w[0] / (6.0 * (g2 - g1));
+  /*
+  Imaginal
+  */
+  w[1] = 2.0 * g2 * (-3.0 + g2 * g2) * (atan(g2) - atan(g1))
+         + (1.0 - 3.0 * g2 * g2) * log((1.0 + g2 * g2) / (1.0 + g1 * g1));
+  w[1] = 4.0 * g2 + w[1] / (g2 - g1);
+  w[1] = 1.0 + w[1] / (g2 - g1);
+  w[1] = w[1] / (2.0 * (g2 - g1) * (g2 - g1));
+}
diff --git a/src/c/libtetrabz_polstat.c b/src/c/libtetrabz_polstat.c
new file mode 100644 (file)
index 0000000..a2c9786
--- /dev/null
@@ -0,0 +1,627 @@
+//
+// Copyright (c) 2014 Mitsuaki Kawamura
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+// 
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+//
+#include "libtetrabz_common.h"
+#include <math.h>
+#include <stdio.h>
+
+/*
+ Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
+ for 0<x<1, 0<y<1-x, 0<z<1-x-y
+*/
+double libtetrabz_polstat_1234(
+    double g1,
+    double g2,
+    double g3,
+    double g4,
+    double lng1,
+    double lng2,
+    double lng3,
+    double lng4)
+{
+  /*
+  1, Different each other
+  */
+  double w2, w3, w4, w;
+  w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 / (g2 - g1);
+  w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 / (g3 - g1);
+  w4 = ((lng4 - lng1) / (g4 - g1) * g4 - 1.0) * g4 / (g4 - g1);
+  w2 = ((w2 - w3) * g2) / (g2 - g3);
+  w4 = ((w4 - w3) * g4) / (g4 - g3);
+  w = (w4 - w2) / (g4 - g2);
+  return w;
+}
+
+
+double libtetrabz_polstat_1231(
+    double g1,
+    double g2,
+    double g3,
+    double lng1,
+    double lng2,
+    double lng3
+) {
+  /*
+  2, g4 = g1
+  */
+  double w2, w3, w;
+  w2 = ((lng2 - lng1) / (g2 - g1) * g2 - 1.0) * g2 * g2 / (g2 - g1) - g1 / 2.0;
+  w2 = w2 / (g2 - g1);
+  w3 = ((lng3 - lng1) / (g3 - g1) * g3 - 1.0) * g3 * g3 / (g3 - g1) - g1 / 2.0;
+  w3 = w3 / (g3 - g1);
+  w = (w3 - w2) / (g3 - g2);
+  return w;
+}
+
+double libtetrabz_polstat_1233(
+    double g1,
+    double g2,
+    double g3,
+    double lng1,
+    double lng2,
+    double lng3
+) {
+  /*
+  # 3, g4 = g3
+  */
+  double w2, w3, w;
+  w2 = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
+  w2 = (g2 * w2) / (g2 - g1);
+  w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
+  w3 = (g3 * w3) / (g3 - g1);
+  w2 = (w3 - w2) / (g3 - g2);
+  w3 = (lng3 - lng1) / (g3 - g1) * g3 - 1.0;
+  w3 = 1.0 - (2.0 * w3 * g1) / (g3 - g1);
+  w3 = w3 / (g3 - g1);
+  w = (g3 * w3 - g2 * w2) / (g3 - g2);
+  return w;
+}
+
+double libtetrabz_polstat_1221(
+    double g1,
+    double g2,
+    double lng1,
+    double lng2
+) {
+  /*
+  4, g4 = g1 and g3 = g2
+  */
+  double w;
+  w = 1.0 - (lng2 - lng1) / (g2 - g1) * g1;
+  w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
+  w = -1.0 + (3.0 * g2 * w) / (g2 - g1);
+  w = w / (2.0 * (g2 - g1));
+  return w;
+}
+
+double libtetrabz_polstat_1222(
+    double g1,
+    double g2,
+    double lng1,
+    double lng2
+) {
+  /*
+  5, g4 = g3 = g2
+  */
+  double w;
+  w = (lng2 - lng1) / (g2 - g1) * g2 - 1.0;
+  w = (2.0 * g1 * w) / (g2 - g1) - 1.0;
+  w = (3.0 * g1 * w) / (g2 - g1) + 1.0;
+  w = w / (2.0 * (g2 - g1));
+  return w;
+}
+
+double libtetrabz_polstat_1211(
+    double g1,
+    double g2,
+    double lng1,
+    double lng2
+) {
+  /*
+  6, g4 = g3 = g1
+  */
+  double w;
+  w = -1.0 + (lng2 - lng1) / (g2 - g1) * g2;
+  w = -1.0 + (2.0 * g2 * w) / (g2 - g1);
+  w = -1.0 + (3.0 * g2 * w) / (2.0 * (g2 - g1));
+  w = w / (3.0 * (g2 - g1));
+  return w;
+}
+
+void libtetrabz_polstat3(
+    double *de,
+    double *w1
+)
+{
+  /*
+  Tetrahedron method for delta(om - ep + e)
+  */
+  int ii, indx[4];
+  double thr, thr2, e[4], ln[4];
+
+  for (ii = 0; ii < 3; ii++) e[ii] = de[ii];
+  eig_sort(4, e, indx);
+
+  thr = e[3] * 1.0e-3;
+  thr2 = 1.0e-8;
+
+  for(ii=0;ii<4;ii++){
+    if (e[ii] < thr2) {
+      if (ii == 3) {
+        printf("  Nesting # ");
+      }
+      ln[ii] = 0.0;
+      e[ii] = 0.0;
+    }
+    else{
+      ln[ii] = log(e[ii]);
+    }
+  }
+
+  if (fabs(e[3] - e[2]) < thr) {
+    if (fabs(e[3] - e[1]) < thr) {
+      if (fabs(e[3] - e[0]) < thr) {
+        /*
+        e[3] = e[2] = e[1] = e[0]
+        */
+        w1[indx[3]] = 0.25 / e[3];
+        w1[indx[2]] = w1[indx[3]];
+        w1[indx[1]] = w1[indx[3]];
+        w1[indx[0]] = w1[indx[3]];
+      }
+      else {
+        /*
+        e[3] = e[2] = e[1]
+        */
+        w1[indx[3]] = libtetrabz_polstat_1211(e[3], e[0], ln[3], ln[0]);
+        w1[indx[2]] = w1[indx[3]];
+        w1[indx[1]] = w1[indx[3]];
+        w1[indx[0]] = libtetrabz_polstat_1222(e[0], e[3], ln[0], ln[3]);
+        if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
+          printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
+          printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
+          printf("weighting 4=3=2");
+        }
+      }
+    }
+    else if (fabs(e[1] - e[0]) < thr) {
+      /*
+      e[3] = e[2], e[1] = e[0]
+      */
+      w1[indx[3]] = libtetrabz_polstat_1221(e[3], e[1], ln[3], ln[1]);
+      w1[indx[2]] = w1[indx[3]];
+      w1[indx[1]] = libtetrabz_polstat_1221(e[1], e[3], ln[1], ln[3]);
+      w1[indx[0]] = w1[indx[1]];
+
+      if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
+        printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
+        printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
+        printf("weighting 4=3 2=1");
+      }
+    }
+    else {
+      /*
+      e[3] = e[2]
+      */
+      w1[indx[3]] = libtetrabz_polstat_1231(e[3], e[0], e[1], ln[3], ln[0], ln[1]);
+      w1[indx[2]] = w1[indx[3]];
+      w1[indx[1]] = libtetrabz_polstat_1233(e[1], e[0], e[3], ln[1], ln[0], ln[3]);
+      w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[1], e[3], ln[0], ln[1], ln[3]);
+
+      if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
+        printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
+        printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
+        printf("weighting 4=3");
+      }
+    }
+  }
+    else if(fabs(e[2] - e[1]) < thr){
+      if (fabs(e[2] - e[0]) < thr) {
+        /*
+        e[2] = e[1] = e[0]
+        */
+        w1[indx[3]] = libtetrabz_polstat_1222(e[3], e[2], ln[3], ln[2]);
+        w1[indx[2]] = libtetrabz_polstat_1211(e[2], e[3], ln[2], ln[3]);
+        w1[indx[1]] = w1[indx[2]];
+        w1[indx[0]] = w1[indx[2]];
+
+        if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
+          printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
+          printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
+          printf("weighting 3=2=1");
+        }
+      }
+      else {
+        /*
+        e[2] = e[1]
+        */
+        w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[0], e[2], ln[3], ln[0], ln[2]);
+        w1[indx[2]] = libtetrabz_polstat_1231(e[2], e[0], e[3], ln[2], ln[0], ln[3]);
+        w1[indx[1]] = w1[indx[2]];
+        w1[indx[0]] = libtetrabz_polstat_1233(e[0], e[3], e[2], ln[0], ln[3], ln[2]);
+
+        if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
+          printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
+          printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
+          printf("weighting 3=2");
+        }
+      }
+    }
+    else if(fabs(e[1] - e[0]) < thr) {
+      /*
+      e[1] = e[0]
+      */
+      w1[indx[3]] = libtetrabz_polstat_1233(e[3], e[2], e[1], ln[3], ln[2], ln[1]);
+      w1[indx[2]] = libtetrabz_polstat_1233(e[2], e[3], e[1], ln[2], ln[3], ln[1]);
+      w1[indx[1]] = libtetrabz_polstat_1231(e[1], e[2], e[3], ln[1], ln[2], ln[3]);
+      w1[indx[0]] = w1[indx[1]];
+
+      if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
+        printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
+        printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
+        printf("weighting 2=1");
+      }
+    }
+    else {
+    /*
+    Different each other.
+    */
+    w1[indx[3]] = libtetrabz_polstat_1234(e[3], e[0], e[1], e[2], ln[3], ln[0], ln[1], ln[2]);
+    w1[indx[2]] = libtetrabz_polstat_1234(e[2], e[0], e[1], e[3], ln[2], ln[0], ln[1], ln[3]);
+    w1[indx[1]] = libtetrabz_polstat_1234(e[1], e[0], e[2], e[3], ln[1], ln[0], ln[2], ln[3]);
+    w1[indx[0]] = libtetrabz_polstat_1234(e[0], e[1], e[2], e[3], ln[0], ln[1], ln[2], ln[3]);
+
+    if (w1[0] < 0.0 || w1[1] < 0.0 || w1[2] < 0.0 || w1[3] < 0.0) {
+      printf("%f %f %f %f\n", e[0], e[1], e[2], e[3]);
+      printf("%f %f %f %f\n", w1[indx[0]], w1[indx[1]], w1[indx[2]], w1[indx[3]]);
+      printf("weighting");
+    }
+  }
+}
+
+void libtetrabz_polstat2(
+    int nb,
+    double *ei1,
+    double **ej1,
+    double **w1
+) {
+  /*
+  Tetrahedron method for theta( - E2)
+  :param ei1:
+  :param ej1:
+  :return:
+  */
+  int ii, jj, ib, indx[4];
+  double v, thr, e[4], de[4], w2[4], tsmall[4][4];
+
+  thr = 1.0e-8;
+
+  for (ib = 0; ib < nb; ib++) {
+
+    for (ii = 0; ii < 4; ii++) e[ii] = -ej1[ii][ib];
+    eig_sort(4, e, indx);
+
+    if ((e[0] <= 0.0 && 0.0 < e[1]) || (e[0] < 0.0 && 0.0 <= e[1])) {
+
+      libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polstat3(de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
+      }
+    }
+    else if ((e[1] <= 0.0 && 0.0 < e[2]) || (e[1] < 0.0 && 0.0 <= e[2])) {
+
+      libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polstat3(de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polstat3(de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polstat3(de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
+      }
+    }
+    else if ((e[2] <= 0.0 && 0.0 < e[3]) || (e[2] < 0.0 && 0.0 <= e[3])) {
+
+      libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polstat3(de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polstat3(de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
+      }
+
+      libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+
+      if (v > thr) {
+        for (ii = 0; ii < 4; ii++) {
+          de[ii] = 0.0;
+          for (jj = 0; jj < 4; jj++)
+            de[ii] = tsmall[ii][jj] * (ej1[indx[jj]][ib] - ei1[indx[jj]]);
+        }
+        libtetrabz_polstat3(de, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jj = 0; jj < 4; jj++)
+            w1[ib][indx[ii]] += v * w2[jj] * tsmall[jj][ii];
+      }
+    }
+    else if (e[3] <= 0.0) {
+      for (ii = 0; ii < 4; ii++)
+        de[ii] = ej1[indx[ii]][ib] - ei1[indx[ii]];
+      libtetrabz_polstat3(de, w2);
+      for (ii = 0; ii < 4; ii++) w1[ib][ii] += w2[ii];
+    }
+  }
+}
+
+void polstat(
+    int *ng,
+    int nk,
+    int nb,
+    double **bvec,
+    double **eig1,
+    double **eig2,
+    double ***wght
+)
+{
+    /*
+    Compute Static polarization function
+    */
+  int it, ik, ib, ii, jj, jb, **ikv, indx[4];
+  double wlsm[4][20], **ei1, **ej1, *ei2, **ej2, e[4], ***w1, **w2, v, tsmall[4][4], thr;
+
+  thr = 1.0e-10;
+
+  libtetrabz_initialize(ng, bvec, wlsm, ikv);
+
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        wght[ik][ib][jb] = 0.0;
+
+  for(it = 0; it < 6*nk; it++) {
+
+    for (ii = 0; ii < 4; ii++) {
+      for (ib = 0; ib < nb; ib++) {
+        ei1[ii][ib] = 0.0;
+        ej1[ii][ib] = 0.0;
+        for (jj = 0; jj < 20; jj++) {
+          ei1[ii][ib] += wlsm[ii][jj] * eig1[ikv[it][jj]][ib];
+          ej1[ii][ib] += wlsm[ii][jj] * eig2[ikv[it][jj]][ib];
+        }
+      }
+    }
+
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        for (ii = 0; ii < 4; ii++)
+          w1[ib][jb][ii] = 0.0;
+
+    for (ib = 0; ib < nb; ib++) {
+
+      for (ii = 0; ii < 4; ii++) e[ii] = ei1[ii][ib];
+      eig_sort(4, e, indx);
+
+      if (e[0] <= 0.0 && 0.0 < e[1]) {
+
+        libtetrabz_tsmall_a1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_polstat2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[1] <= 0.0 && 0.0 < e[2]) {
+
+        libtetrabz_tsmall_b1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_polstat2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_polstat2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_b3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_polstat2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[2] <= 0.0 && 0.0 < e[3]) {
+
+        libtetrabz_tsmall_c1(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_polstat2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c2(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_polstat2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+
+        libtetrabz_tsmall_c3(e, 0.0, &v, tsmall);
+
+        if (v > thr) {
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++) {
+              ei2[ii] = tsmall[ii][jj] * ej1[indx[jj]][ib];
+              for (jb = 0; jb < nb; jb++)
+                ej2[ii][jb] = tsmall[ii][jj] * ej1[indx[jj]][jb];
+            }
+          libtetrabz_polstat2(nb, ei2, ej2, w2);
+          for (ii = 0; ii < 4; ii++)
+            for (jj = 0; jj < 4; jj++)
+              for (jb = 0; jb < nb; jb++)
+                w1[ib][jb][indx[ii]] += v * w2[jb][jj] * tsmall[jj][ii];
+        }
+      }
+      else if (e[3] <= 0.0) {
+        for (ii = 0; ii < 4; ii++) {
+          ei2[ii] = ej1[ii][ib];
+          for (jb = 0; jb < nb; jb++)
+            ej2[ii][jb] = ej1[ii][jb];
+        }
+        libtetrabz_polstat2(nb, ei2, ej2, w2);
+        for (ii = 0; ii < 4; ii++)
+          for (jb = 0; jb < nb; jb++)
+            w1[ib][jb][ii] += v * w2[jb][ii];
+      }
+      else {
+        continue;
+      }
+    }
+    for (ii = 0; ii < 20; ii++)
+      for (jj = 0; jj < 4; jj++)
+        for (ib = 0; ib < nb; ib++)
+          for (jb = 0; jb < nb; jb++)
+            wght[ikv[it][ii]][ib][jb] += w1[ib][jb][jj] * wlsm[jj][ii];
+  }
+  for (ik = 0; ik < nk; ik++)
+    for (ib = 0; ib < nb; ib++)
+      for (jb = 0; jb < nb; jb++)
+        wght[ik][ib][jb] /= (6.0 * (double) nk);
+}