#ifndef GENGWED_H
#define GENGWED_H

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

/* 
 * gengwed.h
 * Electrodynamic constants, basic relations, general definitions 
 * concerning 2D optical wave propagation and eigenmode analysis 
 */

#include "matrix.h"
#include "structure.h"

/* compiler xlC on IBM SP needs this statement ?!? */
#undef HZ

#define PI	3.14159265358979323846
#define PI_2	1.57079632679489661923
#define PI_4	0.78539816339744830962
#define PIM2    6.28318530717958647693

/* vacuum speed of light [um/fs] */
#define	SOLC0       2.99792458E-1
/* 1/sqrt(mu0) [sqrt(A*um/(V*fs))] */
#define INVSQRTMU0  2.82094792E-2 
/* 1/sqrt(ep0) [sqrt(V*um/(A*fs))] */
#define	INVSQRTEP0 10.62736593 

/* vacuum permittivity [10^(-15) A s/(V um) = A fs/(V um) ] */
#define EP0 8.854187817E-3
/* vacuum permeability [10^(-15) V s/(A um) = V fs/(A um) ] */
#define MU0 1.2566370614E+3 

/* vacuum wavenumber [1/um], function of the vacuum wavelength [um] */
inline double val_k0(double lambda) { return 6.283185307179586/lambda; }; 
/* 1/omega/epsilon_0 [V um / A], function of the vacuum wavelength [um] */
inline double val_invomep0(double lambda) { return 59.95849160*lambda; };
/* 1/omega/mu_0 [A um / V], function of the vacuum wavelength [um] */
inline double val_invommu0(double lambda) { return (4.224638618e-4)*lambda; }; 
/* angular frequency [radiants/fs], function of the vacuum wavelength [um] */
inline double val_omega(double lambda) { return 1.883651567/lambda; };
/* corresponding period in time [fs], function of the vacuum wavelength [um] */
inline double val_period_T(double lambda) { return lambda*3.335640952; }; 
/* vacuum wavelength [um], function of the angular frequency [radiants/fs] */
inline double val_lambda(double omega) { return 1.883651567/omega; };

/* ----------------------------------------------------------------------- */

/* passive cavities: translation between characterizing quantities;
   given: complex eigenfrequency omega [rad/fs], or its imaginary part alpha */

/* Q factor */
double theQ(Complex omega); 

/* cavity lifetime [fs] */
double tauC(double alpha);

/* outgoing radiation, power vs. frequency, FWHM [rad/fs] */
double FWHM_omega(double alpha);

/* outgoing radiation, power vs. wavelength, FWHM [mum] */
double FWHM_lambda(Complex omega);

/* ----------------------------------------------------------------------- */

/* Cartesian and cylindrical field components */ 
enum Fcomp {EX, EY, EZ, HX, HY, HZ, SX, SY, SZ, ER, EP, HR, HP, SR, SP, ET, HT, ST};

/* Fcomp -> field character 'E', 'H', 'S' */
char fldchr(Fcomp cp);

/* Fcomp -> comp character 'x', 'y', 'z', 'r', 'p', 't' */
char cpchr(Fcomp cp);

/* 'EHS''xyrzpt' -> Fcomp */
Fcomp fcomp(char eh, char cp);

/* field "strength" identifiers, for plotting purposes  
   mE: |(Ex, Ey, Ez)|, vectorial electric field, absolute value
   mH: |(Hx, Hy, Hz)|, vectorial magnetic field, absolute value
   qE: |(Ex, Ey, Ez)|^2, vectorial electric field, absolute square
   qH: |(Hx, Hy, Hz)|^2, vectorial magnetic field, absolute square
   mW: electromagnetic energy density */
enum FSid {mE, mH, qE, qH, mW};

/* FSid -> character 'E', 'H', 'W' */
char idchr(FSid id);

/* FSid -> character 'm', 'q' */
char mqchr(FSid id);

/* polarization parameter */ 
enum Polarization {TE, TM};

/* Polarization -> 'e', 'm' */
char polchr(Polarization p);

/* Polarization -> 'E', 'M' */
char polCHR(Polarization p);

/* Polarization -> principal field component */
Fcomp principalcomp(Polarization p);

/* Directions of propagation */
enum Propdir {FORW, BACK};

/* Boundary types for planar mode analysis */
enum Boundary_type {OPEN, LIMD, LIMN};

/* Types of planar modes */
enum Propagation_type {PROPAG, EVANESC};

/* evaluation of integrals: choice of principal 
   field or its derivative */
enum FldorDer {FLD, DER};

/* ----------------------------------------------------------------------- */

/* Attributes for field output  
   MOD: absolute value 
   ORG: the original field 
   SQR: the squared field
   REP: real part 
   IMP: imaginary part */
