     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);
     

    
    std::cout << "computing dsDIR" << std::endl;
    dg::geo::Fieldaligned<dg::aProductGeometry3d, dg::IDMatrix, dg::DVec>  dsFA( mag, g3d, dg::DIR, dg::DIR, dg::geo::FullLimiter(), rk4eps, mx, my, true, true, integrate);
    std::cout << "computing dsNEU" << std::endl;
    dg::geo::Fieldaligned<dg::aProductGeometry3d, dg::IDMatrix, dg::DVec> dsNUFA( mag, g3d,dg::NEU, dg::NEU, dg::geo::FullLimiter(), rk4eps, mx, my, true, true, integrate);

    dg::geo::DS<dg::aProductGeometry3d, dg::IDMatrix, dg::DMatrix, dg::DVec> ds ( dsFA, dg::not_normed, dg::centered), 
        dsNU ( dsNUFA, dg::not_normed, dg::centered);
    
    dg::DVec function = dg::evaluate( funcNEU, g3d) ,
                        temp( function),

    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;
