     dg::CurvatureR curvatureR(gp);
     dg::CurvatureZ curvatureZ(gp);
     dg::GradLnB gradLnB(gp);
-    for (unsigned i=0;i<5;i++)
+  
+    for (unsigned k=1;k<3;k++) //n iterator
     {
-        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);
+    for (unsigned i=0;i<6;i++) //Nxy iterator
+    {
+
+
+    std::stringstream ss1,ss2;
+    ss1 << "dzerr1n" <<k*n<<"Nxy"<<pow(2,i)* Nx<<".txt";
+    ss2 << "dzerr2n" <<k*n<<"Nxy"<<pow(2,i)* Nx<<".txt";
+    std::string dzerr1fn = ss1.str();
+    std::string dzerr2fn = ss2.str();
+//     std::cout << dzerr1fn;
+    std::ofstream dzerrfile1((char *) dzerr1fn.c_str());
+    std::ofstream dzerrfile2((char *) dzerr2fn.c_str());
+    for (unsigned zz=0;zz<5;zz++) //Nz iterator
+    {
+        std::cout << "n = " << k*n << " Nx = " <<pow(2,i)* Nx << " Ny = " <<pow(2,i)* Ny << " Nz = "<<pow(2,zz)* Nz <<"\n";
+        dg::Grid3d<double> g3d( Rmin,Rmax, Zmin,Zmax, 0, 2.*M_PI,k*n,pow(2,i)* Nx,pow(2,i)* Ny, pow(2,zz)*Nz);
         const dg::DVec w3d = dg::create::w3d( g3d);
 
 

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

int main()
{

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