diff options
author | Ian Jauslin <ian.jauslin@roma1.infn.it> | 2016-05-24 13:39:23 +0000 |
---|---|---|
committer | Ian Jauslin <ian.jauslin@roma1.infn.it> | 2016-05-24 13:39:23 +0000 |
commit | fa9b6f2b9bcb80778e63ef2aa4e17c7573de0015 (patch) | |
tree | 92b740d0736c9ed6f5bda051c224c8bb7196bb03 /src/ss_integral_double.c |
Diffstat (limited to 'src/ss_integral_double.c')
-rw-r--r-- | src/ss_integral_double.c | 437 |
1 files changed, 437 insertions, 0 deletions
diff --git a/src/ss_integral_double.c b/src/ss_integral_double.c new file mode 100644 index 0000000..00ab637 --- /dev/null +++ b/src/ss_integral_double.c @@ -0,0 +1,437 @@ +/* +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 "ss_integral_double.h" +#include <math.h> + +#define PI 3.1415926535897932385L +#define SQRT3 1.7320508075688772935L + +// compute the integral +int ss_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + long double val; + + // compute A_FF, A_RF and A_RR + ssFF_integrate_double(out, I, params, abcissa, weights, threads); + ssRF_integrate_double(&val, I, params, abcissa, weights, threads); + *out+=2*val; + ssRR_integrate_double(&val, I, params, abcissa, weights, threads); + *out+=val; + + return(0); +} + + +// compute the integral A_FF +int ssFF_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + ssFF_argsint_rho_double args; + int ret; + array_pthread_t thread_ids; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + + if(threads==1){ + ret=integrate_gauss_ldouble(out, &ssFF_integrand_rho_double, 0., 1., abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_ldouble(out, &ssFF_integrand_rho_double, 0., 1., abcissa, weights, &args, threads, &thread_ids); + array_pthread_t_free(thread_ids); + } + + // extra factor + *out*=6; + + return(ret); +} + +// integrand of the integral over rho in A_FF +int ssFF_integrand_rho_double(long double* out, long double rho, void* args){ + ssFF_argsint_theta_double nargs; + int ret; + ssFF_argsint_rho_double* argument=(ssFF_argsint_rho_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.rho=rho; + + ret=integrate_gauss_ldouble(out, &ssFF_integrand_theta_double, -PI/6, PI/2, argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over theta in A_FF +int ssFF_integrand_theta_double(long double* out, long double theta, void* args){ + ssFF_argsint_psi_double nargs; + int ret; + ssFF_argsint_theta_double* argument=(ssFF_argsint_theta_double*)args; + long double val; + int a; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.rho=argument->rho; + nargs.theta=theta; + + // split the integral into four parts to improve the precision + *out=0; + for(a=0;a<4;a++){ + ret=integrate_gauss_ldouble(&val, &ssFF_integrand_psi_double, PI/4*(a-2), PI/4*(a-1), argument->abcissa, argument->weights, &nargs); + if(ret<0){ + break; + } + *out+=val; + } + + return(ret); +} + +// integrand of the integral over psi in A_FF +int ssFF_integrand_psi_double(long double* out, long double psi, void* args){ + ssFF_argsint_z_double nargs; + int ret; + ssFF_argsint_psi_double* argument=(ssFF_argsint_psi_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.rho=argument->rho; + nargs.theta=argument->theta; + nargs.psi=psi; + + ret=integrate_gauss_ldouble(out, &ssFF_integrand_z_double, 0., 1., argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over z in A_FF +int ssFF_integrand_z_double(long double* out, long double z, void* args){ + ssFF_argsint_z_double* argument=(ssFF_argsint_z_double*)args; + long double p1,p2,q1,q2,F1,F2; + + // r and varphi (store in p1 and p2) + p1=(1-z)*(1+z*(1-cosl(argument->psi))/(1+cosl(argument->psi))); + p2=acosl(cos(2*argument->psi)-p1/2*sin(2*argument->psi)*sin(2*argument->psi)); + if(argument->psi<0.){ + p2=-p2; + } + + // jacobian + *out=4*argument->rho*argument->rho*argument->rho*p1*(1+p1*cosl(2*argument->psi))/(1+cosl(argument->psi))/sqrtl(1+p1*cosl(argument->psi)*cosl(argument->psi)); + + // cutoffs + *out*=ss_cutoff_double(argument->rho)*ss_cutoff_double(argument->rho*p1); + + q1=SQRT3+argument->rho*p1*cosl(argument->theta+p2); + q2=argument->params.omega*(1+argument->rho*p1*sinl(argument->theta+p2)); + p1=SQRT3+argument->rho*cosl(argument->theta); + p2=argument->params.omega*(1+argument->rho*sinl(argument->theta)); + F1=p1+q1-SQRT3; + F2=p2+q2-argument->params.omega; + + *out*=(*(argument->I))(p1, p2, q1, q2, F1, F2, argument->params.t1, argument->params.t2, argument->params.phi, argument->params.W); + + return(0); +} + + +// compute the integral A_RF +int ssRF_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + ssRF_argsint_theta_double args; + int ret; + long double val; + int a; + array_pthread_t thread_ids; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + + // split the integral in three parts to improve the precision + *out=0; + for(a=0;a<3;a++){ + if(threads==1){ + ret=integrate_gauss_ldouble(&val, &ssRF_integrand_theta_double, (4*a-1)*PI/6, (4*a+3)*PI/6, abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_ldouble(&val, &ssRF_integrand_theta_double, (4*a-1)*PI/6, (4*a+3)*PI/6, abcissa, weights, &args, threads, &thread_ids); + array_pthread_t_free(thread_ids); + } + if(ret<0){ + break; + } + + *out+=val; + } + + // extra factor + *out*=3; + + return(ret); +} + +// integrand of the integral over theta in A_RF +int ssRF_integrand_theta_double(long double* out, long double theta, void* args){ + ssRF_argsint_varphi_double nargs; + int ret; + ssRF_argsint_theta_double* argument=(ssRF_argsint_theta_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.theta=theta; + + ret=integrate_gauss_ldouble(out, &ssRF_integrand_varphi_double, -PI/6, PI/2, argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over varphi in A_RF +int ssRF_integrand_varphi_double(long double* out, long double varphi, void* args){ + ssRF_argsint_r_double nargs; + int ret; + ssRF_argsint_varphi_double* argument=(ssRF_argsint_varphi_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.theta=argument->theta; + nargs.varphi=varphi; + + ret=integrate_gauss_ldouble(out, &ssRF_integrand_r_double, 0., 1., argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over r in A_RF +int ssRF_integrand_r_double(long double* out, long double r, void* args){ + ssRF_argsint_rho_double nargs; + int ret; + ssRF_argsint_r_double* argument=(ssRF_argsint_r_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.theta=argument->theta; + nargs.varphi=argument->varphi; + nargs.r=r; + + ret=integrate_gauss_ldouble(out, &ssRF_integrand_rho_double, 0., ss_R_double(nargs.varphi, 0), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over rho in A_RF +int ssRF_integrand_rho_double(long double* out, long double rho, void* args){ + int nu; + ssRF_argsint_rho_double* argument=(ssRF_argsint_rho_double*)args; + long double p1,p2,q1,q2,F1,F2; + + *out=0.; + for(nu=-1;nu<=1;nu=nu+2){ + q1=SQRT3+argument->r*cosl(argument->theta); + q2=argument->params.omega*(1+argument->r*sinl(argument->theta)); + F1=SQRT3+rho*cosl(argument->varphi); + F2=nu*(1+rho*sinl(argument->varphi)); + p1=SQRT3+F1-q1; + p2=argument->params.omega+F2-q2; + + *out+=argument->r*rho*(*(argument->I))(p1, p2, q1, q2, F1, F2, argument->params.t1, argument->params.t2, argument->params.phi, argument->params.W)*(1-ss_cutoff_double(ss_norm_double(F1-q1,F2-q2)))*ss_cutoff_double(argument->r); + } + + return(0); +} + + +// compute the integral A_RR +int ssRR_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + ssRR_argsint_theta_double args; + int ret; + int a; + long double lower, upper; + long double val; + array_pthread_t thread_ids; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + + // split the integral in three parts + *out=0; + for(a=0;a<3;a++){ + lower=(4*a-1)*PI/6; + upper=(4*a+3)*PI/6; + + // save sector + args.sector_theta=a; + + if(threads==1){ + ret=integrate_gauss_ldouble(&val, &ssRR_integrand_theta_double, lower, upper, abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_ldouble(&val, &ssRR_integrand_theta_double, lower, upper, abcissa, weights, &args, threads, &thread_ids); + array_pthread_t_free(thread_ids); + } + if(ret<0){ + break; + } + + *out+=val; + } + // extra factor + *out*=3; + + return(ret); +} + +// integrand of the integral over theta in A_RR +int ssRR_integrand_theta_double(long double* out, long double theta, void* args){ + ssRR_argsint_varphi_double nargs; + int ret; + ssRR_argsint_theta_double* argument=(ssRR_argsint_theta_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.sector_theta=argument->sector_theta; + nargs.theta=theta; + + ret=integrate_gauss_ldouble(out, &ssRR_integrand_varphi_double, -PI/6, PI/2, argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over varphi in A_RR +int ssRR_integrand_varphi_double(long double* out, long double varphi, void* args){ + ssRR_argsint_r_double nargs; + int ret; + ssRR_argsint_varphi_double* argument=(ssRR_argsint_varphi_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.theta=argument->theta; + nargs.varphi=varphi; + + ret=integrate_gauss_ldouble(out, &ssRR_integrand_r_double, 0., ss_R_double(nargs.theta, argument->sector_theta), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over r in A_RR +int ssRR_integrand_r_double(long double* out, long double r, void* args){ + ssRR_argsint_rho_double nargs; + int ret; + ssRR_argsint_r_double* argument=(ssRR_argsint_r_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.theta=argument->theta; + nargs.varphi=argument->varphi; + nargs.r=r; + + ret=integrate_gauss_ldouble(out, &ssRR_integrand_rho_double, 0., ss_R_double(nargs.varphi, 0), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over rho in A_RR +int ssRR_integrand_rho_double(long double* out, long double rho, void* args){ + int eta, nu; + ssRR_argsint_rho_double* argument=(ssRR_argsint_rho_double*)args; + long double p1,p2,q1,q2,F1,F2; + + *out=0.; + for(eta=-1;eta<=1;eta=eta+2){ + for(nu=-1;nu<=1;nu=nu+2){ + q1=SQRT3+argument->r*cosl(argument->theta); + q2=eta*(1+argument->r*sinl(argument->theta)); + F1=SQRT3+rho*cosl(argument->varphi); + F2=nu*(1+rho*sinl(argument->varphi)); + p1=SQRT3+F1-q1; + p2=argument->params.omega+F2-q2; + + *out+=argument->r*rho*(*(argument->I))(p1, p2, q1, q2, F1, F2, argument->params.t1, argument->params.t2, argument->params.phi, argument->params.W)*(1-ss_cutoff_double(ss_norm_double(F1-q1,F2-q2)))*(1-ss_cutoff_double(ss_norm_double(q1-SQRT3,q2-argument->params.omega))); + } + } + + return(0); +} + + +// R(theta) +long double ss_R_double(long double theta, int sector){ + return(1./cosl(theta+(8*sector-1)*PI/6)); +} + +// cutoff function +long double ss_cutoff_double(long double x){ + long double e1,e2; + + if(x<=0.5){ + return(1.); + } + if(x>=1){ + return(0.); + } + + e1=expl(-0.5/(1-x)); + e2=expl(-0.5/(x-0.5)); + return(e1/(e2+e1*(1-e2))); +} + +// periodic norm +long double ss_norm_double(long double k1, long double k2){ + long double n1, n2, n3; + long double tmp1, tmp2; + long double t1, t2; + + // split R^2 into equilateral triangles and find which triangle (k1,k2) is in + n1=floorl(k1/SQRT3)+1; + n2=floorl(k1/2/SQRT3-k2/2)+1; + n3=floorl(k1/2/SQRT3+k2/2)+1; + + // translate (k1,k2) by -a*(SQRT3,3)-b(SQRT3,-3) + // a-b + tmp1=-floorl((n2-n3+1)/3); + // b + tmp2=floorl((n1-tmp1)/2); + // a + tmp1=tmp1+tmp2; + t1=k1-SQRT3*(tmp1+tmp2); + t2=k2-3*(tmp1-tmp2); + + return(sqrtl(t1*t1+t2*t2)); +} |