 #include <iostream>
 #include <iomanip>
+#include <sstream>
+#include <fstream>
 #include <vector>
+#include "timer.cuh"
+#include <string>
+
 
 // #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 << "Type n, Nx, Ny, Nz\n";
    std::cin >> n>> Nx>>Ny>>Nz;
//     n=1;Nx=100;Ny=100;Nz=15;
//     Nxh=Nx/2;
    
//     Nyh=Ny/2;
//     Rmin=R_0-1;
//     Zmin=-Nyh;
//     Rmax=R_0+1; 
//     Zmax=Nyh;
    std::cout << "The grid parameters" <<"\n";
    std::cout  << Rmin<<"  " << Rmax <<"  " << Zmin <<"  " <<Zmax <<"\n";
//     std::cin >> n>> Nx>>Ny>>Nz;
    dg::Grid3d<double> g3d( Rmin,Rmax, Zmin,Zmax, 0, 2.*M_PI, n, Nx, Ny, Nz);
    const dg::DVec w3d = dg::create::w3d( g3d);
    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::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 i=0;i<5;i++)
    {
        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);
        const dg::DVec w3d = dg::create::w3d( g3d);


        std::cout <<"-----(1) test with testfunction" << "\n";
        dg::TestFunction func(psip);
        dg::DeriTestFunction derifunc(psip,psipR,psipZ,ipol,invB);
        std::cout << "Construct parallel  derivative\n";
        dg::Timer t;
        t.tic();
        dg::DZ<dg::DVec> dz( field, g3d); 
        t.toc();
        std::cout << "Creation of parallel Derivative took "<<t.diff()<<"s\n";

    //     std::cout << "Evaluate functions on the grid\n";
        dg::DVec function = dg::evaluate( func, g3d),dzfunc(function);
        dg::DVec diff(g3d.size());

        const dg::DVec solution = dg::evaluate( derifunc, g3d);
        std::cout << "compute parallel z derivative\n";
        t.tic();
        dz( function, dzfunc);
        t.toc();
        std::cout << "Compuation of parallel Derivative took "<<t.diff()<<"s\n";    
        //Compute norm of computed dz
        double normdz = dg::blas2::dot( w3d, dzfunc);
        std::cout << "Norm dz "<<sqrt( normdz)<<"\n";
        //Compute norm of analytical solution
        double normsol = dg::blas2::dot( w3d,solution);
        std::cout << "Norm solution "<<sqrt( normsol)<<"\n";
        dg::blas1::axpby( 1., solution, -1., dzfunc,diff);
        double normdiff = dg::blas2::dot( w3d, diff);
        std::cout << "Norm diff "<<sqrt( normdiff)<<"\n";
        double reldiff=sqrt( dg::blas2::dot( diff, w3d, diff)/normsol );
        std::cout << "Relative Difference Is "<< reldiff<<"\n";

        std::cout <<"-----(2) test with gradlnb" << "\n";    
        //Evaluate analyitcal gradlnb on grid
        dg::DVec gradLnBsolution = dg::evaluate( gradLnB, g3d);
        //compute Ln(B) on grid
        dg::DVec lnBongrid = dg::evaluate( lnB, g3d);
        dg::DVec dzlnBongrid(g3d.size());
        dg::DVec diff2(g3d.size());
        std::cout << "compute parallel z derivative\n";
        t.tic();
        dz(lnBongrid,dzlnBongrid);
        t.toc();
        std::cout << "Compuation of parallel Derivative took "<<t.diff()<<"s\n";    
        //Compute norm of computed dz
        double normdz2 = dg::blas2::dot( w3d, dzlnBongrid);
        std::cout << "Norm dz "<<sqrt( normdz2)<<"\n";
        //Compute norm of analytical gradlnb solution
        double normsol2 = dg::blas2::dot( w3d,gradLnBsolution);
        std::cout << "Norm solution "<<sqrt( normsol2)<<"\n";
        //Compute difference between the solutions
        dg::blas1::axpby( 1., gradLnBsolution , -1., dzlnBongrid,diff2); //diff = gradlnB - dz(ln(B))
        double normdiff2=dg::blas2::dot( w3d, diff2);
        std::cout << "Norm diff "<<sqrt( normdiff2)<<"\n"; //=sqrt(gradlnB - dz(ln(B)))
        double reldiff2 =sqrt( dg::blas2::dot( diff2, w3d, diff2)/normsol2 );
        std::cout << "Relative Difference Is "<<reldiff2 <<"\n";
        dzerrfile1 << reldiff  << std::endl;
        dzerrfile2 << reldiff2 << std::endl;
    }
    dzerrfile1.close();
    dzerrfile2.close();
    }

    return 0;
}
