UserGuide
                        RectDisp  1.4
                       Vadim A. Markel
                      vmarkel@upenn.edu

Reference: V.A.Markel, M.Schöbinger and K.Hollaus, A fast method to
compute dispersion diagrams of three-dimensional photonic crystals
with rectangular geometry, Comp.Phys.Comm. 279, 108441 (2022).

https://doi.org/10.1016/j.cpc.2022.108441

CONTENTS
0. Introduction
1. Environment, Language and Compiling, Parallelism
2. Input and Output
3. Parameter File RectDisp.par
4. Defining Permittivities
5. Output File disp.dat
6. Units
7. Application Notes
8. Demo Directory
9. Limitations
0. Introduction
===============

0.1. Codes in this package are provided as is without any implicit or
explicit warranty. Execution of the codes depends on several
parameters. It is the responsibility of the user to select the correct
values for these parameters, verify convergence, etc.

0.2. Codes are free to use to any non-profit organization and
individual researcher working for a nonprofit organization. Commercial
use is prohibited.

0.3. If these codes are used in any scientific publication,
attribution and/or notification of the author would be appreciated.

0.4. This UserGuide is not exhaustive; in case of questions or if you
need help, please contact the author at the e-mail above.

0.5. Current limitations:

     (a) Only propagation in the crystallographic planes XZ or YZ or
     strictly along Z-axis are currently supported. This means that
     either of three cases are allowed (see Point 5.2 for definition
     of parameters fx, fy): 
       (i)  fx=fy=0
       (ii) fx=/=0 and fy=0
       (ii) fx=0   and fy=/=0
       
     (b) Several useful options for the root-searching algorithms are
     currently not adjustable through parameter selection.

     (c) Dispersion formulas for the inclusion and host permittivities
     are not adjustable through parameters. Instead, one can edit the
     files cfun_i.f and cfun_h.f directly to change the laws of
     dispersion. The form of these files is very
     self-explanatory. Considering that there are many ways in which
     dispersion can be defined, it is too tedious to control through
     input parameters. An alternative approach that allows full
     flexibility is to read the permittivities from file by using
     PERM_i='FILE' or PERM_h='FILE' (see Section 3, description of
     Parameters [2] and [3]).


1. Environment, Language and Compiling, Parallelism
===================================================

1.1. The codes have been developed and tested under Linux. Codes were
not tested under Windows or Mac OS, but are expected to work with
proper compilation.

1.2. Codes do not use any external libraries.

1.3. Codes are written in modern Fortran. Compilation was tested using
Intel's Fortran (ifort) and Gnu Fortran (gfortran). Intel's compiler
is preferred and yields better performance and easier/more efficient
parallelization.

1.4. To compile the codes using ifort, simply change to the root
directory ~/RectDisp/ and type "make all <ENTER>" or "make RectDisp
<ENTER>". To clean the directory and re-compile from scratch, type
"make clean <ENTER>" and then "make all <ENTER>".

1.5. To compile using gfortran, copy the file Makefile.gnu to Makefile
using the command "cp Makefile.gnu Makefile <ENTER>". Then type "make
all <ENTER>".

1.6. In case of successful compilation, the executable will be placed
in ~/RectDisp/bin/. From ~/RectDisp/work/, run the code by typing
"../bin/RectDisp.exe <ENTER>" 

1.7. If compiled with ifort, the executable will be parallelalized
automatically, and the resulting executable will be scalable. However,
scaling is much better for large values of Lx, Ly, Lz. For L=8, it is
not recommended using more than 2 or 3 threads. For L=64, 16 threads
can provide good performance.

1.8 Users can control the number of threads by setting the value of
the environmental variable OMP_NUM_THREADS with an appropriate shell
command. For example, under csh or tcsh, type something like "setenv
OMP_NUM_THREADS 16 <ENTER>" (16 here is an example). The computer
should have 16 physical (not logical) threads or else there will be no
expected scaling.

