#ifndef FIM_H #define FIM_H /* * METRIC --- Mode expansion modeling in integrated optics / photonics * http://metric.computational-photonics.eu/ */ /* * fim.h * guided modes of circular multi-step-index optical fibers */ /* ------------------------------------------------------------------------ */ /* a fiber mode */ class FIMode { public: // fiber radius double R; // center position double x0; double y0; // refractive index profile; positions relative to radius R Waveguide wg; // the fiber structure OvlStruct fib; // angular order of the FIM int aol; // in case of aol ==0 : pid = 'e' (TE) or pid = 'm' (TM) // for aol >= 1: pid = 'h' (hybrid) char pid; // vacuum wavenumber double k0; // angular frequency * vacuum permittivity double omep0; // angular frequency * vacuum permeability double ommu0; // propagation constant double beta; // effective mode index double neff; // normalized effective permittivity double npcB; // FIM field Complex field(Fcomp cp, double x, double y) const; Complex field(Fcomp cp, double r) const; void emfield(double x, double y, Complex &ex, Complex &ey, Complex &ez, Complex &hx, Complex &hy, Complex &hz) const; void emfield_r(double x, double y, Complex &er, Complex &et, Complex &ez, Complex &hr, Complex &ht, Complex &hz) const; // ... values on a rectangular mesh Cmatrix field(Fcomp cp, Rect disp, int npx, int npy) const; // directional interference: mix with a*(reversed mode) Complex field_ri(Fcomp cp, double x, double y, Complex a) const; void emfield_ri(double x, double y, Complex a, Complex &er, Complex &et, Complex &ez, Complex &hr, Complex &ht, Complex &hz) const; // mode order: for each aol, an index counting (found) modes // with increasing effective index, starting with 1 int ord; // token of the mode char ids[20]; char sids[20]; // set mode id with mode count o void classify(int o); // set normalized effective permttivity B void setnpcB(); // power transported in axial direction double pfSz(double r) const; double power() const; // normalize to unit power void normalize(); // test all interface conditions, 0: reject, 1: acceptable int check(); // evaluation of confinement factors: // numerically integrate product of field components // (conjugate of first field) over full circle x radial layer l Complex lfp(Fcomp cpc, Fcomp cpn, double r) const; Complex numquad(Fcomp cpc, Fcomp cpn, int l) const; // translate the FIM: shift the origin by dx, dy void translate(double dx, double dy); // direction of angular momentum: FORW (anticlockwise, default), BACK (clockwise) Propdir dir; // reverse the FIM: clockwise propagation <-> anticlockwise propagation void reverse(); // ----------------------------------------------------------- // visualization: output to MATLAB .m files // evaluation of the field pattern in the display window // xbeg <= x <= xend, zbeg <= z <=zend // npx, npz: number of points in output mesh // ext0, ext1: filename id characters for the .m file // profile plot corresponding to component cp (static, t=0) // cp: EX, EY, EZ, HX, HY, HZ, ER, ET, HR, HT // foa: MOD, SQR, REP, IMP void plot(Fcomp cp, Afo foa, Rect disp, int npx, int npy, char ext0, char ext1) const; // export full field data ("everything") into a viewer MATLAB file void viewer(Rect disp, int npx, int npy, char ext0, char ext1) const; // piecewise internal representation by Bessel functions, // coefficients Dvector cA; Dvector cB; Dvector cC; Dvector cD; // radial wavenumber, layer-wise Dvector kappa; // evaluate transverse resonance condition for effective index n, // f=1: set cA, cB, cC, cD accordingly double travres(double n, int f); double tr_error; // remaining error // trap a root in the transverse resonance condition (bisection) double converge(double n0, double n1); // radial profile & derivative, principal components void profile(double r, double &e, double &de, double &h, double &dh) const; // interior and exterior crop radii <-> mode normalization double cRi; double cRe; // use discretized versions (1) of the radial profile int disc_prof; // on radial interval [ir < wg.hx(.) < re] with nr steps void discretize(double ri, double re, int nr); private: // storage space for discretized profiles Ivector n_rad; Dmatrix v_rad; Dmatrix v_prfE, v_prfH; Dmatrix v_dprE, v_dprH; Dmatrix d_prfE, d_prfH; Dmatrix d_dprE, d_dprH; double d_ri; double d_re; // transverse resonance evaluation double travres(); double travresTE(); double travresTM(); double cBe, cDe; int fillc; }; /* ------------------------------------------------------------------------ */ /* an array of FIMs */ class FIModeArray { public: // number of modes included int num; // initialize FIModeArray(); // destroy ~FIModeArray(); // copy constructor FIModeArray(const FIModeArray& ma); // assignment FIModeArray& operator=(const FIModeArray& s); // free allocated memory void clear(); // subscripting FIMode& operator() (int i); FIMode operator() (int i) const; // add a mode void add(FIMode m); // delete a mode entry void remove(int i); // add an entire FIModeArray nma void merge(FIModeArray& nma); // sort the array by FIM attenuation: lowest first void sort(); // prune the array: remove modes with identical orders l,m, void prune(); private: FIMode *mvec; }; /* ------------------------------------------------------------------------ FIM mode analysis procedures ------------------------------------------------------------------------ */ /* solver procedures */ /* small radius/wavelength, to avoid singluarity at the origin */ extern double FIMp_R_TINY; /* minimum effective index distance from local refractive indices */ extern double FIMsP_MinNDiff; /* convergence to TRC roots: effective index accuray */ extern double FIMsP_NeffAcc; /* minimum effective index distance for modes of the same angular order to be considered different */ extern double FIMsP_NeffSep; /* for a normalized mode: acceptable violation of interface conditions */ extern double FIMsP_IFCerr; /* for the fiber with radial layering w, radius r, centered at (xc, yc), find guided FIMs at vacuum wavelength w.lambda quiet = 0, 1, 2: all, less, no log output */ /* ... modes of all angular mode orders */ int fimsolve(Waveguide w, double r, double xc, double yc, FIModeArray& fima, int quiet); /* ... modes with angular orders between and including lmin and lmax */ int fimsolve(Waveguide w, double r, double xc, double yc, int lmin, int lmax, FIModeArray& fima, int quiet); /* internal: identify guided FIMs of angular order l | l==0: TE (p='e') or TM (p='m') */ int fimsolveproc(Waveguide w, double r, double xc, double yc, int l, char p, FIModeArray& fima, int quiet); /* Bessel functions & derivatives, integer order, real argument, wrapper functions */ double BfJ(int o, double a); double BfY(int o, double a); double BfK(int o, double a); double BfI(int o, double a); double dBfJ(int o, double a); double dBfY(int o, double a); double dBfK(int o, double a); double dBfI(int o, double a); /* error: FIMCylEvalError != 0 */ extern int FIMCylEvalError; #endif // FIM_H