1 files changed, 107 insertions, 0 deletions
diff --git a/src/hh_integral_double.c b/src/hh_integral_double.c
new file mode 100644
index 0000000..2776849
--- /dev/null
+++ b/src/hh_integral_double.c
@@ -0,0 +1,107 @@
+/*
+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 "hh_integral_double.h"
+
+#include <math.h>
+
+#define PI    3.1415926535897932385L
+
+// compute the integral
+int hh_integrate_double(long double* out, hh_params_double params, array_ldouble abcissa, array_ldouble weights){
+  hh_argsint1_double args;
+  int ret;
+
+  args.params=params;
+  args.abcissa=abcissa;
+  args.weights=weights;
+
+  ret=integrate_gauss_ldouble(out, &hh_integrand1_double, -PI/6, PI/6, abcissa, weights, &args);
+
+  return(ret);
+}
+
+// integrand of the integral over theta
+int hh_integrand1_double(long double* out, long double theta, void* args){
+  hh_argsint2_double nargs;
+  int ret;
+  hh_argsint1_double* argument=(hh_argsint1_double*)args;
+
+  nargs.params=argument->params;
+  nargs.theta=theta;
+
+  ret=integrate_gauss_ldouble(out, &hh_integrand2_double, 0, 1./cos(theta-PI/6), argument->abcissa, argument->weights, &nargs);
+
+  return(ret);
+}
+
+// integrand of the integral over rho
+int hh_integrand2_double(long double* out, long double rho, void* args){
+  hh_argsint2_double* argument=(hh_argsint2_double*)args;
+  long double m, xi, alpha2, O2;
+
+  alpha2=-2*sinl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))+cosl(PI/sqrtl(3.)*rho*cosl(argument->theta)));
+  O2=1+4*cosl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))-cosl(PI/sqrtl(3.)*rho*cosl(argument->theta)));
+  m=argument->params.W-2*argument->params.t2*argument->params.sinphi*alpha2;
+  xi=sqrtl(m*m+argument->params.t1*argument->params.t1*O2);
+  *out=rho*m/xi;
+
+  return(0);
+}
+
+
+// derivative
+int hh_d_integrate_double(long double* out, hh_params_double params, array_ldouble abcissa, array_ldouble weights){
+  hh_argsint1_double args;
+  int ret;
+
+  args.params=params;
+  args.abcissa=abcissa;
+  args.weights=weights;
+
+  ret=integrate_gauss_ldouble(out, &hh_d_integrand1_double, -PI/6, PI/6, abcissa, weights, &args);
+
+  return(ret);
+}
+
+// derivative of the integrand of the integral over theta
+int hh_d_integrand1_double(long double* out, long double theta, void* args){
+  hh_argsint2_double nargs;
+  int ret;
+  hh_argsint1_double* argument=(hh_argsint1_double*)args;
+
+  nargs.params=argument->params;
+  nargs.theta=theta;
+
+  ret=integrate_gauss_ldouble(out, &hh_d_integrand2_double, 0, 1./cos(theta-PI/6), argument->abcissa, argument->weights, &nargs);
+
+  return(ret);
+}
+
+// derivative of the integrand of the integral over rho
+int hh_d_integrand2_double(long double* out, long double rho, void* args){
+  hh_argsint2_double* argument=(hh_argsint2_double*)args;
+  long double m, xi, alpha2, O2;
+
+  alpha2=-2*sinl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))+cosl(PI/sqrtl(3.)*rho*cosl(argument->theta)));
+  O2=1+4*cosl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))-cosl(PI/sqrtl(3.)*rho*cosl(argument->theta)));
+  m=argument->params.W-2*argument->params.t2*argument->params.sinphi*alpha2;
+  xi=sqrtl(m*m+argument->params.t1*argument->params.t1*O2);
+  *out=rho/xi*(1-m*m/xi/xi);
+
+  return(0);
+}
+