dispersion diagrams
My research interest is Photonic Crystals and periodic dispersive structures. As you know, dispersion diagram is one of the most important properties of periodic structures which is noticed extremely.
To calculate dispersion diagrams by FDTD, many solutions in papers are considered. one of the popular methods is FDTD to obtain dispersion diagrams in periodic structures, finite or infinite. Dispersion diagram defines relationship between wave number and frequency.
Could you help me to understand how dispersion diagrams of a periodic structure can be calulated? of course, the calculation methods(PWM,FDTD,FEM,...) are not important in this part.
I have soem experience to simulate periodic structures by FDTD. But,I have encountered some problems to obtain dispersion diagram. That's why I dont know exactly how I can relate fourier transform results to wave number while I have got samples in different points for many time steps.
Any suggestion would be helpful.
This is a good tutorial on what a dispersion diagram is and how to generate it:
http://www.emtalk.com/tut_2.htm
I also have such problem.Does anyone konw how to calculate dispersion diagrams by FDTD in matlab.
Hi, essentially a dispersion diagram is a plot of the wavelength in the dispersive material, lambda, against the wavelength in free space, lambda0.
In a periodic system it becomes simply a plot of phase (0->pi) vs. frequency.
Outside of this range the system is degenerate
omega(phase)=omega(phase+N*pi), for N=±integer.
Conventionally (in photonic crystal literature), this is plotted as a graph with wavevector as the abscissa and frequency as the ordinate.
The frequency, omega, of the wave is constant in all media.
The wavevector, k, changes as k=omega*n/c, where c is the speed of light in vacuum and n is the phase index of the wave.
The phase velocity v_p=omega/k=c/n is the gradient of a line passing through the origin, intersecting the point (k,omega).
The local slope of the dispersion diagram yields the group velocity, v_g=(omega)/dk.
Often, due to the scalability of Maxwell's equations, these units are normalised.
omega -> u (normalised frequency dimensionless), u=omega*a/2/pi/c=a/lambda0. Where a is period.
k -> is expressed in units of pi/a. k=2*pi*n/lambda0=2*pi/lambda
For the simplest implementations of these techniques:
In FDTD and PWM the wavevector is the independent variable and the frequency the eigenvalue.
In TMM, EME and FEM the frequency is the independent variable and wavevector is the eigenvalue.
The periodic boundary condition is extremely simple to implement in a 1D FDTD system and follows:
P. Harms, R. Mittra and W. Ko, "Implementation of the periodic boundary condition in the finite-difference time-domain algorithm for FSS structures,"
IEEE Trans. on Antennas and Propagation, vol. 42, no. 9, pp. 1317-1324, September 1994.
GRID STRUCTURE
h1...........................hN+1
|_|_|_|_|_|_|_|_|_|_|
..|_|_|_|_|_|_|_|_|_|
.e1.........................eN
..|---..period=a... ---|
..0........................(N-1)dx
The first and last e values sit on the periodic boundary (at real space positions z=0, and z=(N-1)dx) and the first and last h values (sit at z=-dx/2 and z=(N-1)dx+dx/2) outside the period.
ALGORITHM
Firstly we deal with the case of k=0.
The fields at periodic intervals are identical.
1: e1:eN are calculated from diff(h1:hN+1), and in turn
2: h2:hN-1 are calculated from diff(e1:eN),
3: the values outside the period are then set hN+1=h2 and h1=hN
this is equivalent to a loop, i.e. joining the two ends of the simulation together.
Secondly we deal with the case of k=pi/a.
The fields at periodic intervals are opposites. The Periodic boundaries are now updated using,
3: hN+1=-h2 and h1=-hN
The fields in the periodic system, at k=0 are perfectly in phase with the lattice.
The fields at k=pi/a are pi out of phase with the lattice.
These two separate k points can be calculated using a single real field.
For values of k between 0 and pi, we require two fields which are pi out of phase in order to correctly apply the phase shift at the boundary.
This is obviously twice as computationally expensive as for the two high symmetry points.
e and h are now complex variables (for ease of storage, and convenience).
The Periodic boundaries are now updated using,
3: hN+1= h2 * exp(j*k*a), h1 = hN * exp(-j*k*a)
These reduce to single field form at k=0 and pi, (their computation is unnecessarily expensive in 2D and 3D).
I've listed this as Matlab code, (I've tested it also in Octave) as the notation is concise.
We divide the period into N Ex and N Hy components, as the HN+1 form above is a little overzealous.
We study an infinitesimally small periodic index perturbation, ie. fill the space with air (n=1).
The resulting folding is an artefact of the model.
Note:
h1.......................hN
|_|_|_|_|_|_|_|_|_|_
..|_|_|_|_|_|_|_|_|_|
.e1........................eN
The fields are stored at various points as a function of time, for each k, and then FFTing them gives us the frequencies.
For the source choose a Gaussian modulated cosine.
Launch and store points can preferentially excite/detect different modes with varying efficiency.
So to plot the dispersion diagram using FDTD,
a: Set phase (in discrete steps, 0->pi) <--?
b: iterate steps 1,2,3 for sufficient time __|
c: FFT stored field values and plot.
If a particular resonance occurs at say frequency ii in the fft spectrum.
Its normalised frequency is u=a*[ii-1]./(dt*(timemax-1)).
Range u=a*[0:timemax-1]./(dt*(timemax-1))
Resolution is limited by length of simulation.
A more advanced method for determining the frequency eigenvalues is the filter diagonalisation method (FDM).
A free (linux based) FDTD software package and FDM program are available at:
http://ab-initio.mit.edu/wiki/index.php/Meep
http://ab-initio.mit.edu/wiki/index.php/Harminv
Hope this helps!
Hey,
u can use MoM method, attached is a paper regarding that
Yash