I've been building a home-brew FDTD code to calculate band diagrams and I've been trying to get the diagram for a hexagonal lattice. Lumerical's website claimed if you have two sources that are a lattice vector apart out of phase by $i\textbf{k} \cdot \textbf{r}$ where k is the Bloch vector for the simulation and r is the lattice vector. However I've tried implementing this solution and still get band folding. Does anyone know how to do this or where there's more literature on the subject?

The dimensions of my simulation are 173x100 and my sources are given by:

Dz(15, 25) += ezInc(Time, 15.0, 25.0, Kx, Ky);

Dz(102, 75) += ezInc(Time, 102.0, 75.0, Kx, Ky);

Dz( 5, 13) += ezInc(Time, 5.0, 13.0, Kx, Ky);

Dz(92, 63) += ezInc(Time, 92.0, 63.0, Kx, Ky);

Dz(10, 10) += ezInc(Time, 10.0, 10.0, Kx, Ky);

Dz(97, 60) += ezInc(Time, 97.0, 60.0, Kx, Ky);

Dz(18, 30) += ezInc(Time, 18.0, 30.0, Kx, Ky);

Dz(105, 80) += ezInc(Time, 105.0, 80.0, Kx, Ky);

double ezInc(double time, double locX, double locY, double Kx, double Ky) {

double arg, dx, delta, omega, je, dt;

dx = 1e-8;

delta = (Kx*locX + Ky*locY)*dx;

// delta = 0;

je = 100;

dt = dx/(3e8*sqrt(2));

omega = 2.0*3e8*M_PI/(je*dx);

arg = (time*dt - dt*80)/16/dt;

arg = arg * arg;

return exp(-0.5*arg)*sin(omega*time*dt - delta);

}