1.9. To parallelize under gfortran, uncomment the line "lopt =
-fopenmp". Then go to the subdirectory ~/RectDisp/sub/ and edit the
file mult_W.f. The user will need to insert appropriate OpenMP
directives to parallelize the loops. Scalability was not verified and
is not guaranteed with gfortran, and parallelization effectiveness
will depend on User's expertise.


2. Input and Output
===================

2.1. The following file is required in the working directory:
RectDisp.par. An example of this file is placed in
~/RectDisp/work/. This file contains several parameters to control the
problem set-up and execution. The parameters must be readable and in
some sense correct. All parameters are checked for sanity at
run-time. If some meaningless parameters are entered (example: a
negative number of frequency points), the code will report an error
and quit.

2.2. If PERM_i='FILE', then the input file eps_i.dat is required in
the working directory. If PERM_h='FILE' the input file eps_h.dat is
required in the working directory. The files, when required, must have
the correct format (see Section 4).

2.3. If PERM_i='DISP', then the file eps_i.dat will be written to the
working directory.  If PERM_h='DISP', then the file eps_h.dat will be
written to the working directory. These files can be used in the
future as input files.

2.4. The PC dispersion data are written to the file disp.dat. This is
the main result of the codes.

2.5. Convergence data are written to the file conv.dat. These data are
for development and debugging and users should not display or
interpret the numbers in conv.dat.


3. Parameter file RectDisp.par
==============================

3.1. IMPORTANT: when editing this file:
 
   a) Do not use tabs, use single spaces only

   b) The text after ! symbols are comments. Do not move ! symbols;
   all actual parameters should fit in the space before the ! symbols. 

3.2. The file RectDisp.par contains parameters that will control the
problem set-up and execution of the codes. The example file has
comments that should be self-explanatory. The following is a more
in-depth description. Numbers below correspond to the sequential
parameter numbers in ~/RectDisp/work/RectDisp.par.

[1] Method : character*8 : The method to solve the nonlinear equation.
            'rEPS' -- hybrid method to be used with purely real
                      permittivities 
            'fixP' -- fixed-point iteration
            'ratA' -- rational approximation

    These options cover all the possibilities implemented so
    far. However, some spelling variations are allowed; see file
    ~/RectDisp/inc/values.f for details.

[2] PERM_i : character*8 : The method for defining permittivity of inclusions.
            'DISP'  -- dispersion formula (see Section 4)
            'FIXED' -- fixed values read from RectDisp.par
            'FILE'  -- read values from file perm_i.dat
    
    These options cover all the possibilities implemented so
    far. However, some spelling variations are allowed; see file
    ~/RectDisp/inc/values.f for details.


[3] PERM_h : character*8 : The method for defining permittivity of host.
            'DISP'  -- dispersion formula (see Section 4)
            'FIXED' -- fixed values read from RectDisp.par
            'FILE'  -- read values from file perm_h.dat

    These options cover all the possibilities implemented so
    far. However, some spelling variations are allowed; see file
    ~/RectDisp/inc/values.f for details.

[4] ceps_i_fix : complex*16 : Fixed permittivity of inclusions to be
                        used if PERM_i='FIXED'

	Enter in the form (2.0, 0.20) where the first and second
        number are the real and imaginary parts of eps_i


[5] ceps_h_fix : complex*16 : Fixed permittivity of the host to be
                        used if PERM_h='FIXED'

	Enter in the form (2.0, 0.20) where the first and second
        number are the real and imaginary parts of eps_h

[6] Pol : character*8 : Polarization label. 

       Is used to select the polarization mode as described in the
       paper. Following values can be given but are not always
       consistent with other parameters (code will let you know if
       not):

       's', 'p', 'x', 'y', 'xz', 'yz'

       For example, if propagation is in XZ plane, the labels 'xz' and
       'p' are equivalent. However, if propagation is strictly along
       Z, then the labels 's' or 'p' are ambiguous; give 'x' or 'y',
       which can yield different results in anisotropic PCs.  For
       spelling variations, see ~/RectDisp/inv/values.f

