diff options
Diffstat (limited to 'src/hh_integral.c')
| -rw-r--r-- | src/hh_integral.c | 276 |
1 files changed, 276 insertions, 0 deletions
diff --git a/src/hh_integral.c b/src/hh_integral.c new file mode 100644 index 0000000..f8eda10 --- /dev/null +++ b/src/hh_integral.c @@ -0,0 +1,276 @@ +/* +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); +} |