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