[7] hx, hy, hz : real*8 : Dimensions of unit cell in arbitrary units.
                          See Section 6 for more information about units.

[8] Lx, Ly, Lz : integer*4 : Size of the reciprocal lattice grid. 

        This is the most important set of parameters controlling
        convergence and performance. The total number of plane waves
        used for the basis is (2*Lx+1)(2*Ly+1)(2*Lz+1). If Lx=Ly=Lz=L,
        increasing L by the factor of 2 will increase the computation
        time by the factor of 16. Memory requirements will increase by
        the factor of 8. However, the codes use very little memory and are
        more likely to become slow before users run out of memory.

[9] ax, ay, az : real*8 :  Size of inclusion in same units are
                           Parameters hx, hy, hz. 

        IMPORTANT: the inclusion size is (2*ax)*(2*ay)*(2*az). Therefore,
        it is required that 2*ax <= hx. If 2*ax=hx, the medium is
        continuous in the X-direction and it is sufficient to set
        Lx=0; same for other directions. See Point 7.4.

[10] rk_min, rk_max : real*8 : Minimum and maximum normalized frequencies
                               for a scan. 

        Here rk = k*hz/pi, where k=omega/c is the free-space wave
        number. Thus, these quantities are units-independent and
        dimensionless. Not used if either PERM_i='FILE' or
        PERM_h='FILE'

[11] fkx, fky : real*8 : Dimensionless parameters defining projection
                         of the Bloch wave number onto the XY
                         plane. 

       At any frequency used, we have kx = k*fkx, ky=k*fky where
       k=omega/c. The Bloch wave vector is of the form (kx, ky, qz).
       At every frequency considered, kx and ky are defined from the above
       condition and qz is computed.

[12] Np : integer*4 : Number of frequency points. 

       If PERM_i='FILE' or PERM_h='FILE', the actual number of
       frequency points can be smaller than Np and defined by the
       number(s) of lines in eps_i.dat and eps_h.dat. But the number
       of frequency points will never be larger than Np: if any of the
       above files has more lines, only the first Np lines will be
       read and used.

[13] Nord : integer*4 : Order of continued fraction. 

       Value of 16 is almost always sufficient. Can be set however to
       a generous value (say, 64) without significant loss of
       performance. The code will determine if the expansion converges
       at smaller Nord and will not perform calculations that are not
       needed.

[14] eps : real*8 : Precision with which the equation should be
                    satisfied, i.e., for solving F(qz)=0 the stop
                    condition is |F(qz)|<=eps. 

       Usually the value eps=0.001 or eps=0.0001 is sufficient for
       good precision. Making eps smaller might not improve the
       precision of obtained solutions due to other imprecision
       involved.


[15] itm_1, itm_2, itm_3 : integer*4 : Maximum numbers of iterations used 
       in iterative solution of the nonlinear equation. 

       For Method='fixP' or Method='ratA', only itm_1 is used. For
       Method='rEPS', all three parameters are used. Recommended
       values are 50,25,25, but, in most cases, this is excessive. The
       parameters are supplied so that the code does not loop
       indefinitely in the case when the solution can not be found.

[16] Nqth : integer*4 : Minimum distance from the edge of first
                        Brillouin zone (FBZ) in units of 2*pi/hz for which
                        to re-check precision and run additional
                        iterations of solver as needed.

       If a solution qz is found such that |Re(q)| > pi/hz, the code will
       translate it to FBZ as 
          qz <-- qz - nq*2*pi/hz 
       where 
          nq = floor[ (Re(qz) + pi/hz) / (2*pi/hz) ]
       Theoretically, if there is a root at qz, there is also a root
       at any qz + n*2*pi/hz, where n is an integer. But due to the
       truncation used in numerical algebra, this property does not
       hold precisely. Therefore, if |nq| > Nqth, the code will
       recheck precision and, if necessary, run additional iterations
       to make sure that the precision goal is met.

          Nqth=0 --- Precision goal will be enforced for all points in
                     disp.dat. 
          Nqth=1 --- Some points in disp.dat may have not met precision
                     goal. However, the difference with Nqth=0 is most
                     often small and not visually discernible. Faster
                     execution compared to Nqth=0
          Nqth=2 --- Even more potential to lose precision but can be
                     not different from Nqth=1
          Nqth>2 --- A warning will be issued about possible loss of
                     precision

       Under many circumstances, Nqth=1 provides the best balance
       between performance and precision. If some points look off,
       change to Nqth=0 (but a fix is not guaranteed). Selecting
       Nqth>2 is not recommended but this option is available.

       By default, Nqth=1 is set in all parameter files supplied with
       this package; it is recommended to change this setting only with
       full understanding of what this entails.


