1 files changed, 89 insertions, 0 deletions
diff --git a/src/zz_integral_double.c b/src/zz_integral_double.c
new file mode 100644
index 0000000..9330ad1
--- /dev/null
+++ b/src/zz_integral_double.c
@@ -0,0 +1,89 @@
+/*
+Copyright 2016 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "zz_integral_double.h"
+#include <math.h>
+
+#define PI    3.1415926535897932385L
+#define SQRT3 1.7320508075688772935L
+
+// I for z1-z2
+long double zz_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W){
+  long double xip,xiq,xiF,mp,mq,mF,zetap,zetaq,zetaF;
+
+  // alpha1 (store in xi)
+  xip=zz_alpha1_double(p1,p2);
+  xiq=zz_alpha1_double(q1,q2);
+  xiF=zz_alpha1_double(F1,F2);
+
+  // alpha2 (store in m)
+  mp=zz_alpha2_double(p1,p2);
+  mq=zz_alpha2_double(q1,q2);
+  mF=zz_alpha2_double(F1,F2);
+
+  zetap=2*t2*cosl(phi)*xip;
+  zetaq=2*t2*cosl(phi)*xiq;
+  zetaF=2*t2*cosl(phi)*xiF;
+
+  mp=W-2*t2*sinl(phi)*mp;
+  mq=W-2*t2*sinl(phi)*mq;
+  mF=W-2*t2*sinl(phi)*mF;
+
+  xip=sqrtl(mp*mp+2*t1*t1*xip);
+  xiq=sqrtl(mq*mq+2*t1*t1*xiq);
+  xiF=sqrtl(mF*mF+2*t1*t1*xiF);
+
+  return((xip+xiq+xiF)*(mp/xip+mq/xiq-mF/xiF-mp*mq*mF/xip/xiq/xiF)*(zetap+zetaq-zetaF)/(((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF))*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF)))/108);
+}
+
+// I for z1+z2
+long double ZZ_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W){
+  long double xip,xiq,xiF,mp,mq,mF,zetap,zetaq,zetaF;
+
+  // alpha1 (store in xi)
+  xip=zz_alpha1_double(p1,p2);
+  xiq=zz_alpha1_double(q1,q2);
+  xiF=zz_alpha1_double(F1,F2);
+
+  // alpha2 (store in m)
+  mp=zz_alpha2_double(p1,p2);
+  mq=zz_alpha2_double(q1,q2);
+  mF=zz_alpha2_double(F1,F2);
+
+  zetap=2*t2*cosl(phi)*xip;
+  zetaq=2*t2*cosl(phi)*xiq;
+  zetaF=2*t2*cosl(phi)*xiF;
+
+  mp=W-2*t2*sinl(phi)*mp;
+  mq=W-2*t2*sinl(phi)*mq;
+  mF=W-2*t2*sinl(phi)*mF;
+
+  xip=sqrtl(mp*mp+2*t1*t1*xip);
+  xiq=sqrtl(mq*mq+2*t1*t1*xiq);
+  xiF=sqrtl(mF*mF+2*t1*t1*xiF);
+
+  return((1-mp/xip*mF/xiF-mq/xiq*mF/xiF+mp/xip*mq/xiq)*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)+(xip+xiq+xiF)*(xip+xiq+xiF))/(((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF))*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF)))/216);
+}
+
+// bar alpha_1(k1,k2)
+long double zz_alpha1_double(long double k1, long double k2){
+  return(2*cosl(PI/3*k2)*(cosl(PI/SQRT3*k1)+cosl(PI/3*k2))+0.5);
+}
+// bar alpha_2(k1,k2)
+long double zz_alpha2_double(long double k1, long double k2){
+  return(2*sinl(PI/3*k2)*(cosl(PI/SQRT3*k1)-cosl(PI/3*k2)));
+}
+