<commit_before>#pragma once

#include "dg/topology/grid.h"
#include "dg/topology/functions.h"
#include "dg/topology/interpolation.h"
#include "dg/topology/derivatives.h"
#include "dg/topology/geometry.h"
#include "dg/functors.h"
#include "dg/runge_kutta.h"
#include "dg/nullstelle.h"
#include "fluxfunctions.h"
#include "ribeiro.h"


namespace dg
{
namespace geo
{

///@cond
namespace flux
{
namespace detail
{

//This leightweights struct and its methods finds the initial R and Z values and the coresponding f(\psi) as
//good as it can, i.e. until machine precision is reached
//Note that f(psi) = 1/q(psi) (The safety factor)
struct Fpsi
{

    //firstline = 0 -> conformal, firstline = 1 -> equalarc
    Fpsi( const CylindricalFunctorsLvl2& psip, const CylindricalFunctorsLvl1& ipol, double x0, double y0, bool verbose = false):
        psip_(psip), fieldRZYT_(psip, ipol, x0, y0), fieldRZtau_(psip),m_verbose(verbose)
    {
        //Find O-point
        double R_O = x0, Z_O = y0;
        dg::geo::findOpoint( psip, R_O, Z_O);
        //define angle with respect to O-point
        fieldRZYT_ = dg::geo::flux::FieldRZYT(psip, ipol, R_O, Z_O);
        X_init = x0, Y_init = y0;
        while( fabs( psip.dfx()(X_init, Y_init)) <= 1e-10 && fabs( psip.dfy()( X_init, Y_init)) <= 1e-10)
            X_init +=  1.;
    }
    //finds the starting points for the integration in y direction
    void find_initial( double psi, double& R_0, double& Z_0)
    {
        unsigned N = 50;
        std::array<double, 2> begin2d{ {0,0} }, end2d(begin2d), end2d_old(begin2d);
        if(m_verbose)std::cout << "In init function\n";
        begin2d[0] = end2d[0] = end2d_old[0] = X_init;
        begin2d[1] = end2d[1] = end2d_old[1] = Y_init;
        double eps = 1e10, eps_old = 2e10;
        while( (eps < eps_old || eps > 1e-7) && eps > 1e-14)
        {
            eps_old = eps; end2d_old = end2d;
            N*=2; dg::stepperRK( "Feagin-17-8-10",  fieldRZtau_, psip_.f()(X_init, Y_init), begin2d, psi, end2d, N);
            eps = sqrt( (end2d[0]-end2d_old[0])*(end2d[0]-end2d_old[0]) + (end2d[1]-end2d_old[1])*(end2d[1]-end2d_old[1]));
        }
        X_init = R_0 = end2d_old[0], Y_init = Z_0 = end2d_old[1];
        if(m_verbose)std::cout << "In init function error: psi(R,Z)-psi0: "<<psip_.f()(X_init, Y_init)-psi<<"\n";
    }

    //compute f for a given psi between psi0 and psi1
    double construct_f( double psi, double& R_0, double& Z_0)
    {
        find_initial( psi, R_0, Z_0);
        std::array<double,3> begin{ {0,0,0} }, end(begin), end_old(begin);
        begin[0] = R_0, begin[1] = Z_0;
        double eps = 1e10, eps_old = 2e10;
        unsigned N = 50;
        while( (eps < eps_old || eps > 1e-7)&& eps > 1e-14)
        {
            eps_old = eps, end_old = end; N*=2;
            dg::stepperRK( "Feagin-17-8-10",  fieldRZYT_, 0., begin, 2*M_PI, end, N);
            eps = sqrt( (end[0]-begin[0])*(end[0]-begin[0]) + (end[1]-begin[1])*(end[1]-begin[1]));
            if(m_verbose)std::cout << "\t error "<<eps<<" with "<<N<<" steps\t";
            //Attention: if this succeeds on the first attempt end_old won't be updated
        }
        if(m_verbose)std::cout << "\t error "<<eps<<" with "<<N<<" steps\t";
        if(m_verbose)std::cout <<end_old[2] << " "<<end[2] <<"\n";
        double f_psi = 2.*M_PI/end[2]; //this actually is 1/q the safety factor
        return f_psi;
    }

