 //dg::Timer t;
 //t.tic();
 //t.toc();
-// boost::timer t;
-// double duration;
-//
-// t.restart();
-// duration = t.elapsed();
 
 struct Heaviside2d
 { Heaviside2d( double sigma):sigma2_(sigma*sigma), x_(0), y_(0) {}

      err_out = nc_put_vara_double(species_wgrp_id[j], species_wspatial_id[7 + add_to_spatial_idx], start0d, count0d, &posY_max);
      velX_max = (posX_max - posX_max_old[j]) /deltaT;
      velY_max = (posY_max - posY_max_old[j]) /deltaT;
      if (i>0)
      { posX_max_old[j] = posX_max;
        posY_max_old[j] = posY_max;
      }
      err_out = nc_put_vara_double(species_wgrp_id[j], species_wspatial_id[8 + add_to_spatial_idx], start0d, count0d, &velX_max);
      err_out = nc_put_vara_double(species_wgrp_id[j], species_wspatial_id[9 + add_to_spatial_idx], start0d, count0d, &velY_max);
      //compactness
      if (i==0)
      { heavi.set_origin( posX_max_hs, posY_max_hs );
        heavy = dg::evaluate( heavi, g2d);
        normalize = dg::blas2::dot( heavy, w2d, ntilde[j]);

    { err_out = nc_put_vara_double(err_time_wgrp_id, err_time_wval_id[j], start0d, count0d, &energy[j]);
    }
    double dcn = cn[0]+cn[1]+cn[2];
    dcn = fabs(dcn/cn[0]);
    std::cout << "Rel. error charge conservation " << dcn << "\n";
    //double dvn = +p.a[0]*mass_[0]*velX[0]+p.a[1]*mass_[1]*velX[1]+p.a[2]*mass_[2]*velX[2];
    //std::cout << "Abs. error flux   conservation " << velX[0]<< " "<<velX[1]<<" "<<velX[2]<<" "<<fabs(dvn) << "\n";
    double dvn = +p.a[0]*NvelX[0]+p.a[1]*NvelX[1]+p.a[2]*NvelX[2];
    dvn = fabs(dvn/(p.a[0]*NvelX[0]));
    std::cout << "Rel. error flux   conservation " << dvn << "\n";
    err_out = nc_put_vara_double(err_time_wgrp_id, err_time_wval_id[5], start0d, count0d, &dcn);
    err_out = nc_put_vara_double(err_time_wgrp_id, err_time_wval_id[5], start0d, count0d, &dvn);
  }
  err_out = nc_close(ncid_out);
  err_in = nc_close(ncid_in);