4. Defining Permittivities
==========================

4.1. There are three methods to define the permittivities:

     a) Fixed values from RectDisp.par

     b) Dispersion formula coded in subroutines cfun_i.f and cfun_h.f

     c) Reading values from input file(s)

4.2. The functions coded in cfun_i.f and cfun_h.f can not be tuned by
changing parameters. Edit these files directly if you wish to change
the dispersion formula. This should be very easy and self-explanatory.

4.3. If either fixed permittivities or dispersion formulas are used
for both host and inclusions, then the frequencies for the scan are
determined using the parameters rk_min, rk_max and Np from the
parameter file RectDisp.par, according to:

   drk = (rk_max-rk_min)/(Np-1)
   for i=1,Np : rk = rk_min + drk*(i-1)

Here rk = k*hz/pi, where k=omega/c. Thus rk is by definition
dimensionless.

4.4. Dispersion formulas in cfun_i.f and cfun_h.f compute complex
permittivities eps_i and eps_h as functions of the dimensionless
frequency rk. The formula is

eps(rk) =  1 + w0*w0*(eps(0)-1) / (w0*w0 - rk*rk - i*g*w0*rk)

where w0 is the dimensionless resonance frequency, w0=omega_0*hz/pi, i
is the complex unit, eps(0) is the static permittivity (the limit when
rk-->0) and g=gamma/omega_0 is the dimensionless parameter specifying
the width of the resonance. Values of w0, eps(0) and g are specified
in the files cfun_i.f and cfun_h.f for the inclusions and the host
separately. All these constants are dimensionless.

4.5. If at least one of the permittivities is read from file, the
parameters rk_min and rk_max are not used to determine the frequency
samples. Specifically, the following approach is used:

a) If PERM_i='FILE' and PERM_h='FILE' (both permittivities are read
from files):

   i)   The code will attempt to read Np first lines from both files. If
   one of the files has Mp<Np lines, only Mp lines will be read from
   both files.

   ii)  Each line read consists of three numbers: rk, eps_re, eps_im

   iii) The values of rk read from similarly-numbered lines in
   eps_i.dat and eps_h.dat must match exactly; otherwise the code will
   quit. For example, assume that we read:
    
   From k-th line of eps_i.dat: x1, y1, z1
   From k-th line of eps_i.dat: x2, y2, z2

   The condition x1==x2 is checked; if it does not hold, the code
   quits. If the condition holds, then:
      rk    <-- x1
      eps_i <-- y1 + i*z1
      eps_h <-- y2 + i*z2

b) If PERM_i='FILE' xor PERM_h='FILE' (exclusively, only one of the
permittivities is read from files):

   i)   The normalized frequency is read from the first column of this
   file
 
   ii)  First Np lines from the file are read; or Mp lines if Mp<Np

   iii) Permittivity that is read from file is defined as above; the
   other permittivity is computed either by dispersion formula (using rk read
   in Point (i)) or is constant depending on whether it is defined as
   'DISP' or 'FIXED'.

4.6. Format of files eps_i.dat, eps_h.dat is the following. Each file
consists of three columns of numbers either comma-separated or not.

   a) The first column is the value of normalized frequency
   rk=k*hz/pi, k=omega/c

   b) The second column is the real value of the permittivity (at that
   frequency)

   c) The third column is the imaginary value of the permittivity (at
   that frequency)


