/*
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);
}