#ifndef MAXWEQ_H #define MAXWEQ_H /* * WMM - Wave-matching method for mode analysis of dielectric waveguides * https://wmm.computational-photonics.eu/ */ /* * maxweq.h * relations concerning the electromagnetic fields for waveguide modes */ /* 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; }; /* field components */ enum Fcomp {EX, EY, EZ, HX, HY, HZ, SZ}; /* Fcomp -> field character */ char fldchr(Fcomp cp); /* Fcomp -> comp character */ char cpchr(Fcomp cp); /* 'EHS''xyz' -> Fcomp */ Fcomp fcomp(char eh, char cp); /* polarization parameter */ enum Polarization {QTE, QTM, VEC}; /* Polarization -> 't', 'v' */ char poltochar1(Polarization p); /* Polarization -> 'e', 'm', 'c' */ char poltochar2(Polarization p); /* Polarization -> principal field component */ Fcomp principalcomp(Polarization p); /* symmetry parameter */ enum Symmetry {SYM, ASY, NOS}; /* Symmetry -> 'n', 's', 'a' */ char symtochar(Symmetry s); /* 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}; /* derivation levels */ enum Derlev {FLD=0, DXF=1, DYF=2, DXY=3, DXX=4, DYY=5}; /* vectorial mode formulation */ enum Vecform {EXEY, HXHY, EZHZ}; /* ##################################################################### */ /* representation of a field component as a sum of derivations of one or two basic fields, depending on the polarization and formulation of the vectorial mode problem */ /* variable factors */ enum Vfac {SFID, SFBETA, SFINVBETA, SFEPS, SFIBMK0E, SFK0EPS, SFK0, SFINVOMEPSEP0, SFINVOMMU0}; /* ... evaluate */ double evalvfac(Vfac f, double b, double k0, double n); /* the global representations ... */ typedef struct { int numc; /* Number of contributions; 1, 2, or 3 */ Vfac glsfid; /* global factors: variabel */ double glfac; /* fixed */ struct { int c; /* basic component 0/1 */ Derlev fod; /* its derivation level */ Vfac sfid; /* variable factor */ double fac; /* fixed factor */ } ctr[3]; } Fieldrep; extern Fieldrep qteEHrep[6]; // QTE extern Fieldrep qtmEHrep[6]; // QTM extern Fieldrep vefEHrep[6]; // vectorial EXEY formulation extern Fieldrep vhfEHrep[6]; // vectorial HXHY formulation extern Fieldrep vzfEHrep[6]; // vectorial EZHZ formulation /* Fcomp fc is mapped to EHrep[iFcomp(fc)] */ int iFcomp(Fcomp fc); /* Fcomp fc is mapped to a Fieldrep, depending on polarization an vec. form. */ Fieldrep getfieldrep(Fcomp cp, Polarization p, Vecform v); /* initialize global field representations */ void initfieldreps(); class Maxweq { public: // initialize Maxweq(); }; extern Maxweq MINIT; /* ##################################################################### */ /* typical initial field: linearly polarized Gaussian beam with circular profile, propagating in a homogeneous medium in z-direction */ class Gaussianbeam { public: double l; // vacuum wavelength, given in um double ng; // refractive index double x0; // location of intensity maximum double y0; double rg; // beam radius double p0; // total power in W double alpha; // polarization angle versus TE // initialize Gaussianbeam(); Gaussianbeam(double l, double ng, double x0, double y0, double rg, double p0, double alpha); // real fields in A/um or V/um double field(Fcomp cp, double x, double y); }; /* overlap 0.5*intint fld1_EX*fld2_HX - fld1_EY*fld2_HY dxdy of real fields, integrals are evaluated as sums over step-functions on an equidistant mesh with (numx+1)*(numy+1) points on rectangle r */ double overlap(double (*fld1)(Fcomp cp, double x, double y), double (*fld2)(Fcomp cp, double x, double y), Rect r, int numpx, int numpy); /* ##################################################################### */ /* permittivity perturbation */ class Perturbation { public: // part of x-y-plane with nonzero perturbation Rect r; // entries of the permittivity perturbation tensor, complex 3x3 matrix Cmatrix e; // initialize Perturbation(); Perturbation(Rect r); Perturbation(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); Perturbation(Cmatrix p); // relevant quantities for isotropic absorbing part // and Hermitian anisotropic part Perturbation(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) */ Perturbation 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] */ Perturbation 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, */ Perturbation magopt_polar(double sFrot, double n, double lambda, Rect rl); /* for equatorial orientation of the magnetization, vec(M) || y-axis, */ Perturbation magopt_equat(double sFrot, double n, double lambda, Rect rl); /* and for longitudinal orientation of the magnetization, vec(M) || z-axis */ Perturbation magopt_long(double sFrot, double n, double lambda, Rect rl); /* ##################################################################### */ #endif // MAXWEQ_H