On Thu, 19 Apr 2007, secondlw wrote:
> Dear Steven and others, How can I get the reflection phase of a
> structure if I use a monochromatic plane wave ? If I set meep to
> calculate with complex field true will help me? Thanks in advance! BR
> Huang
First, your question is not necessarily well-defined.
You only have a well-defined phase if you can isolate the reflection of a
single propagating mode. For example, this is true for structures that
are uniform in the directions transverse to the incident direction (e.g.
multilayer films), or more generally for periodic structures with
periodicity shorter than the diffraction limit so that there are no
non-specular reflections. It is also true if you have a single-mode
waveguide, in which case if you look far enough away you will see just the
one guided mode in the waveguide.
I will assume that this is true in the remainder of the email. Otherwise,
you have to define more precisely what you mean by "the" phase.
In this case, it is quite straightforward to get the reflected phase, and
even the reflected phase as a function of frequency, because when you have
a single mode the complex reflection coefficient can be determined by the
field at a single point (*any* point where the field is non-zero and you
have a single mode).
Just do something like the following:
1) As your source, input a short Gaussian pulse.
2) Pick a point where you want to measure the reflected field. This
should be a point *outside* the PML, and *far* enough from the structure
so that any evanescent fields etc. have died away. Also, put your point
in *front* of the source so that you can also use it to measure
transmission in your normalization run, below. Also, pick a field
component to look at -- it can be any field component that is non-zero.
Define a couple of variables to record this information so you can refer
to it later and easily change it
(define-param R-comp Ez) ; field component to measure reflection
(define-param R-pt (vector3 1 2 3)) ; point to measure reflection
3) For your run command, do something like:
(run-until
400
(at-beginning output-epsilon)
(in-volume (volume (center R-pt) (size 0 0 0))
(to-appended "R" (output-component R-comp))))
(print "dt = " (meep-fields-dt-get fields) "\n")
What this will do is to output an HDF5 file R.h5 containing a
1-dimensional data set given by the field R-comp at R-pt as a function of
time.
You should replace "400" with however long it takes the fields to decay.
In principle you could use (stop-when-decayed R-comp R-pt 1e-3), but I put
in a fixed time so that both your simulation and the normalization run
(below) take the same number of time steps.
The last line is to print out the time-step size, delta t, which will be
needed below.
4) Then, you will run the simulation *twice*, once with your structure,
outputting R.h5, and once with *no* scattering structure, outputting
R0.h5. The latter file is used for normalization.
5) Read both files into some program, e.g. Matlab or your favorite
spreadsheet (via h5totxt). e.g. in Matlab:
R = hdf5read("R.h5", "ez");
R0 = hdf5read("R0.h5", "ez");
Also, paste in the value of "dt" printed out above, into a variable dt.
6) Now, to compute the reflected phase as a function of frequency, you
just need to compute the Fourier transform of R - R0. In Matlab:
N = length(R);
Rangle = arg(fft(R - R0));
Rangle = Rangle(1:N/2); % frequencies above Nyquist are redundant
freq = [0:length(Rangle)-1] / (N * dt);
plot(freq, Rangle)
The final command should plot the phase angle (in radians) vs. frequency,
for the usual Meep frequency units.
I'm just typing the above commands off the top of my head, so I may have
made a few mistakes, but hopefully you get the idea.
Cordially,
Steven G. Johnson
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
