#ifndef STRUCTURE_H
#define STRUCTURE_H

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

/* 
 * structure.h
 * Multilayer slab waveguides, waveguide sequences,
 * rectangular optical circuits;
 * cross section (2D) of 3D rectangular waveguides (VEIM solver)
 */

#include <stdio.h>
#include "matrix.h"

#define AWAY           1.0e+20
#define HDIST          1.0e-8
#define LIMBDSEP       0.001
#define COMPTOL_IV     1.0e-7
#define COMPTOL_N      1.0e-7
#define COMPTOL_HX     1.0e-7
#define COMPTOL_LAMBDA 1.0e-7
#define MINREFIND      1.0
#define MAXREFIND      10.0

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

/* a normalized parameter that specifies whether localized fields related 
   to pieces of a waveguide are used to establish approximate modal 
   basis sets */
extern double DECOUPLED;

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

/* interval from the waveguide cross section */
class Interval 
{
public:
	// lower boundary
	double x0;
	// upper boundary 
	double x1;
	// initialize
	Interval();
	Interval(double xa, double xb);
	// output
	void write(FILE *dat) const;
	// input
	void read(FILE *dat);
	// compare the present object to i, return 1 if equal, 0 otherwise
	int equal(const Interval& i) const;
	// length
	inline double len() { return x1-x0; }
	// restrict iv to the part inside the present object, 
	// return value: 0: no overlap, 1: something from iv remains
	int clip(Interval& iv) const;
	// test whether x is between x0 and x1
	inline int in(double x) { return (x0 <= x && x <= x1) ? 1 : 0; }
};

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

/* waveguide geometry + permittivity profile; vacuum wavelength ... */
class Waveguide 
{
public:
	// number of inner layers
	int nx;
	// positions of the dielectric interfaces
	Dvector hx;

	// get interval index corresponding to position x  
	int layeridx(double x) const;
	// get layer boundaries corresponding to index l  
	Interval layer(int l) const;
	// get layer boundaries corresponding to position x  
	Interval layer(double x) const;

	// vacuum wavelength 
	double lambda;
	// refractive index profile
	Dvector n;
	// permittivity on layer l
	double eps(int l) const;
	// permittivity at position x
	double eps(double x) const;

	// lower bound for effective mode permittivity, default value
	double defaultepseffmin() const;
	// upper bound for effective mode permittivity, default value
	double defaultepseffmax() const;

	// translate: hx -> hx+dx
	void translate(double dx);

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

	// initialize
	Waveguide();
	Waveguide(int vnx); 
	Waveguide(int vnx, Dvector vhx, double l, Dvector n);

	// free allocated memory
	void strip();

	// compare the present object to wg, return 1 if equal, 0 otherwise
	int equal(const Waveguide& wg) const;

	// test for constant refractive index, return 1, if const., 0 otherwise
	int constmed() const;

	// test matching of interface positions with another waveguide,
	// all boundaries of the current object should appear in w
	int bdmatch(Waveguide w) const;

	// remove unnecessary boundaries from waveguide,
	// down to a minimum of minnx inner layers
	void optimize(int minnx);
	// ... down to a minimum of one inner layer
	void optimize();

	// check consistency of the waveguide object: error&exit, if not ok
	void consistency() const;

	// check the symmetry of the waveguide, return value l 
	// l == 0: no symmetry 
	// l >= 1: symmetric structure with respect to layer l
	int checksymmetry() const;

	// split the present waveguide into two, at layer sl 
	// the pieces are stored in wg1 and wg2
	void split(int sl, Waveguide& wg1, Waveguide& wg2) const;
	
	// expand the waveguide towards a structure 
	// that is symmetric with respect to xs; xs > hx(nx) 
	Waveguide expand(double xs);

	// mirror the waveguide at xs; xs > hx(nx) 
	Waveguide mirror(double xs);

	// determine a limit for the effective permittivity, beyond which 
	// the inner layer l optically decouples the neighbouring structures
	double decoupepseff(int l);

	// determine a limit for the effective permittivity, beyond which 
	// the lower boundary at position p becomes irrelevant 
	double lbdecepseff(double p);

	// determine a limit for the effective permittivity, beyond which 
	// the upper boundary at position p becomes irrelevant 
	double ubdecepseff(double p);
	
