to attract a considerable attention of the scientific community due to their unique properties and potential technological applications [@sak]. To a large extent, PhCs applications are based on the photonic bandgap effect and involve sophisticated defect structures, cavities and wave guides in the periodic crystal host. The full dynamics of the electromagnetic field in such structures may be studied in principle rigorously by direct numerical solution of Maxwell’s equations using the finite difference time domain (FDTD) method [@Taflove2000]. For the calculation of stationary modes, however, the numerical effort may be substantially reduced by using the Galerkin method [@Fletcher1984], i.e., by expanding the electromagnetic field in terms of an appropriate orthonormal set of basis functions which renders the stationary electromagnetic wave equation as a discrete matrix eigenvalue problem. In this case the proper choice of the set of basis functions is crucial in order to obtain an accurate description while keeping the dimension of the eigenvalue problem minimal.

While for extended wave problems an expansion in terms of Bloch functions, the eigenmodes of the unperturbed PhC, is natural [@Chigrin2009; @Kremers2009a; @Kremers2009], for the description of defect structures the use of Wannier functions as a basis set [@LEUNG1993; @Lidorikis98; @Albert2000; @Albert2002] is 