#ifndef SLAMS_H
#define SLAMS_H

/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/*
 * slams.h
 * Mode solver routines, 
 * guided mode analysis & rigorous spectral dicretization
 */

/* mode solver options, default values, alter global setting, if necessary */

/* tolerance for fixing the propagation constants */
extern double SLAMS_BQTOL;
double bqTol(double bq);

/* minimum distance at which propagation constants are considered different */
extern double SLAMS_BQSEP;
double bqSep(double bq);

/* tolerance for the maximum relative error in the continuity of the 
   basic field component for accepting a field as a valid mode */
extern double SLAMS_DISCERRTOL;

/* mode spectrum computation, BEP/QUEP field discretization: 
   minimum difference between a (vaguely defined) "mode propagation angle" 
   and the transverse axis. If modes close to the transition from propagating 
   to evanescent type with almost vanishing propagation constant occur, 
   the spectral discretization procedure issues a warning message 
   ((!) in quiet mode), since this can lead to highly erroneous results. 
   Workaround: a slight change of the computational window width */
extern double SLAMS_MODEANGLELIM; 

/* --- mode analysis, helper functions --------------------------------- */
int trap(ModeArray& ma, Mode& m,
         double bqmin, double bqmax, double sf, int sm,
         Dvector& bqf, Ivector& found, Dmatrix& susp, int& nsc, int lfsc,
         int verb);
int trap(ModeArray& ma, Mode& m,
         double bqmin, double bqmax, double sf, int sm, int lfsc,
         Dvector& bqf, Ivector& found, int verb);
int recognizemode(Mode& m, int lop, Ivector& found, Dvector& bqf,
                  ModeArray& ma, int verb);

/* --- mode analysis --------------------------------------------------- */

/* determine all modes with squared propagation constants between 
   bqmin and bqmax, for the problem specified by wg, pol, 
   boundaries of type bdtb at bpb, of type bdtt at bpt; 
   verb: 0: no control output, 1: progress, 2: all */
void findmodes(Waveguide wg, Polarization pol, 
               Boundary_type bdtb, double bpb, Boundary_type bdtt, double bpt,
               double bqmin, double bqmax, ModeArray& ma, int verb);
/* nstd: 1: solve nonstandard  effective VEIMS problem,
   ec: VEIMS effective permittivity */
void findmodes(Waveguide wg, Polarization pol,
               Boundary_type bdtb, double bpb,
               Boundary_type bdtt, double bpt,
               double bqmin, double bqmax, ModeArray& ma, int verb,
               int nstd, Dvector ec);

/* homogeneous medium, constant refractive index:
   determine mnm modes up to order mnm-1, for the problem specified by
   wg, pol, with boundaries of type bdt at bpb and bpt;
   verb: 0: no control output, 1: progress, 2: all */
void freespacemodes(Waveguide wg, Polarization pol,
                    Boundary_type bdt, double bpb, double bpt,
                    int mnm, ModeArray& ma, int verb);

/* --- mode analysis, simpler driver routines -------------------------- */

/* complete guided mode analysis, returns the number of found modes 
   wg: the structure under consideration
   pol: polarization type, TE or TM
   ma: the modes found (output)
   quiet == 1: suppress log output */
int modeanalysis(Waveguide wg, Polarization pol, ModeArray& ma, int quiet); 
int modeanalysis(Waveguide wg, Polarization pol, ModeArray& ma); 

/* compute the discretized mode spectrum on an interval with artificial 
   boundaries, returns the number of found modes 
   wg: the structure under consideration
   pol: polarization type, TE or TM
   cw: the computational window
   bt: the boundary type, one of LIMD or LIMN
   mnm: highest order mode (+1) to be included, determine up to mnm eigenmodes 
   ma: the modes found (output)
   quiet: == 1: suppress log output */
int modespectrum(Waveguide wg, Polarization pol, Interval cw, Boundary_type bt,
                 int mnm, ModeArray& ma, int quiet);
int modespectrum(Waveguide wg, Polarization pol, Interval cw, Boundary_type bt,
                 int mnm, ModeArray& ma);

/* compute the discretized mode spectrum on an interval with Dirichlet  
   boundaries, returns the number of found modes 
   wg: the structure under consideration
   pol: polarization type, TE or TM
   cw: the computational window
   mnm: highest order mode (+1) to be included, determine up to mnm eigenmodes 
   ma: the modes found (output)
   quiet: == 1: suppress log output */
int modespectrum(Waveguide wg, Polarization pol, 
                 Interval cw, int mnm, ModeArray& ma, int quiet);
int modespectrum(Waveguide wg, Polarization pol, 
                 Interval cw, int mnm, ModeArray& ma);

#endif // SLAMS_H