enum Afo {MOD, ORG, SQR, REP, IMP};

/* convert complex field values f to real numbers, 
   depending on the attribute foa */
double  realize(Complex f, Afo foa);
Dvector realize(Cvector f, Afo foa);
Dmatrix realize(Cmatrix f, Afo foa);

/* Afo -> 'm', 'o', 's', 'r', 'i' */
char afochr(Afo foa);

/* Afo + Fcomp -> '|Hx|', 'Re Ey', 'Im Hz'; filled into desc */
int afocpdesc(Afo foa, Fcomp cp, char *desc);
/* FSid -> '|E|', '|E|^2', 'W', etc.; filled into desc */
int iddesc(FSid id, char *desc);

/* ----------------------------------------------------------------------- */

/* permittivity perturbations, planar waveguides, 1D cross sections */
class Perturbation
{
public:
	// part of cross section with nonzero perturbation, perturbed region
	Interval pr;          

	// entries of the permittivity perturbation tensor, complex 3x3 matrix 
	Cmatrix e;

	// initialize 
	Perturbation();
	Perturbation(Interval iv);
	Perturbation(Interval iv, 
	             Complex p_xx, Complex p_xy, Complex p_xz,	
	             Complex p_yx, Complex p_yy, Complex p_yz,	
	             Complex p_zx, Complex p_zy, Complex p_zz);	
	// relevant quantities for isotropic absorbing part 
        // and Hermitian anisotropic part  
	Perturbation(Interval iv, 
	             double rxx, double ryy, double rzz, 
                     double rxy, double ixz, double iyz, double iabs);

	// input
	void read(FILE *dat);
	//output
	void write(FILE *dat);
};

/* ----------------------------------------------------------------------- */

/* special perturbations */

/* a change of the local refractive index:  
   r: the perturbed region on the waveguide cross section 
   n: uniform refractive index of the material in r
   deltan: the shift of the refractive index 
   The perturbation leads to an isotropic real first order permittivity shift 
   of 2*n*deltan */
Perturbation refindshift(double n, double deltan, Interval r);

/* isotropic absorbtion: 
   r: the lossy region on the waveguide cross section 
   alpha: intensity attenuation [1/micrometer] of the material in r,
          the bulk material damps the intensity of a plane wave according to
          exp(-alpha z), where z is the propagation distance
   n: refractive index of the material in r
   lambda: vacuum wavelength [micrometer]
   The attenuation leads to a complex refractive index n - i alpha lambda/4/pi,
   or a permittivity perturbation of - i n alpha lambda /(2 pi) */
Perturbation attenuation(double alpha, double n, double lambda, Interval r);

/* magnetooptic permittivity contribution:
   r: the magnetooptic region on the waveguide cross section 
   sFrot: specific Faraday rotation of the material [degrees/cm] (!)
   theta, phi: orientation of the magnetization in spherical coordinates,
               vec(M) = M(sin(theta)cos(phi), sin(theta)sin(phi), cos(theta)) 
   n: refractive index of the material in r
   lambda: vacuum wavelength [micrometer] */
Perturbation magopt(double sFrot, double theta, double phi, 
                    double n, double lambda, Interval r);

/* as above, for polar orientation of the magnetization, vec(M) || x-axis, */
Perturbation magopt_polar(double sFrot, double n, double lambda, Interval r);

/* for equatorial orientation of the magnetization, vec(M) || y-axis, */
Perturbation magopt_equat(double sFrot, double n, double lambda, Interval r);

/* and for longitudinal orientation of the magnetization, vec(M) || z-axis */
Perturbation magopt_long(double sFrot, double n, double lambda, Interval r);

/* ---------------------------------------------------------------------- */

/* permittivity perturbation, channel waveguides, 2D cross section */
class ChWgPert
{
public:
	// part of x-y-plane with nonzero perturbation 	
	Rect r;          

	// entries of the permittivity perturbation tensor, complex 3x3 matrix 
	Cmatrix e;

	// initialize 
	ChWgPert();
	ChWgPert(Rect r);
	ChWgPert(Rect r, 
	       Complex p_xx, Complex p_xy, Complex p_xz,	
	       Complex p_yx, Complex p_yy, Complex p_yz,	
	       Complex p_zx, Complex p_zy, Complex p_zz);	
	ChWgPert(Cmatrix p);
	// relevant quantities for isotropic absorbing part 
        // and Hermitian anisotropic part  
	ChWgPert(Rect r, double rxx, double ryy, double rzz, 
			     double rxy, double ixz, double iyz, double iabs);

	// input
	void read(FILE *dat);
	//output
	void write(FILE *dat);
};

/* ##################################################################### */

/* special perturbations */