	// refractive index profile -> MATLAB .m file  
	// disp: output region on the x axis; ext0, ext1: filename id chars 
	void plot(Interval disp, char ext0, char ext1) const;
	void plot(char ext0, char ext1) const;

	// special != 0: n stores permittivity values in place of refractive
	// indices, consistency checks for ref.ind. and perm. are disabled;
	// a misuse of the Waveguide container for problems with negative 
	// permittivity
	int special;
};

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

/* rectangular 2D optical circuit, 
   geometry & permittivity profile, vacuum wavelength;
   (so far) meant as a help to define a SegWgStruct */

class Circuit 
{
public:
	// number of inner horizontal layers
	int nx;
	// positions (x) of layer interfaces 
	Dvector hx;
	// number of inner vertical slices
	int nz;
	// positions (z) of slice interfaces 
	Dvector hz;
	
	// vacuum wavelength
	double lambda;
	// refractive indices 
	Dmatrix n;

	// initialize
	Circuit(); 
	Circuit(int vnx, int vnz); 

	// get rectangle index (l, m) corresponding to position (x,z)  
	void rectidx(double x, double z, int& l, int& m) const;

	// refractive index profile -> MATLAB .m file  
	// ext0, ext1: filename id chars 
	// xbeg, xend, zbeg, zend: display window in the x-z-plane
	void plot(double xbeg, double xend, double zbeg, double zend,
	          char ext0, char ext1) const;
	void plot(char ext0, char ext1) const;

	// special != 0: n stores permittivity values in place of refractive
	// indices, consistency checks for ref.ind. and perm. are disabled;
	// a misuse of the Circuit container for problems with negative 
	// permittivity
	int special;
};

/* port interfaces of a Circuit */
enum CBorder {RIGHT, TOP, LEFT, BOTTOM};

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

/* segmented planar waveguide: geometry + permittivity profile; 
                               vacuum wavelength ... */

class SegWgStruct 
{
public:
	// number of inner segments
	// 0: waveguide junction / end facet; -1: homogeneous waveguide
	int nz;
	// positions of the segment boundaries
	Dvector hz;

	// initialize
	SegWgStruct(); 
	SegWgStruct(int vnz); 
	SegWgStruct(Circuit rc); 
	SegWgStruct(Circuit rc, int opt); // opt = 1: merge equal segments & layers, default
	// destroy
	~SegWgStruct(); 
	// copy constructor
	SegWgStruct(const SegWgStruct& swg);
	// assignment
	SegWgStruct& operator=(const SegWgStruct& swq);

	// free allocated memory
	void clear();

	// set all segments to wg
	void init(Waveguide wg);

	// subscripting
	inline Waveguide& operator() (int s) { return wgvec[s]; }
	inline Waveguide  operator() (int s) const { return wgvec[s]; }

	// get segment index corresponding to position z
	int segidx(double z) const;
	// get segment boundaries corresponding to index s  
	Interval segment(int s) const;
	// get segment boundaries corresponding to position z  
	Interval segment(double z) const;

	// get refractive index at position x, z
	double n(double x, double z) const;

	// get permittivity at position x, z
	double eps(double x, double z) const;

	// vacuum wavelength: should be unique in all segments 
	double lambda();

	// check consistency of the SegWgStruct object: error&exit, if not ok
	void consistency() const;

	// translate: wg(s).hx -> wg(s).hx+dx for all segments
	void xtranslate(double dx);
	// translate: hz -> hz+dz
	void ztranslate(double dz);

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

	// for the QUEP simulation framework:
	// construct the corresponding vertical segmented waveguide sequence 
	SegWgStruct rotate() const;
	// enlarge the sequence by identical input and output segments
	SegWgStruct expand(Interval w) const;
	// remove the additional outer segments 
	SegWgStruct shrink() const;
	// combine equal adjacent segments 
	void optimize();

	// a default window around the "inner" structure, 
	// e.g. for purposes of plotting 
	void defaultwindow(Interval& wx, Interval& wz) const;

	// refractive index profile -> MATLAB .m file  
	// ext0, ext1: filename id chars 
	// xbeg, xend, zbeg, zend: display window in the x-z-plane
	void plot(double xbeg, double xend, double zbeg, double zend,
	          char ext0, char ext1) const;
	void plot(char ext0, char ext1) const;

