 #include "interpolation.cuh"
 #include "draw/host_window.h"
 
-
-
-
-
-
 int main()
 {
     double Rmin,Zmin,Rmax,Zmax;
 //     double A,R_0,a, elongation;
 //     double psipmin,psipmax;
     unsigned n, Nx, Ny, Nz;
-    
 
     std::vector<double> c(13);
     //read and store geom data

#include <iostream>
#include <iomanip>
#include <vector>

// #include "typedefs.cuh"
#include "xspacelib.cuh"

    Zmin=-1.1*gp.a*gp.elongation;
    Rmax=gp.R_0+1.1*gp.a; 
    Zmax=1.1*gp.a*gp.elongation;
    std::cout << "The grid parameters" <<"\n";
    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;
    for (unsigned k=1;k<4;k++)
    {
    std::stringstream ss1,ss2;
    ss1 << "dzerr1n" <<k*n<<"Nz"<<Nz<<".txt";
    ss2 << "dzerr2n" <<k*n<<"Nz"<<Nz<<".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());
    dg::Field field(gp);
    dg::Psip psip(gp.R_0,gp.A,gp.c);
    dg::PsipR psipR(gp.R_0,gp.A,gp.c);
