         const dg::DVec v3d = dg::create::inv_volume( g3d);
 
         std::cout << "Computing ds..." << std::endl;
-        //DFA dsFA( solovev::ConformalField( gp, g3d.x(), g3d.f_x()), g3d, rk4eps, dg::NoLimiter(), dg::DIR); 
-        //dg::DS<DFA, dg::DMatrix, dg::DVec> ds( dsFA, solovev::ConformalField(gp, g3d.x(), g3d.f_x()), dg::not_normed, dg::centered, false);
-        DFA dsFA( solovev::OrthogonalField( gp, g2d, g2d.g()), g3d, rk4eps, dg::NoLimiter(), dg::DIR); 
-        dg::DS<DFA, dg::DMatrix, dg::DVec> ds( dsFA, solovev::OrthogonalField(gp, g2d, g2d.g()), dg::not_normed, dg::centered, false);
+        //DFA dsFA( conformal::Field( gp, g3d.x(), g3d.f_x()), g3d, rk4eps, dg::NoLimiter(), dg::DIR); 
+        //dg::DS<DFA, dg::DMatrix, dg::DVec> ds( dsFA, conformal::Field(gp, g3d.x(), g3d.f_x()), dg::not_normed, dg::centered, false);
+        DFA dsFA( orthogonal::Field( gp, g2d, g2d.f2_xy()), g3d, rk4eps, dg::NoLimiter(), dg::DIR); 
+        dg::DS<DFA, dg::DMatrix, dg::DVec> ds( dsFA, orthogonal::Field(gp, g2d, g2d.f2_xy()), dg::not_normed, dg::centered, false);
 
         std::cout << "ds constructed!" << std::endl;
         dg::DVec function = dg::pullback( funcNEU, g3d), derivative(function), temp(function),

        //Nzn = unsigned(Nz*pow(2,i));
        //Nxn = (unsigned)ceil(Nx*pow(2,(double)(i*2./n)));
        //Nyn = (unsigned)ceil(Ny*pow(2,(double)(i*2./n)));
        //solovev::ConformalRingGrid3d<dg::DVec> g3d(gp, psi_0, psi_1, n, Nxn, Nyn,Nzn, dg::DIR);
        //solovev::ConformalRingGrid2d<dg::DVec> g2d = g3d.perp_grid();
        solovev::OrthogonalRingGrid3d<dg::DVec> g3d(gp, psi_0, psi_1, n, Nxn, Nyn,Nzn, dg::DIR);
        solovev::OrthogonalRingGrid2d<dg::DVec> g2d = g3d.perp_grid();
        g3d.display();
        //g2d.display();
        std::cout << "NR = " << Nxn << std::endl;

        dg::DMatrix dx = dg::create::dx( g3d);
        dg::DMatrix dy = dg::create::dy( g3d);
        dg::DMatrix dz = dg::create::dz( g3d);
        dg::HVec hbx, hby;
        dg::blas1::transfer( function, hbx);
        dg::blas1::transfer( function, hby);
        dg::geo::pushForwardPerp( bhatR, bhatZ, hbx, hby, g3d);
        const dg::DVec bx(hbx), by(hby);
        dg::blas2::symv( dx, lnB, derivative);