	// convert this to a Circuit object
	Circuit circuit() const;
	
	// modelling of periodically corrugated waveguides:
	// extend this periodically over the display window zbeg, zend
	SegWgStruct unwrap(double zbeg, double zend) const;
	// remove outermost equal slices
	SegWgStruct reduce() const;

	// special() != 0: n stores permittivity values in place of refractive
	// indices, consistency checks for ref.ind. and perm. are disabled;
	// a misuse of the Circuit container for problems with negative 
	// permittivity
	int special() const;

private:
	Waveguide *wgvec;
};

/* vertical interfaces of a segment */
enum SBorder {SBLEFT, SBRIGHT};

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

/* for the VEIMS solver:
   rectangular section on the waveguide cross section plane */
class Rect 
{
public:
	// lower left point
	double x0;
	double y0;
	// upper right point
	double x1;
	double y1;
	// initialize
	Rect();
	Rect(double xa, double ya, double xb, double yb);
	// output
	void write(FILE *dat);
	// input
	void read(FILE *dat);
};
extern Rect XYplane;

/* for the VEIMS solver:
   cross section (2D) of a 3D rectangular waveguide, 
   geometry & permittivity profile, vacuum wavelength;
   (so far) meant as a help to define a SegWgCrs */
class WgCrs 
{
public:
	// number of inner horizontal layers
	int nx;
	// positions (x) of layer interfaces 
	Dvector hx;
	// number of inner vertical slices
	int ny;
	// positions (y) of slice interfaces 
	Dvector hy;
	
	// vacuum wavelength
	double lambda;
	// refractive indices 
	Dmatrix n;
	// permittivity on rectangle l, m
	double eps(int l, int m) const;
	// permittivity at position x, y
	double eps(double x, double y) const;

	// translate: hx -> hx+dx, hy -> hy+dy
	void translate(double dx, double dy);

	// initialize
	WgCrs(); 
	WgCrs(int vnx, int vny); 
	WgCrs(int vnx, Dvector vhx,
	      int vny, Dvector vhy,
	      double l, Dmatrix n);

	// free allocated memory
	void strip();

	// get rectangle index (l, m) corresponding to position (x,y)  
	void rectidx(double x, double y, int& l, int& m) const;
	// get rectangle boundaries corresponding to index (l,m)  
	Rect rectbounds(int l, int m) const;
	// get rectangle boundaries corresponding to position (x,y)  
	Rect rectbounds(double x, double y) const;
	// test neighbourhood of two rectangles
	int testconnect(int l0, int m0, int l1, int m1) const;

	// test matching of boundary positions with another waveguide
	int bdmatch(WgCrs w) const; 

	// refractive index profile -> MATLAB .m file  
	// ext0, ext1: filename id chars 
	// xbeg, xend, ybeg, yend: display window in the x-y-plane
	void plot(double xbeg, double xend, double ybeg, double yend,
	          char ext0, char ext1) const;
	void plot(Rect disp, char ext0, char ext1) const;
	void plot(char ext0, char ext1) const;
};

/* for the VEIMS solver:
   cross section (2D) of a 3D rectangular waveguide, 
   defined as a sequence of (1D) Waveguides */

class SegWgCrs
{
public:
	// number of inner segments, ny >=1
	int ny;
	// positions of the segment boundaries
	Dvector hy;

	// initialize
	SegWgCrs(); 
	SegWgCrs(int vny); 
	SegWgCrs(WgCrs cr); 
	// destroy
	~SegWgCrs(); 
	// copy constructor
	SegWgCrs(const SegWgCrs& swc);
	// assignment
	SegWgCrs& operator=(const SegWgCrs& swc);

	// free allocated memory
	void clear();

	// set all segments to wg
	void init(Waveguide wg);

	// subscripting
	inline Waveguide& operator() (int s) { return wgvec[s]; }
	inline Waveguide  operator() (int s) const { return wgvec[s]; }

	// input from FILE dat 
	void read(FILE *dat);

	// output to FILE dat
	void write(FILE *dat) const;

	// get segment index corresponding to position y
	int segidx(double y) const;
	// get segment boundaries corresponding to index s  
	Interval segment(int s) const;
	// get segment boundaries corresponding to position y  
	Interval segment(double y) const;

