         container chi( x[0].size()), omega( x[0].size());
         for( unsigned i=0; i<2; i++)
             thrust::transform( x[i].begin(), x[i].end(), expy[i].begin(), dg::EXP<double>());
-        dg::blas1::pointwiseDot( expy[0], x[2], omega);
-        dg::blas1::pointwiseDot( expy[1], x[3], chi);
-        dg::blas1::axpby( -1., omega, 1., chi); //-N_eU_e + N_iU_i
+        dg::blas1::pointwiseDot( expy[0], x[2], omega); //N_e U_e 
+        dg::blas1::pointwiseDot( expy[1], x[3], chi); //N_i U_i
+        dg::blas1::axpby( -1., omega, 1., chi); //-N_e U_e + N_i U_i
         dg::blas1::pointwiseDivide( chi, expy[0], omega);//J_par/N_e
-        dg::blas1::pointwiseDivide( chi, expy[1], chi); //J_par/N_i
-
-        dg::blas1::axpby( -p.c/p.mu[0]/p.eps_hat, omega, 1., y[2]);
-        dg::blas1::axpby( -p.c/p.mu[1]/p.eps_hat, chi, 1., y[3]);
-        //cut contributions to boundary 
+        dg::blas1::pointwiseDivide( chi, expy[0], chi); //J_par/N_i    now //J_par/N_e  //n_e instead of n_i
+
+        dg::blas1::axpby( -p.c/p.mu[0]/p.eps_hat, omega, 1., y[2]);  // dtU_e =- C/hat(mu)_e J_par/N_e
+        dg::blas1::axpby( -p.c/p.mu[1]/p.eps_hat, chi, 1., y[3]);    // dtU_e =- C/hat(mu)_i J_par/N_i   //n_e instead of n_i
+        //cut contributions to boundary now with damping on all 4 quantities
         for( unsigned i=0; i<y.size(); i++)
             dg::blas1::pointwiseDot( pupil_, y[i], y[i]);
     }

        gp(gp),
        w3d_( dg::create::w3d(g)), v3d_(dg::create::v3d(g)),
        temp( g.size()),
        pupil_( dg::evaluate( solovev::Pupil( gp), g))
    {
        LaplacianM_perp = dg::create::laplacianM_perp( g, dg::normed, dg::symmetric);
        //LaplacianM_para = dg::create::laplacianM_parallel( g, dg::PER);
    }
    void operator()( const std::vector<container>& x, std::vector<container>& y)
    {

// 
//     dg::blas1::axpby( -p.c/p.mu[0]/p.eps_hat, omega, 1., yp[2]);
//     dg::blas1::axpby( -p.c/p.mu[1]/p.eps_hat, chi, 1., yp[3]);

    //add parallel diffusion with naive implementation
    for( unsigned i=0; i<4; i++)
    {
        dz(dzy[i], omega); //dz (dz (N,U))
        dg::blas1::axpby( p.nu_parallel, omega, 1., yp[i]);                     //dt(lnN,U) = dt(lnN,U) + dz (dz (lnN,U))
    }
    //add particle source to dtN
    for( unsigned i=0; i<2; i++)
    {
        dg::blas1::pointwiseDivide( source, expy[i], omega); //source/N
        dg::blas1::axpby( 1., omega, 1, yp[i]  );                               //dtlnN = dtlnN + source/N
    }
    //cut boundary terms 
    for( unsigned i=0; i<2; i++)
    {
//         dg::blas1::pointwiseDivide( damping, expy[i], omega); 
//         dg::blas1::axpby( -1., omega, 1., yp[i]);
//         dg::blas1::pointwiseDot( damping, y[i], omega); 
//         dg::blas1::axpby( -1., omega, 1., yp[i]);
    }
    for( unsigned i=0; i<4; i++) //damping and pupil
    {
        dg::blas1::pointwiseDot( damping, yp[i], yp[i]); 
//         dg::blas1::pointwiseDot( iris, yp[i], yp[i]);
    }
}

template< class container>
void Feltor< container>::curve( const container& src, container& target)
{
    dg::blas2::gemv( arakawa.dx(), src, target); //d_R src
    dg::blas2::gemv( arakawa.dy(), src, omega);  //d_Z src
    dg::blas1::pointwiseDot( curvR, target, target); // C^R d_R src
    dg::blas1::pointwiseDot( curvZ, omega, omega);   // C^Z d_Z src
    dg::blas1::axpby( 1., omega, 1., target ); // (C^R d_R + C^Z d_Z) src
    //Curvature_phi = b^phi
}

template< class container>
