XGC1
elliptics.hpp
Go to the documentation of this file.
1 #ifndef ELLIPTICS_HPP
2 #define ELLIPTICS_HPP
3 // This function is for evaluating elliptic integrals
4
5 #include <Kokkos_Core.hpp>
6
7 //> ellip_agm & ellip_agm_v
8 // To calculate complete elliptic integrals
9 // Using arithmetic-geometric mean (agm) and modified agm (magm)
10 // ellip_agm by E.S. Yoon - Mar. 09, 2015 : original code
11 // ellip_agm_v by Nathan Wichmann - June, 2015 : vectorized ellip_agm
12 // Note that ellip_agm_v is NOT for OpenACC
13 KOKKOS_INLINE_FUNCTION void ellip_agm(double k,double& ellip_k,double& ellip_e,int& vpic_ierr){
14  const int n_order=10;
15  const double e_tol=1e-8;
16
17  int order = 10;
18  double mx_n = 1.0;
19  double my_n = 1.0 - k ; //beta^2
20  double mz_n = 0.0;
21  double x_n = 1.0;
22  double y_n = sqrt(my_n); //beta
23
24  for (int i=0; i<n_order; i++){
25  //agm
26  double x_np1=(x_n+y_n)*0.5;
27  y_n=sqrt(x_n*y_n);
28  x_n=x_np1; //update results
29
30  //magm
31  double mx_np1=(mx_n+my_n)*0.5;
32  double rd=sqrt((mx_n-mz_n)*(my_n-mz_n));
33  double my_np1=mz_n+rd;
34  mz_n=mz_n-rd;
35  mx_n=mx_np1;
36  my_n=my_np1;
37
38  if((abs(x_n-y_n) < e_tol) && (abs(mx_n-my_n) < e_tol)){
39  order = i+1;
40  break;
41  }
42  }
43
44  ellip_k=0.5*PI/x_n;
45  ellip_e=mx_n*ellip_k;
46  vpic_ierr=order;
47 }
48
49 #endif
KOKKOS_INLINE_FUNCTION void ellip_agm(double k, double &ellip_k, double &ellip_e, int &vpic_ierr)
Definition: elliptics.hpp:13
constexpr double PI
Definition: constants.hpp:8