/* isotropic absorbtion: 
   rl: the lossy area on the waveguide cross section 
   alpha: intensity attenuation of the material in rl, [1/micrometer]
          the material damps the intensity of a plane wave according to
          exp(-alpha z)
   n: refractive index of the material in rl
   lambda: vacuum wavelength [micrometer]
   The attenuation leads to a complex refractive index n - i alpha lambda/4/pi,
   or a permittivity perturbation of - i n alpha lambda /(2 pi) */
ChWgPert attenuation(double alpha, double n, double lambda, Rect rl);

/* magnetooptic permittivity contribution:
   rl: the magnetooptic area on the waveguide cross section 
   sFrot: specific Faraday rotation of the material [degrees/cm] (!)
   theta, phi: orientation of the magnetization in spherical coordinates,
               vec(M) = M(sin(theta)cos(phi), sin(theta)sin(phi), cos(theta)) 
   n: refractive index of the material in rl
   lambda: vacuum wavelength [micrometer] */
ChWgPert magopt(double sFrot, double theta, double phi, 
                    double n, double lambda, Rect rl);

/* as above, for polar orientation of the magnetization, vec(M) || x-axis, */
ChWgPert magopt_polar(double sFrot, double n, double lambda, Rect rl);

/* for equatorial orientation of the magnetization, vec(M) || y-axis, */
ChWgPert magopt_equat(double sFrot, double n, double lambda, Rect rl);

/* and for longitudinal orientation of the magnetization, vec(M) || z-axis */
ChWgPert magopt_long(double sFrot, double n, double lambda, Rect rl);

/* ----------------------------------------------------------------------- */

/* convert a complex permittivity eps into a complex refractive index n
   such that eps = n*n, eps != 0 required; 
   eps = eps'-i eps'', n = n'-i n'', sign convention: n' > 0 */
Complex cepston(Complex e);

/* ----------------------------------------------------------------------- */

/* Gaussian beams, input specification for the Helmholtz solvers,
   paraxial representation of moderately tilted beams, Gaussian profile 
   prescribed on the x-axis, propagation in the positive z-direction 
   with specific wavelength in a homogeneous dielectric medium */
class Gaussianbeam
{
public:
        Polarization pol; // TE, TM
        double w;         // 1/e-half-width (field) of the beam, and  
        double x0;        // beam center at the initial z-position
        double theta;     // tilt angle (degrees!) with respect to the z-axis 
        double lambda;    // vacuum wavelength,
        double ng;        // background refractive index. 

        // initialize 
        Gaussianbeam(Polarization tp, double tw, double tx0, double ttheta, 
                     double tlambda, double tng);

        // field profile at the initial z-position 
        Complex field(Fcomp cp, double x) const;

private: 
	// principal field component
	Complex phi(double x) const;
	// ... its z-derivative, according to the paraxial approximation
	Complex dzphi(double x) const;
	// ... its x-derivative, according to the paraxial approximation
	Complex dxphi(double x) const;
	// internal values
	double amp;   // amplitude, normalization
	double kncth; // k ng cos theta 
	double knsth; // k ng sin theta 
	Complex tef;  // i/omega/mu0
	Complex tmf;  // i/omega/epsilon0/ng/ng
};

/*
 * helper functions describing a scalar Gaussian beam propagating along 
 * the z-axis, half-width w
 */

Complex gbm_phi(double kn, double w, double x, double z);
Complex gbm_dxphi(double kn, double w, double x, double z);
Complex gbm_dzphi(double kn, double w, double x, double z);
Complex gbm_dxxphi(double kn, double w, double x, double z);
Complex gbm_dzzphi(double kn, double w, double x, double z);

/* ----------------------------------------------------------------------- */

/* a signed square root */
double sgn_sqrt(double x);

/* ----------------------------------------------------------------------- */

/* 
 * for computing parameter scans of resonant structures, narrow peaks in an
 * otherwise regular background:
 * given a range r of arguments and an initial stepsize h0, determine the 
 * next value of a series or arguments x (not ordered!, x is being extended 
 * by one value), such that initially appearing minima and maxima in the curve
 * x, y are resolved with a minimum stepsize of hm upon completion 
 * return value != 0: finished; 
 * ext_type: -1: resolve minima in y(x),
 * ext_type:  1: resolve maxima in y(x), default: resolve both;
 * minydev: minimum absolute deviation in y for an extremum to be identified
 * rough, inefficient, ...
 */
int next(Interval r, double h0, double hm, Dvector &x, Dvector y, int ext_type,
         double minydev);
int next(Interval r, double h0, double hm, Dvector &x, Dvector y, int ext_type);
int next(Interval r, double h0, double hm, Dvector &x, Dvector y);

#endif // GENGWED_H