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