     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
  */ 

  PsipR psipR_;
  PsipZ psipZ_;  
};
/**
 * @brief d B / d R
 */ 

//             yp[1][i] = -y[0][i]*psipR_(y[0][i],y[1][i],y[2][i])/ipol_(y[0][i],y[1][i],y[2][i]) ;             //dZ/dphi = -R/I Psip_Z
//             yp[2][i] =  y[0][i]*y[0][i]/invB_(y[0][i],y[1][i],y[2][i])/ipol_(y[0][i],y[1][i],y[2][i]);       //ds/dphi =  R^2 B/I
            yp[2][i] =  y[0][i]*y[0][i]/invB_(y[0][i],y[1][i])/ipol_(y[0][i],y[1][i]);       //ds/dphi =  R^2 B/I
            yp[0][i] = y[0][i]*psipZ_(y[0][i],y[1][i])/ipol_(y[0][i],y[1][i]);              //dR/dphi =  R/I Psip_Z
            yp[1][i] = -y[0][i]*psipR_(y[0][i],y[1][i])/ipol_(y[0][i],y[1][i]) ;             //dZ/dphi = -R/I Psip_Z
        }
    }
