         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]);
     }

    curvZ( dg::evaluate(solovev::CurvatureZ(gp), g)),
    gradlnB( dg::evaluate(solovev::GradLnB(gp) , g)),
    iris( dg::evaluate( solovev::Pupil( gp), g)),
    source( dg::evaluate( dg::Gaussian( p.R_0, 0, p.b, p.b, p.amp_source, 0), g)),
    damping( dg::evaluate( solovev::Damping(gp ), g)), 
    phi( 2, chi), curvphi( phi), dzphi(phi), expy(phi),  
    dzy( 4, chi), curvy(dzy),
    A (dg::create::laplacianM_perp( g, dg::not_normed, dg::symmetric)),
    dz(solovev::Field(gp), g),
    arakawa( g), 
    pol(     g), 
    invert_pol( omega, omega.size(), p.eps_pol), 

        {
            dg::blas2::gemv( LaplacianM_perp, x[i], temp);
            dg::blas2::gemv( LaplacianM_perp, temp, y[i]);
            dg::blas1::axpby( -p.nu_perp, y[i], 0., y[i]); // - nu_perp lapl_RZ (lapl_RZ (lnN,U)) //factor MISSING!?!
//             dg::blas2::gemv( LaplacianM_para, x[i], temp);
//             dg::blas1::axpby(  -p.nu_parallel, temp, 1., y[i]);
        }
        //add parallel resistivity
        std::vector<container>  expy(2);
