     }     
     double operator()( double R, double Z)
     {
-       return invB_(R,Z)/R/R*(bR_(R,Z) *psipZ_(R,Z) - bZ_(R,Z)* psipR_(R,Z)) ;
+       return invB_(R,Z)*invB_(R,Z)/R*(bR_(R,Z) *psipZ_(R,Z) - bZ_(R,Z)* psipR_(R,Z)) ;
     }
     double operator()( double R, double Z, double phi)
     {
-       return invB_(R,Z,phi)/R/R*(bR_(R,Z,phi) *psipZ_(R,Z,phi) - bZ_(R,Z,phi)* psipR_(R,Z,phi)) ;
+       return invB_(R,Z,phi)* invB_(R,Z,phi)/R*(bR_(R,Z,phi) *psipZ_(R,Z,phi) - bZ_(R,Z,phi)* psipR_(R,Z,phi)) ;
     }
     private:
 //     InvB invB_; 


#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
//6 analytical quantities
namespace dg

    double lnN_inner_;
    Psip psip_;
};
/**
 * @brief returns zonal flow field
 */ 
struct ZonalFlow
{
    ZonalFlow(GeomParameters gp, double k_psi):
        gp_(gp),
        k_psi_(k_psi),
        psip_(Psip(gp.R_0,gp.A,gp.c)) {
    }
    double operator() (double R, double Z) 
    {
      if (psip_(R,Z)<0.) return (1.+cos(2.*M_PI*psip_(R,Z)*k_psi_));
      return 1.;
      
    }
    double operator() (double R, double Z,double phi) 
    {
        if (psip_(R,Z,phi)<0.) return (1.+cos(2.*M_PI*psip_(R,Z,phi)*k_psi_));
        return 1.;
    }
    private:
    GeomParameters gp_;
    double k_psi_;
    Psip psip_;
};

/**
 * @brief testfunction to test the parallel derivative
 */ 
