Dear meep-users,
To validate my simulation results of a dispersive medium I calculated the
propagation factor (and thus refractive index). It seems trivial but even using
the dispersive material given in meep tutorial (http://ab-
initio.mit.edu/wiki/index.php/Meep_Tutorial/Material_dispersion) I do not get
even close to the exact results.
What I do is, propagate a broadband pulse in a 1D grid filled with a dispersive
media and take Fourier transform of E-field at 2 distant points in the grid.
I compare the results with an analytical wave of E0 exp j(wt - kz). At the 2
points the wave would be,
E1 = E0 exp j(wt - k z1)
E2 = E0 exp j(wt - k z2)
E2/E1 = exp j k L ; L = z2-z1;
k = beta + j alpha = w/c * n;
E2/E1 = exp (- alpha L) * exp (j beta L).
Ideally I should be getting the same results from FDTD and analytical
propagation isn't it? Increasing resolution too does not converge.
Really appreciate your help. I have attached my .ctl file and matlab code.
Thanks
Philip
----------------------- ctl file -------------------------------
(set-param! resolution 100)
(set-param! Courant .5)
(define dpml 2.0)
(define lenx 6.0)
(define sx (+ (* 2 dpml) lenx))
(define-param fcen 1.0)
(define-param df 2.0)
(set! geometry-lattice (make lattice (size sx no-size no-size)))
(define dielectric-media (make dielectric (epsilon 2)))
(define dispersive-media (make dielectric (epsilon 2.25)
(E-polarizations
(make polarizability
(omega 1.1) (gamma 1e-5) (sigma 0.5))
(make polarizability
(omega 0.5) (gamma 0.1) (sigma 2e-5))
))
)
;(set! default-material dielectric-media)
(set! default-material dispersive-media)
(set! sources
(list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ez)
(center -2.5 0)
)
)
)
(set! pml-layers (list (make pml (thickness dpml) (direction X))))
(define dt (/ Courant resolution))
(run-sources+ 10
(stop-when-fields-decayed 10 Ez (vector3 2.5 0) 1e-3)
(to-appended "eps" (at-beginning output-epsilon))
(to-appended "ez" (at-every dt output-efield-z))
)
----------------------- ctl file -------------------------------
----------------------- matlab code ----------------------------
function hdf5_spectrum()
fname = 'D:\meep\test1-ez.h5';
hinfo = hdf5info(fname);
dset = hdf5read(hinfo.GroupHierarchy.Datasets(1));
resolution = 100;
Courant = 0.5;
p1 = floor(size(dset, 2)/2);
p2 = p1 + floor(100);
fcen = 1;
df = 1.99;
f = [fcen-df/2:df/1000:fcen+df/2];
dt = Courant/resolution;
dx = 1/resolution;
% E1 = E0 exp j(wt - k z1)
% E2 = E0 exp j(wt - k z2)
% E2/E1 = exp j k L ; L = z2-z1;
% k = beta + j alpha = w/c * n;
% E2/E1 = exp (- alpha L) * exp (j beta L)
E1 = calc_fft(f, dset(:, p1), dt);
E2 = calc_fft(f, dset(:, p2), dt);
L = (p2-p1) * dx;
E2E1 = E2./E1;
ampl = abs(E2E1);
alpha = -log(ampl)./L;
beta = (unwrap(angle(E2)) - unwrap(angle(E1)))/L;
w = 2*pi*f;
n1 = beta./w;
n11 = alpha./w;
n = n1 + j * n11;
er = n.^2;
figure(1);
plot(f, real(er));
title('Re[e_r]');
figure(2);
plot(f, imag(er));
title('Im[e_r]');
end
function y = calc_fft(f, x, dt)
y = zeros(1, numel(f));
for n=1:numel(x)
y = y + x(n) .* exp(j*2*pi*f*dt*n);
end
end
----------------------- matlab code ----------------------------
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
