Ian Jauslin

Ian Jauslin

Home Publications Talks Software Teaching Contact
summaryrefslogtreecommitdiff
path: root/src/ss_integral_double.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/ss_integral_double.c')
-rw-r--r--src/ss_integral_double.c437
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));
+}
generated by cgit v1.2.3-54-g00ecf (git 2.44.0) at 2024-04-28 04:41:23 +0000