	// get refractive index at position x, y
	double n(double x, double y) const;

	// get permittivity at position x, y
	double eps(double x, double y) const;

	// vacuum wavelength: should be unique in all segments 
	double lambda();

	// check consistency of the SegWgCrs object: error&exit, if not ok
	void consistency() const;

	// translate: wg(s).hx -> wg(s).hx+dx for all segments
	void xtranslate(double dx);
	// translate: hy -> hy+dy
	void ytranslate(double dy);

	// construct the corresponding vertical segmented waveguide sequence 
	SegWgCrs rotate() const;

	// a default window around the "inner" structure, 
	// e.g. for purposes of plotting 
	void defaultwindow(Interval& wx, Interval& wy) const;

	// refractive index profile -> MATLAB .m file  
	// ext0, ext1: filename id chars 
	// xbeg, xend, ybeg, yend: display window in the x-y-plane
	void plot(double xbeg, double xend, double ybeg, double yend,
	          char ext0, char ext1) const;
	void plot(Rect disp, char ext0, char ext1) const;
	void plot(char ext0, char ext1) const;

	// convert this to a WgCrs object
	WgCrs wgcrs() const;

private:
	Waveguide *wgvec;
};

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

/* more general specification of a dielectric structure
   as a stack of possibly overlapping geometric shapes with constant 
   permittivity: */

/* individual shapes */
enum OvlType {HLAY, VLAY, RING, DISK, RECT};
class Overlay 
{
public:

	// type of the shape
	OvlType type;
	
	// initialize:
	// generic
	Overlay();
	// horizontal (HLAY) or vertical layer (VLAY), 
	// refr.ind. nc, thickness d, center at x=ofs (HLAY) or z=ofs (VLAY)
	Overlay(OvlType t, double nc, double d, double ofs);
	// RING cavity, refr.ind. nc, radius r, core width d, origin at ox, oz
	Overlay(OvlType t, double nc, double r, double d, double ox, double oz);
	// DISK cavity, refr.ind. nc, radius r, origin at ox, oz
	Overlay(OvlType t, double nc, double r, double ox, double oz);
	// RECTangular region def. by Intervals rx, rz, refr.ind. nc
	Overlay(OvlType t, double nc, Interval rx, Interval rz);

	// the interior permittivity 
	double eps;
	// and refractive index
	double n;

	// decide whether a position falls inside the shape (1) or not (0) 
	int inside(double x, double z);

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

	/// the internal representation
	double xo, zo;
	double ubd, lbd;
};

/* the entire stack */
class OvlStruct  
{
public:
	// number of overlays
	// 0: homogeneous permittivity
	int no;

	// background-permittivity
	double bgeps;
	// and refractive index
	double bgn;

	// initialize 
	OvlStruct(); 
	// initialize with background-permittivity b, wavelength l
	OvlStruct(double b, double l); 
	// destroy
	~OvlStruct(); 
	// copy constructor
	OvlStruct(const OvlStruct& swg);
	// assignment
	OvlStruct& operator=(const OvlStruct& swq);

	// free allocated memory
	void clear();

	// place an Overlay on the top of the stack
	void place(Overlay o);

	// subscripting
	inline Overlay& operator() (int s) { return ovlvec[s]; }
	inline Overlay  operator() (int s) const { return ovlvec[s]; }

	// get refractive index at position x, z
	double n(double x, double z) const;

	// get permittivity at position x, z
	double eps(double x, double z) const;

	// vacuum wavelength
	double lambda;

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

	// for purely rectangular structure: 
	// convert this to a SegWgStruct ...
	SegWgStruct segwgstruct() const;
	// ... or to a Circuit object; 
	Circuit circuit() const;

	// a default window around the "inner" structure, 
	// e.g. for purposes of plotting 
	void defaultwindow(Interval& wx, Interval& wz) const;

	// refractive index profile -> MATLAB .m file  
	// ext0, ext1: filename id chars 
	// xbeg, xend, zbeg, zend: display window in the x-z-plane
	void plot(double xbeg, double xend, double zbeg, double zend,
	          char ext0, char ext1) const;
	void plot(char ext0, char ext1) const;

private:
	Overlay *ovlvec;
};

/* ----------------------------------------------------------------------- */
#endif  // STRUCTURE_H