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>

//> 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& vpic_ierr){
    const int n_order=10;
    const double e_tol=1e-8;

    int 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;
    vpic_ierr=order;
}

#endif