    double operator()( double psi)
    {
        double R_0, Z_0;
        return construct_f( psi, R_0, Z_0);
    }
    double f_prime( double psi)
    {
        //compute fprime
        double deltaPsi = fabs(psi)/100.;
        double fofpsi[4];<commit_msg>Make Fpsi more resilient to nan<commit_after>#pragma once

#include "dg/topology/grid.h"
#include "dg/topology/functions.h"
#include "dg/topology/interpolation.h"
#include "dg/topology/derivatives.h"
#include "dg/topology/geometry.h"
#include "dg/functors.h"
#include "dg/runge_kutta.h"
#include "dg/nullstelle.h"
#include "fluxfunctions.h"
#include "ribeiro.h"


namespace dg
{
namespace geo
{

///@cond
namespace flux
{
namespace detail
{

//This leightweights struct and its methods finds the initial R and Z values and the coresponding f(\psi) as
//good as it can, i.e. until machine precision is reached
//Note that f(psi) = 1/q(psi) (The safety factor)
struct Fpsi
{

    //firstline = 0 -> conformal, firstline = 1 -> equalarc
    Fpsi( const CylindricalFunctorsLvl2& psip, const CylindricalFunctorsLvl1& ipol, double x0, double y0, bool verbose = false):
        psip_(psip), fieldRZYT_(psip, ipol, x0, y0), fieldRZtau_(psip),m_verbose(verbose)
    {
        //Find O-point
        double R_O = x0, Z_O = y0;
        dg::geo::findOpoint( psip, R_O, Z_O);
        //define angle with respect to O-point
        fieldRZYT_ = dg::geo::flux::FieldRZYT(psip, ipol, R_O, Z_O);
        X_init = x0, Y_init = y0;
        while( fabs( psip.dfx()(X_init, Y_init)) <= 1e-10 && fabs( psip.dfy()( X_init, Y_init)) <= 1e-10)
            X_init +=  1.;
    }
    //finds the starting points for the integration in y direction
    void find_initial( double psi, double& R_0, double& Z_0)
    {
        unsigned N = 50;
        std::array<double, 2> begin2d{ {0,0} }, end2d(begin2d), end2d_old(begin2d);
        if(m_verbose)std::cout << "In init function\n";
        begin2d[0] = end2d[0] = end2d_old[0] = X_init;
        begin2d[1] = end2d[1] = end2d_old[1] = Y_init;
        double eps = 1e10, eps_old = 2e10;
        while( (eps < eps_old || eps > 1e-7) && eps > 1e-14)
        {
            eps_old = eps; end2d_old = end2d;
            N*=2; dg::stepperRK( "Feagin-17-8-10",  fieldRZtau_, psip_.f()(X_init, Y_init), begin2d, psi, end2d, N);
            eps = sqrt( (end2d[0]-end2d_old[0])*(end2d[0]-end2d_old[0]) + (end2d[1]-end2d_old[1])*(end2d[1]-end2d_old[1]));
        }
        X_init = R_0 = end2d_old[0], Y_init = Z_0 = end2d_old[1];
        if(m_verbose)std::cout << "In init function error: psi(R,Z)-psi0: "<<psip_.f()(X_init, Y_init)-psi<<"\n";
    }

    //compute f for a given psi between psi0 and psi1
    double construct_f( double psi, double& R_0, double& Z_0)
    {
        find_initial( psi, R_0, Z_0);
        std::array<double,3> begin{ {0,0,0} }, end(begin), end_old(begin);
        begin[0] = R_0, begin[1] = Z_0;
        double eps = 1e10, eps_old = 2e10;
        unsigned N = 50;
        unsigned nan_counter = 0;
        while( (eps < eps_old || eps > 1e-7)&& eps > 1e-14)
        {
            eps_old = eps, end_old = end; N*=2;
            dg::stepperRK( "Feagin-17-8-10",  fieldRZYT_, 0., begin, 2*M_PI, end, N);
            eps = sqrt( (end[0]-begin[0])*(end[0]-begin[0]) + (end[1]-begin[1])*(end[1]-begin[1]));
            if(m_verbose)std::cout << "\t error "<<eps<<" with "<<N<<" steps\t";
            if( std::isnan( eps) && nan_counter < 4) eps = 1e10, end = end_old, nan_counter++;
        }
        if(m_verbose)std::cout << "\t error "<<eps<<" with "<<N<<" steps\t";
        if(m_verbose)std::cout <<end_old[2] << " "<<end[2] <<"\n";
        double f_psi = 2.*M_PI/end[2]; //this actually is 1/q the safety factor
        return f_psi;
    }

    double operator()( double psi)
    {
        double R_0, Z_0;
        return construct_f( psi, R_0, Z_0);
    }
    double f_prime( double psi)
    {
        //compute fprime
        double deltaPsi = fabs(psi)/100.;
        double fofpsi[4];
{'fraction_non_alphanumeric': 0.1642538222162089, 'fraction_numerical': 0.03964280882153971, 'mean_word_length': 2.696, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4187802', 'n_tokens_mistral': 3114, 'n_tokens_neox': 2951, 'n_words': 691}