/*
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.h"

#include <stdarg.h>
// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears
#define MPFR_USE_VA_LIST
#include <mpfr.h>
#include <libinum.h>

int hh_integrate(mpfr_t* out, hh_params params, array_mpfr abcissa, array_mpfr weights){
  mpfr_t lower, upper;
  // arguments for first integral
  hh_argsint1 args1;
  int ret;

  args1.params=params;
  args1.abcissa=abcissa;
  args1.weights=weights;

  mpfr_inits(lower, upper, NULL);

  // compute the boundaries of the integral over theta
  // pi/6
  mpfr_const_pi(upper, MPFR_RNDN);
  mpfr_div_ui(upper, upper, 6, MPFR_RNDN);
  // -pi/6
  mpfr_neg(lower, upper, MPFR_RNDN);

  // integrate
  ret=integrate_gauss_mpfr(out, &hh_integrand1, lower, upper, abcissa, weights, &args1);

  mpfr_clears(lower, upper, NULL);

  return(ret);
}

// integrand of the integral over theta
int hh_integrand1(mpfr_t* out, mpfr_t theta, void* args){
  hh_argsint2 args2;
  mpfr_t lower, upper;
  int ret;

  mpfr_inits(upper, lower, NULL);

  // recover parameters
  args2.params=((hh_argsint1*)args)->params;
  args2.theta=theta;

  // boundaries of the integral over rho
  // 1/cos(theta-pi/6)
  mpfr_const_pi(upper, MPFR_RNDN);
  mpfr_div_ui(upper, upper, 6, MPFR_RNDN);
  mpfr_sub(upper, theta, upper, MPFR_RNDN);
  mpfr_cos(upper, upper, MPFR_RNDN);
  mpfr_ui_div(upper, 1, upper, MPFR_RNDN);

  // 0
  mpfr_set_ui(lower, 0, MPFR_RNDN);

  ret=integrate_gauss_mpfr(out, &hh_integrand2, lower, upper, ((hh_argsint1*)args)->abcissa, ((hh_argsint1*)args)->weights, &args2);

  mpfr_clears(upper, lower, NULL);

  return(ret);
}

// integrand of the integral over rho
int hh_integrand2(mpfr_t* out, mpfr_t rho, void* args){
  mpfr_t tmp1, tmp2, tmp3;

  mpfr_inits(tmp1, tmp2, tmp3, NULL);

  // out = Omega^2
  // tmp1 = alpha2
  hh_Omega2_alpha2(*out, tmp1, rho, ((hh_argsint2*)args)->theta, tmp2, tmp3);
  // tmp1 = m
  hh_m(tmp1, ((hh_argsint2*)args)->params.W, ((hh_argsint2*)args)->params.t2, ((hh_argsint2*)args)->params.sinphi, tmp1);
  // out = xi^2
  hh_xi2(*out, *out, tmp1, ((hh_argsint2*)args)->params.t1, tmp2);

  // out = rho/sqrt(xi^2)*m
  mpfr_sqrt(*out, *out, MPFR_RNDN);
  mpfr_div(*out, rho, *out, MPFR_RNDN);
  mpfr_mul(*out, *out, tmp1, MPFR_RNDN);

  mpfr_clears(tmp1, tmp2, tmp3, NULL);

  return(0);
}

// derivative
int hh_d_integrate(mpfr_t* out, hh_params params, array_mpfr abcissa, array_mpfr weights){
  mpfr_t lower, upper;
  // arguments for first integral
  hh_argsint1 args1;
  int ret;

  args1.params=params;
  args1.abcissa=abcissa;
  args1.weights=weights;

  mpfr_inits(lower, upper, NULL);

  // compute the boundaries of the integral over theta
  // pi/6
  mpfr_const_pi(upper, MPFR_RNDN);
  mpfr_div_ui(upper, upper, 6, MPFR_RNDN);
  // -pi/6
  mpfr_neg(lower, upper, MPFR_RNDN);

  // integrate
  ret=integrate_gauss_mpfr(out, &hh_d_integrand1, lower, upper, abcissa, weights, &args1);

  mpfr_clears(lower, upper, NULL);

  return(ret);
}

// derivative of the integrand of the integral over theta
int hh_d_integrand1(mpfr_t* out, mpfr_t theta, void* args){
  hh_argsint2 args2;
  mpfr_t lower, upper;
  int ret;

  mpfr_inits(lower, upper, NULL);

  // recover parameters
  args2.params=((hh_argsint1*)args)->params;
  args2.theta=theta;

  // boundaries of the integral over rho
  // 1/cos(theta-pi/6)
  mpfr_const_pi(upper, MPFR_RNDN);
  mpfr_div_ui(upper, upper, 6, MPFR_RNDN);
  mpfr_sub(upper, theta, upper, MPFR_RNDN);
  mpfr_cos(upper, upper, MPFR_RNDN);
  mpfr_ui_div(upper, 1, upper, MPFR_RNDN);

  // 0
  mpfr_set_ui(lower, 0, MPFR_RNDN);

  ret=integrate_gauss_mpfr(out, &hh_d_integrand2, lower, upper, ((hh_argsint1*)args)->abcissa, ((hh_argsint1*)args)->weights, &args2);

  mpfr_clears(lower, upper, NULL);

  return(ret);
}

// derivative of the integrand of the integral over rho
int hh_d_integrand2(mpfr_t* out, mpfr_t rho, void* args){
  mpfr_t tmp1, tmp2, tmp3;

  mpfr_inits(tmp1, tmp2, tmp3, NULL);

  // out = Omega^2
  // tmp1 = alpha2
  hh_Omega2_alpha2(*out, tmp1, rho, ((hh_argsint2*)args)->theta, tmp2, tmp3);
  // tmp1 = m
  hh_m(tmp1, ((hh_argsint2*)args)->params.W, ((hh_argsint2*)args)->params.t2, ((hh_argsint2*)args)->params.sinphi, tmp1);
  // out = xi^2
  hh_xi2(*out, *out, tmp1, ((hh_argsint2*)args)->params.t1, tmp2);

  // tmp2 = 1-m^2/xi^2
  mpfr_pow_ui(tmp2, tmp1, 2,MPFR_RNDN);
  mpfr_div(tmp2, tmp2, *out, MPFR_RNDN);
  mpfr_ui_sub(tmp2, 1, tmp2, MPFR_RNDN);

  // out = rho/sqrt(xi^2)*(1-m^2/xi^2)
  mpfr_sqrt(*out, *out, MPFR_RNDN);
  mpfr_div(*out, rho, *out, MPFR_RNDN);
  mpfr_mul(*out, *out, tmp2, MPFR_RNDN);

  mpfr_clears(tmp1, tmp2, tmp3, NULL);

  return(0);
}

// Omega^2 and alpha_2
// provide two initialized tmp mpfr_t's
// Omega2 and alpha2 must be initialized
int hh_Omega2_alpha2(mpfr_t Omega2, mpfr_t alpha2, mpfr_t rho, mpfr_t theta, mpfr_t tmp1, mpfr_t tmp2){
  // Omega2 and alpha2 will be used as tmp variables whenever possible

  // Omega2 = pi
  mpfr_const_pi(Omega2, MPFR_RNDN);

  // tmp1 = pi/sqrt(3)*rho
  mpfr_sqrt_ui(tmp1, 3, MPFR_RNDN);
  mpfr_div(tmp1, Omega2, tmp1, MPFR_RNDN);
  mpfr_mul(tmp1, tmp1, rho, MPFR_RNDN);

  // alpha2 = cos(theta)
  mpfr_cos(alpha2, theta, MPFR_RNDN);

  // tmp1 = cos(pi/sqrt(3)*rho*cos(theta))
  mpfr_mul(tmp1, tmp1, alpha2, MPFR_RNDN);
  //// alpha2 free
  mpfr_cos(tmp1, tmp1, MPFR_RNDN);

  // alpha2 = pi/3*(1+rho*sin(theta))
  mpfr_sin(alpha2, theta, MPFR_RNDN);
  mpfr_mul(alpha2, alpha2, rho, MPFR_RNDN);
  mpfr_add_ui(alpha2, alpha2, 1, MPFR_RNDN);
  mpfr_div_ui(alpha2, alpha2, 3, MPFR_RNDN);
  mpfr_mul(alpha2, alpha2, Omega2, MPFR_RNDN);
  //// Omega2 free

  // Omega2 = cos(pi/3*(1+rho*sin(theta)))
  mpfr_cos(Omega2, alpha2, MPFR_RNDN);

  // tmp2 = sin(pi/3*(1+rho*sin(theta)))
  mpfr_sin(tmp2, alpha2, MPFR_RNDN);
  //// alpha2 free

  // alpha2 = -2*sin(pi/3*(1+rho*sin(theta)))*(cos(pi/3*(1+rho*sin(theta)))+cos(pi/sqrt(3)*rho*cos(theta)))
  mpfr_add(alpha2, Omega2, tmp1, MPFR_RNDN);
  mpfr_mul(alpha2, alpha2, tmp2, MPFR_RNDN);
  //// tmp2 free
  mpfr_mul_si(alpha2, alpha2, -2, MPFR_RNDN);

  // tmp1 = cos(pi/3*(1+rho*sin(theta)))-cos(pi/sqrt(3)*rho*cos(theta))
  mpfr_sub(tmp1, Omega2, tmp1, MPFR_RNDN);
  // Omega2 = 1+4*cos(pi/3*(1+rho*sin(theta)))*(cos(pi/3*(1+rho*sin(theta)))-cos(pi/sqrt(3)*rho*cos(theta)))
  mpfr_mul(Omega2, Omega2, tmp1, MPFR_RNDN);
  //// tmp1 free
  mpfr_mul_ui(Omega2, Omega2, 4, MPFR_RNDN);
  mpfr_add_ui(Omega2, Omega2, 1, MPFR_RNDN);

  return(0);
}

// m
// out must be initialized
// out and alpha2 can point to the same number
int hh_m(mpfr_t out, mpfr_t W, mpfr_t t2, mpfr_t sinphi, mpfr_t alpha2){
  // out = W-2*t2*sinphi*alpha2
  mpfr_mul(out, alpha2, sinphi, MPFR_RNDN);
  mpfr_mul(out, out, t2, MPFR_RNDN);
  mpfr_mul_ui(out, out, 2, MPFR_RNDN);
  mpfr_sub(out, W, out, MPFR_RNDN);
  return(0);
}

// xi^2
// provide one initialized tmp mpfr_t
// out must be initialized
// out and Omega2 can point to the same number
// tmp and m can point to the same number
int hh_xi2(mpfr_t out, mpfr_t Omega2, mpfr_t m, mpfr_t t1, mpfr_t tmp){
  // out = t1^2*Omega^2
  mpfr_mul(out, Omega2, t1, MPFR_RNDN);
  mpfr_mul(out, out, t1, MPFR_RNDN);

  // tmp = m^2
  mpfr_pow_ui(tmp, m, 2, MPFR_RNDN);

  // out = m^2+t1^2*Omega^2
  mpfr_add(out, out, tmp, MPFR_RNDN);

  return(0);
}