     const container& polarisation( const std::vector<container>& y);
 
     container chi, omega;
+#if def APAR
+    container apar;
+#endif
     const container binv, curvR, curvZ, gradlnB;
     const container iris, source, damping;
 
     std::vector<container> phi, curvphi, dzphi, expy;
     std::vector<container> dzy, curvy; 
-
+#if def APAR
+    std::vector<container> arakAN,arakAN,arakAphi;
+#endif
     //matrices and solvers
     Matrix A; 
     dg::DZ<container> dz;
     dg::ArakawaX< dg::DMatrix, container>    arakawa; 
     //dg::Polarisation2dX< thrust::host_vector<value_type> > pol; //note the host vector
     dg::Polarisation< dg::DMatrix, container, container > pol; //note the host vector
+#if def APAR
+    dg::Maxwell<Matrix,container,container> maxwell;
+#endif
     dg::Invert<container>       invert_pol;
 
     const container w3d, v3d, one;

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

    const solovev::GeomParameters gp;
    const container w3d_, v3d_;
    container temp;
    const container lapiris_,pupil_;
    dg::DMatrix LaplacianM_perp;
//     dg::DMatrix LaplacianM_para;
};
