     dg::PsipRZ psipRZ(gp.R_0,gp.A,gp.c);  
     dg::Ipol ipol(gp.R_0,gp.A,psip);
     dg::InvB invB(gp.R_0,ipol,psipR,psipZ);
+    dg::LnB lnB(gp.R_0,ipol,psipR,psipZ);
     dg::BR bR(gp.R_0,gp.A,psipR,psipRR,psipZ,psipRZ,invB);
     dg::BZ bZ(gp.R_0,gp.A,psipR,psipZ,psipZZ,psipRZ,invB);
-    dg::CurvatureR curvatureR(invB,bZ);
-    dg::CurvatureZ curvatureZ(invB,bR);
-    dg::GradLnB gradLnB(invB, psipR,psipZ,bR,bZ);
+    dg::CurvatureR curvatureR(gp);
+    dg::CurvatureZ curvatureZ(gp);
+    dg::GradLnB gradLnB(gp);
+    for (unsigned i=0;i<5;i++)
+    {
+        std::cout << "n = " << k*n << " Nx = " <<pow(2,i)* Nx << " Ny = " <<pow(2,i)* Ny << " Nz = "<< Nz <<"\n";
+        dg::Grid3d<double> g3d( Rmin,Rmax, Zmin,Zmax, 0, 2.*M_PI,k*n,pow(2,i)* Nx,pow(2,i)* Ny, Nz);
+        const dg::DVec w3d = dg::create::w3d( g3d);
 
-    dg::TestFunction func(psip);
-    dg::DeriTestFunction deri(psip,psipR,psipZ,ipol,invB);
-    std::cout << "Construct z derivative\n";
-    dg::DZ<dg::DVec> dz( field, g3d); 
-    std::cout << "Evaluate functions on the grid\n";
-    dg::DVec function = dg::evaluate( func, g3d),derivative(function);
-   
-    const dg::DVec solution = dg::evaluate( deri, g3d);
 
-    dz( function, derivative);
-    dg::blas1::axpby( 1., solution, -1., derivative);
-    double norm = dg::blas2::dot( w3d, solution);
-    std::cout << "Norm Solution "<<sqrt( norm)<<"\n";
-//     double norm2 = dg::blas2::dot( w3d, solution);
-    std::cout << "Relative Difference Is "<< sqrt( dg::blas2::dot( derivative, w3d, derivative)/norm )<<"\n";    
+        std::cout <<"-----(1) test with testfunction" << "\n";
+        dg::TestFunction func(psip);
+        dg::DeriTestFunction derifunc(psip,psipR,psipZ,ipol,invB);
+        std::cout << "Construct parallel  derivative\n";
+        dg::Timer t;
+        t.tic();
+        dg::DZ<dg::DVec> dz( field, g3d); 
+        t.toc();
+        std::cout << "Creation of parallel Derivative took "<<t.diff()<<"s\n";
+
+    //     std::cout << "Evaluate functions on the grid\n";
+        dg::DVec function = dg::evaluate( func, g3d),dzfunc(function);
+        dg::DVec diff(g3d.size());
+
+        const dg::DVec solution = dg::evaluate( derifunc, g3d);
+        std::cout << "compute parallel z derivative\n";
+        t.tic();
+        dz( function, dzfunc);
+        t.toc();
+        std::cout << "Compuation of parallel Derivative took "<<t.diff()<<"s\n";    
+        //Compute norm of computed dz
+        double normdz = dg::blas2::dot( w3d, dzfunc);
+        std::cout << "Norm dz "<<sqrt( normdz)<<"\n";
+        //Compute norm of analytical solution
+        double normsol = dg::blas2::dot( w3d,solution);
+        std::cout << "Norm solution "<<sqrt( normsol)<<"\n";
+        dg::blas1::axpby( 1., solution, -1., dzfunc,diff);
+        double normdiff = dg::blas2::dot( w3d, diff);
+        std::cout << "Norm diff "<<sqrt( normdiff)<<"\n";
+        double reldiff=sqrt( dg::blas2::dot( diff, w3d, diff)/normsol );
+        std::cout << "Relative Difference Is "<< reldiff<<"\n";
+
+        std::cout <<"-----(2) test with gradlnb" << "\n";    
+        //Evaluate analyitcal gradlnb on grid
+        dg::DVec gradLnBsolution = dg::evaluate( gradLnB, g3d);
+        //compute Ln(B) on grid
+        dg::DVec lnBongrid = dg::evaluate( lnB, g3d);
+        dg::DVec dzlnBongrid(g3d.size());
+        dg::DVec diff2(g3d.size());
+        std::cout << "compute parallel z derivative\n";
+        t.tic();
+        dz(lnBongrid,dzlnBongrid);
+        t.toc();
+        std::cout << "Compuation of parallel Derivative took "<<t.diff()<<"s\n";    
+        //Compute norm of computed dz
+        double normdz2 = dg::blas2::dot( w3d, dzlnBongrid);
+        std::cout << "Norm dz "<<sqrt( normdz2)<<"\n";
+        //Compute norm of analytical gradlnb solution
+        double normsol2 = dg::blas2::dot( w3d,gradLnBsolution);
+        std::cout << "Norm solution "<<sqrt( normsol2)<<"\n";
+        //Compute difference between the solutions
+        dg::blas1::axpby( 1., gradLnBsolution , -1., dzlnBongrid,diff2); //diff = gradlnB - dz(ln(B))
+        double normdiff2=dg::blas2::dot( w3d, diff2);
+        std::cout << "Norm diff "<<sqrt( normdiff2)<<"\n"; //=sqrt(gradlnB - dz(ln(B)))
+        double reldiff2 =sqrt( dg::blas2::dot( diff2, w3d, diff2)/normsol2 );
+        std::cout << "Relative Difference Is "<<reldiff2 <<"\n";
+        dzerrfile1 << reldiff  << std::endl;
+        dzerrfile2 << reldiff2 << std::endl;
+    }
+    dzerrfile1.close();
+    dzerrfile2.close();
+    }
+
     return 0;
 }

#include "interpolation.cuh"
#include "draw/host_window.h"






int main()
{
    double Rmin,Zmin,Rmax,Zmax;
//     double A,R_0,a, elongation;
//     double psipmin,psipmax;
    unsigned n, Nx, Ny, Nz;
    

    std::vector<double> c(13);
    //read and store geom data

#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <vector>
#include "timer.cuh"
#include <string>


// #include "typedefs.cuh"
#include "xspacelib.cuh"
