     double z0 = 0, z1 = 2.*M_PI;
     std::cout << "Type n, Nx, Ny, Nz\n";    
     std::cin >> n >> Nxn >> Nyn >> Nzn;
-    
+
     dg::CylindricalGrid3d g3d( Rmin,Rmax, Zmin,Zmax, z0, z1,  n,Nxn ,Nyn, Nzn,dg::DIR, dg::DIR, dg::PER);
     dg::Grid2d g2d( Rmin,Rmax, Zmin,Zmax,  n, Nxn ,Nyn);
     


    const dg::DVec gradlnB_ = dg::evaluate(gradlnB, g3d);
  
    dsNU( function, derivative); //ds(f)
    dg::blas1::pointwiseDivide(ones,  inverseB, temp); //B
    ds.centeredDiv( temp, divBT); // dsT B
    ds.centeredDiv( ones, divbT);

    std::cout << "--------------------testing ds" << std::endl;
    double norm = dg::blas2::dot( w3d, solution);
    std::cout << "|| Solution ||   "<<sqrt( norm)<<"\n";
    double err =dg::blas2::dot( w3d, derivative);
    std::cout << "|| Derivative || "<<sqrt( err)<<"\n";
    dg::blas1::axpby( 1., solution, -1., derivative);
    err =dg::blas2::dot( w3d, derivative);
    std::cout << "Relative Difference in DS is "<< sqrt( err/norm )<<"\n"; 

    std::cout << "--------------------testing dsT" << std::endl;
    norm = dg::blas2::dot( w3d, divbsol);
    std::cout << "|| divbsol ||  "<<sqrt( norm)<<"\n";
    err = dg::blas2::dot( w3d, divbT);
    std::cout << "|| divbT  ||   "<<sqrt( err)<<"\n";
    dg::blas1::axpby( 1., divbsol, -1., divbT);
    err = dg::blas2::dot(divbT, w3d,divbT);
    std::cout << "Relative Difference in DST is   "<<sqrt(err/norm)<<"\n";
    std::cout << "-------------------- " << std::endl;
    double normdivBT =dg::blas2::dot(divBT, w3d,divBT);
    std::cout << "|| divB || "<<sqrt( normdivBT)<<"\n";  
    
    return 0;
}

    std::ifstream is("guenther_params.js");
    reader.parse(is,js,false);
    dg::geo::guenther::Parameters gp(js);
    //////////////////////////////////////////////////////////////////////////
    double Rmin=gp.R_0-1.0*gp.a;
    double Zmin=-1.0*gp.a*gp.elongation;
    double Rmax=gp.R_0+1.0*gp.a; 
    double Zmax=1.0*gp.a*gp.elongation;
    /////////////////////////////////////////////initialze fields /////////////////////    
    dg::geo::TokamakMagneticField mag = dg::geo::createGuentherField(gp.R_0, gp.I_0);
    dg::geo::InvB invb(mag);
    dg::geo::GradLnB gradlnB(mag);    
    dg::geo::Divb divb(mag);
    dg::geo::guenther::FuncNeu funcNEU(gp.R_0,gp.I_0);
    dg::geo::guenther::DeriNeu deriNEU(gp.R_0,gp.I_0);
    
    unsigned n=3;
    unsigned Nxn = 20;
