     double psipmax = dg::blas1::reduce( psipog2d, 0. ,thrust::maximum<double>()); //DEPENDS ON GRID RESOLUTION!!
     double volumeXGrid;
     /// -------  Elements for fsa of curvature operators ----------------
-    // This is a numerical test because the fsa of the curvature operators exactly vanishes (at least for I = cte)
+    // This is a numerical test because the fsa of the curvature operators exactly vanishes (at least for I = cte and almost for everything else)!!
     dg::geo::CylindricalVectorLvl0 bhat_ = dg::geo::createBHat( mag);
     dg::geo::CylindricalVectorLvl0 curvB_ = dg::geo::createTrueCurvatureNablaB( mag);
     dg::geo::CylindricalVectorLvl0 curvK_ = dg::geo::createTrueCurvatureKappa( mag);

        dg::blas1::scal( dvdzeta, 4.*M_PI*M_PI);
        dg::Grid1d gX1d( gX2d.x0(), gX2d.x1(), npsi, Npsi, dg::NEU);
        dg::HVec X_psi_vol = dg::integrate( dvdzeta, gX1d);
        map1d.emplace_back( "X_psi_vol", X_psi_vol,
            "Flux volume on X-point grid");
        dg::blas1::pointwiseDivide( X_psip1d, dvdzeta, X_psip1d);

        // h02 factor
        dg::HVec h02 = dg::evaluate( Hoo(mag), grid2d);
        fsa.set_container( h02);
        map1d.emplace_back( "hoo_fsa", dg::evaluate( fsa, grid1d),
           "Flux surface average of novel h02 factor with delta function");



        dg::geo::FluxVolumeIntegral<dg::HVec> fvi( (dg::CartesianGrid2d)grid2d, mag);