5. Output File disp.dat
=======================

5.1. The main result of the computations -- the dispersion points qz
as functions of the normalized frequency rk=k*hz/pi -- are written to
the file disp.dat.

5.2. The particular restriction of the Bloch variety considered by the
codes is the following: For each normalized frequency
rk=k*hz/pi=omega*hz/pi*c, we seek the Bloch wave vector in the form
q=(kx, ky, qz) where kx=fx*k, ky=fy*k and qz is computed by the
codes. Note that k=omega/c, and the normalized projections are
rkx=fx*rk, rky=fy*rk. Dimensionless parameters fx, fy are read from
RectDisp.par. Thus, for each frequency, the components kx and ky of
the Bloch wave vector are pre-determined and not computed. This is
usually explained by the boundary conditions at the z=0 interface,
which must be considered when a plane wave enters a PC slab through
that interface. See the paper for more details of the problem
formulation.

5.3. The file disp.dat contains three columns:

   a) First column is normalized frequency, rk

   b) Second column is Re(qz)*hz/pi -- the normalized real part of the
   z-component of Bloch wave vector

   c) Third column is Im(qz)*hz -- normalized imaginary part of the
   z-component of Bloch wave vector, but without the 1/pi factor

5.4. All numbers in disp.dat are by definition dimensionless. However,
the real and imaginary parts of qz are differently normalized. Re(qz)
is divided by pi/hz -- the right bound of the first Brillouin
zone. The imaginary part is divided by 1/hz. Thus, the quantity
Im(qz)*hz in the third column of disp.dat quantifies the decay of the
Bloch wave upon propagation by one unit cell in the z-direction:
|E(z+hz)| = exp(-Im(qz)*hz)*|E(z)|.

5.5. Data points are written into the disp.dat file only if the
solution at a normalized frequency rk was found. If not, no data point
will be written for this rk. Therefore, the file disp.dat may have
gaps, that is, missing frequency points for some values of rk that one
might otherwise expect, i.e., from the frequency sampling according to
Point 4.3.

5.5. If a solution is found for some rk, then two data points are
always written to disp.dat representing qz and -qz (both are
solutions). For this reason, do not try to connect the data points by
lines in plots. If you need to connect the points, remove first the
unwanted lines (say, with negative real parts Re(qz)) from
disp.dat.


6. Units
========

6.1. The only dimensional parameters that must be specified are hx,
hy, hz and ax, ay, az. These parameters must be given in the same
units. However, these units are arbitrary.

6.2. The parameter rk defined either through rk_min and rk_max or from
input files is dimensionless: rk=k*hz/pi = omega*hz/pi*c =
2*hz/lambda. Therefore, rk depends not only on omega (lambda) but also
on hz. The only place where this must be taken into account is in the
dispersion formulas or, more generally, when establishing
correspondence between rk and eps_i(rk), eps_h(rk).

6.3. For example, the resonance frequency in cfun_i.f is set as
w0=2.0. This means that, for the resonance frequency omega_0 or the
resonance wavelength lambda_0, we have:

             omega_0*hz / pi*c = 2*hz/lambda_0 = w0 = 2.0


7. Application Notes
====================

7.1. The most important parameters that control the performance are
Lx, Ly and Lz. If Lx=Ly=Lz=L, and we increase L by the factor of 2,
the computation time will increase by the factor of 16.

7.2. However, to verify convergence, it makes sense to increase L by
at least the factor of 2. Comparing L=16 and L=17 is not meaningful.

7.3. The parameter Nord has some influence on computation time but is
not very important. You can use a generous value, say Nord=64; the
codes will determine the depth of continued fraction that is
sufficient for convergence and will not perform unnecessary
computations.

