 
     std::vector<double> c(13);
     //read and store geom data
-    std::vector<double> v;
-    try{ v = file::read_input( "geometry_params.txt"); }
+       std::vector<double> v, v2;
+    try{ 
+        if( argc==1)
+        {
+            v = file::read_input( "geometry_params.txt"); 
+            v2 = file::read_input( "../feltor/input.txt");
+
+        }
+        else
+        {
+            v = file::read_input( argv[1]); 
+            v2 = file::read_input( argv[2]);
+        }
+    }
     catch (toefl::Message& m) {  
         m.display(); 
         for( unsigned i = 0; i<v.size(); i++)
             std::cout << v[i] << " ";
             std::cout << std::endl;
         return -1;}
-    for( unsigned i = 1; i<v.size(); i++)
-    std::stringstream title;
+
     //write parameters from file into variables
     const solovev::GeomParameters gp(v);
+    const eule::Parameters p(v2);
+    p.display( std::cout);
     gp.display( std::cout);
 
-    Rmin=gp.R_0-gp.a;
-    Zmin=-gp.a*gp.elongation;
-    Rmax=gp.R_0+gp.a; 
-    Zmax=gp.a*gp.elongation;
+    Rmin=gp.R_0-(p.boxscale)*gp.a;
+    Zmin=-(p.boxscale)*gp.a*gp.elongation;
+    Rmax=gp.R_0+(p.boxscale)*gp.a; 
+    Zmax=(p.boxscale)*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  << 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;
-    
-    
+        
     solovev::Field field(gp);
     solovev::Psip psip(gp.R_0,gp.A,gp.c);
     solovev::PsipR psipR(gp.R_0,gp.A,gp.c);

    double R_0;
}; 

int main()
{
    double Rmin,Zmin,Rmax,Zmax;
//     double A,R_0,a, elongation;

                std::cout << "Rel Diff = "<<reldiff3 <<"\n";
                dzerrfile1 << pow(2,zz)*Nz <<" " << reldiff << std::endl;
                dzerrfile2 << pow(2,zz)*Nz <<" " << reldiff2 << std::endl;
                
                std::cout <<"----(4) test div(B) != 0 "<<"\n";
                dg::DVec bRongrid = dg::evaluate( fieldR, grid);
                dg::DVec bZongrid = dg::evaluate( fieldZ, grid);
                dg::DVec dRbR(g3d.size());
                dg::DVec dZbZ(g3d.size());                
                dg::DVec invRbR(g3d.size());                
//                    //cut boundaries
//                 dg::blas1::pointwiseDot( pupilongrid, bRongrid,bRongrid); 
//                 dg::blas1::pointwiseDot( pupilongrid, bZongrid, bZongrid); 
                
                dg::DVec divB(g3d.size());                
                dg::blas2::gemv( arakawa.dx(), bRongrid, dRbR);
                dg::blas2::gemv( arakawa.dy(), bZongrid, dZbZ);
                dg::blas1::pointwiseDot( invnormrongrid , bRongrid, invRbR);
                dg::blas1::axpby( 1., dRbR   , 1., dZbZ, divB);
                dg::blas1::axpby( 1./gp.R_0, invRbR , 1., divB);
                double normdivB2= dg::blas2::dot( w3d, divB); 
                std::cout << "divB = "<<sqrt( normdivB2)<<"\n";
                
             }
            dzerrfile1.close();
            dzerrfile2.close();
