elliptics.hpp

Go to the documentation of this file.

#ifndef ELLIPTICS_HPP
#define ELLIPTICS_HPP
// This function is for evaluating elliptic integrals
 
#include <Kokkos_Core.hpp>
#include "constants.hpp"
#include "space_settings.hpp"
 
//> ellip_agm & ellip_agm_v
// To calculate complete elliptic integrals
// Using arithmetic-geometric mean (agm) and modified agm (magm)
// ellip_agm by E.S. Yoon - Mar. 09, 2015 : original code
// ellip_agm_v by Nathan Wichmann - June, 2015 : vectorized ellip_agm
// Note that ellip_agm_v is NOT for OpenACC
KOKKOS_INLINE_FUNCTION void ellip_agm(double k,double& ellip_k,double& ellip_e,int& order){
    constexpr int n_order=10;
    constexpr double e_tol=1e-8;
 
    order = 10;
    double mx_n = 1.0;
    double my_n = 1.0 - k ; //beta^2
    double mz_n = 0.0;
    double x_n = 1.0;
    double y_n = sqrt(my_n); //beta
 
    for (int i=0; i<n_order; i++){
        //agm
        double x_np1=(x_n+y_n)*0.5;
        y_n=sqrt(x_n*y_n);
        x_n=x_np1; //update results
 
        //magm
        double mx_np1=(mx_n+my_n)*0.5;
        double rd=sqrt((mx_n-mz_n)*(my_n-mz_n));
        double my_np1=mz_n+rd;
        mz_n=mz_n-rd;
        mx_n=mx_np1;
        my_n=my_np1;
 
        if((abs(x_n-y_n) < e_tol) && (abs(mx_n-my_n) < e_tol)){
            order = i+1;
            break;
        }
    }
 
    ellip_k=0.5*PI/x_n;
    ellip_e=mx_n*ellip_k;
}
 
#endif