     std::cout  << Rmin<<"rho_s " << Rmax <<"rho_s " << Zmin <<"rho_s " <<Zmax <<"rho_s " <<"\n";
     std::cout << "Type n, Nx, Ny, Nz\n";
     std::cin >> n>> Nx>>Ny>>Nz;
-        dg::Field field(gp);
-    dg::Psip psip(gp.R_0,gp.A,gp.c);
-    dg::PsipR psipR(gp.R_0,gp.A,gp.c);
-    dg::PsipRR psipRR(gp.R_0,gp.A,gp.c);  
-    dg::PsipZ psipZ(gp.R_0,gp.A,gp.c);  
-    dg::PsipZZ psipZZ(gp.R_0,gp.A,gp.c);   
-    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(gp);
-    dg::CurvatureZ curvatureZ(gp);
-    dg::GradLnB gradLnB(gp);
-  
-    for (unsigned k=1;k<3;k++) //n iterator
+    
+    
+    solovev::Field field(gp);
+    solovev::Psip psip(gp.R_0,gp.A,gp.c);
+    solovev::PsipR psipR(gp.R_0,gp.A,gp.c);
+    solovev::PsipRR psipRR(gp.R_0,gp.A,gp.c);  
+    solovev::PsipZ psipZ(gp.R_0,gp.A,gp.c);  
+    solovev::PsipZZ psipZZ(gp.R_0,gp.A,gp.c);   
+    solovev::PsipRZ psipRZ(gp.R_0,gp.A,gp.c);  
+    solovev::Ipol ipol(gp.R_0,gp.A,psip);
+    solovev::InvB invB(gp.R_0,ipol,psipR,psipZ);
+    solovev::LnB lnB(gp.R_0,ipol,psipR,psipZ);
+    solovev::BR bR(gp.R_0,gp.A,psipR,psipRR,psipZ,psipRZ,invB);
+    solovev::BZ bZ(gp.R_0,gp.A,psipR,psipZ,psipZZ,psipRZ,invB);
+    solovev::CurvatureR curvatureR(gp);
+    solovev::CurvatureZ curvatureZ(gp);
+    solovev::GradLnB gradLnB(gp);
+    solovev::Pupil pupil(gp);
+    
+    dg::HVec v5(1, 0);
+    std::vector<thrust::host_vector<double> > in(3, v5);
+    std::vector<thrust::host_vector<double> > out(3, v5);
+    in[0][0]=gp.R_0+0.9*gp.a; 
+    in[1][0]=0.9*gp.a*gp.elongation;
+    in[2][0]=0.;
+
+    dg::integrateRK4( field, in, out,  2*M_PI, gp.rk4eps);
+    
+    std::cout << "Check single field by integrating from 0 to 2pi" << "\n";
+    std::cout <<"Rin =  "<< in[0][0] <<" Zin =  "<<in[1][0] <<" sin  = "<<in[2][0]<<"\n";
+    std::cout <<"Rout = "<< out[0][0]<<" Zout = "<<out[1][0]<<" sout = "<<out[2][0]<<"\n";
+
+    for (unsigned k=1;k<2;k++) //n iterator
     {
-        for (unsigned i=0;i<6;i++) //Nxy iterator
+        for (unsigned i=0;i<1;i++) //Nxy iterator
         {
             std::stringstream ss1,ss2;
             ss1 << "dzerr1n" <<k*n<<"Nxy"<<pow(2,i)* Nx<<".txt";

    for( unsigned i = 1; i<v.size(); i++)
    std::stringstream title;
    //write parameters from file into variables
    const dg::GeomParameters gp(v);
    gp.display( std::cout);

    Rmin=gp.R_0-1.1*gp.a;

#include "file/read_input.h"
#include "dg/backend/timer.cuh"
#include "dg/backend/xspacelib.cuh"
#include "geometry.h"
#include "dg/backend/dz.cuh"
#include "dg/algorithm.h"
#include "dg/backend/functions.h"