7.4. Two-dimensional and one-dimensional problems can be solved by
setting ax=0.5*hx or ay=0.5*hy or both. If ax=0.5*hx, the medium is
continuous in the x-direction and, therefore, two-dimensional. In this
case, it is sufficient to set Lx=0. Codes will run much faster. Same
holds in the case ay=0.5*hy. If both conditions hold, then the PC is a
one-dimensional layered medium. One can set Lx=Ly=0 in this case. For
one-dimensional media, analytical solution is available, which can be
used as a check.

7.5. There are three solution Methods that can be set in the parameter
file. Selecting the method depends on the application. Most difficult
case is when both permittivities are real. However, computation of the
pass bands still will work fine (and fast) even in this case with
Method='fixP'. The difficulties arise in the band gaps and evanescent
mode bands. In those spectral regions one can use Method='rEPS'. If
permittivities are complex, using Method='ratA' is a good first
choice.

7.6. Still, some solutions might not be found by the methods so far
implemented. This becomes more of a problem at higher frequencies. The
first two pass bands are usually computed stably, as well as the band
gaps around them.

7.7. The codes have some additional adjustments that can be used to
improve root-finding, but currently we do not provide tuning
capabilities through the parameter file. This capabilities are used
for development and may be added in the future.

7.8. Sometimes, a computed data point looks off, which usually means
that it falls out of an otherwise smooth line. There can be several
reasons for this such as:

   a) Shallow minimum of |F(qz)| (a false solution)

   b) An instability in the continued fraction expansion. Sometimes,
   the continued fraction tails contain a 0/0 uncertainty, and we
   currently flush it to zero. Opening the uncertainty is possible
   analytically, but considering that such occurrences are rare, and
   required programming is somewhat involved, we did not yet implement
   this capability.

   c) A sharp resonance and/or insufficient value of L. This can be
   checked by changing L. Not a frequent occurrence and more typical
   at higher frequencies.

Debugging these instabilities may be difficult. However, when a point
is obviously off, some additional checks can be performed by zooming
in at a narrow frequency interval and running additional simulations
with different parameters, i.e., with larger L or a different Method.

7.9. Codes were preliminarily tested for Drudean metals and appeared
to work, but checks made so far are insufficient. Author makes no
claims regarding applicability to metal inclusions.

7.10. As any numerical method, the method implemented in this package
can be broken. The easiest way to break it is to request computations
at very high frequencies or unrealistically large L. Start from
small-to-moderate frequencies and L and then expand the rand of these
parameters gradually.


8. Demo Directories
===================

8.1. The subdirectory ~/rectDisp/demo/ has several 2D and 3D
examples. These examples were selected so that converging detailed
results can be obtain on a typical laptop in 2 to 3 minutes. 

8.2. To run the examples, navigate the directory tree, which is
structured as ~/RectDisp/demo/nD/ExampleName/, where n=2,3 refers to
dimensionality and ExampleName is something like SmallCube. Read the
README file in each subdirectory. Run the example by typing
"./RectDisp.exe <ENTER>".

8.3. Note that functionality of the demo examples relies on several
symbolic links. Files to which these links point should exist. All
such files are within the directory tree with the root ~/rectDisp/. In
particular, the executable RectDisp.exe should exist in
~/RectDisp/bin/ subdirectory.

8.4. Each subdirectory ~/ExampleName/ has a further subdirectory
~/Data/, where pre-computed data and graphics produced from these data
are placed. Users can reproduce the data by running the example as
explained in Point 8.2. To re-create the graphics from data, type
"./cmd <ENTER>" from the ~/Data/ subdirectory. This assumes that the
original data files are in place, or user-created data files with the
same names are in place.


9. Limitations and Future Development
=====================================

9.1. At higher frequencies, there might exist multiple solutions for
qz at a given normalized frequency rk. This usually happens above the
second band gap. The code currently finds only one solution for a
given frequency and it can therefore "jump" from one branch to
another. The resultant dispersion can look noisy or irregular. In
future development, this issue will be addressed by applying
polynomial approximation to the root problem and finding all solutions
qz in the FBZ.

9.2. More general propagation directions or restrictions of the Bloch
variety will be implemented in the future.