//Script for FreeFem++cs  http://www.ann.jussieu.fr/~lehyaric/ffcs/index.htm
//to compute WG eigen vector (magnetic field) and value (wave vector)
// of axially symmetrical dielectric resonators
// See Chaper 9 the example Whispering gallery modes of teh DOC. 


//Developed by I. S. Grudinin at JPL/Caltech
//Copyright 2012 California Institute of Technology. Government sponsorship acknowledged.
//If this script was helpful in a project that resulted in a publication, please cite our paper:
//``Finite element modeling of coupled optical microdisk resonators for displacement sensing,''
//I. S. Grudinin and Nan Yu (2012)
//    J. Opt. Soc. Am. B

//based on "Traceable 2-D finite-element simulation of the whispering-gallery modes of axisymmetric electromagnetic resonators"
//by M. Oxborrow
//IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES  Volume: 55   Issue: 6   Pages: 1209-1218

real radius=36; //approximate radius of the cavity
real yb=-10,yt=-yb; //window yb=bottom and yt=top coordinates
real xl=radius-5,xr=radius+3; //window xl=left and xr=right coordinates
real angle=asin((yt)/radius); //angle of the sphere segment to model in radians
int Nm=60;  	//number of mesh vertices per border 
real ne=1.46; //n_e-extraordinary refractive index (root of permittivity parallel to z-axis, epara)
real no=1.46; //n_o-ordinary refractive index (root of permittivity orthogonal to z-axis, eorto)
real nm=1;     //refractive index of surrounding medium (air)

int nev=4; // number of eigen values to find

int M=213;	//azimuthal mode order ~ 2Pi*n*R/lambda
real alpha=1; 	//penalty term

//describing the mesh:
border W1l(t=0,1) { x=xl+(radius*cos(angle)-xl)*(1-t); y=yt; label=1;}  // window upper border left
border W1r(t=0,1) { x=xr-(xr-radius*cos(angle))*(t); y=yt; label=1;}  	// window upper border right
border W2(t=0,1) { x=xr; y=yb+(yt-yb)*t; label=1;}  			// window right border
border W3l(t=0,1) { x=xl+(radius*cos(angle)-xl)*(t); y=yb; label=1;}  	// window lower border left
border W3r(t=0,1) { x=xr-(xr-radius*cos(angle))*(1-t); y=yb; label=1;}  // window lower border right
border W4(t=0,1) { x=xl; y=yt-(yt-yb)*t; label=1;} 			// window left border
border S(t=0,1) { x=radius*cos((t-0.5)*2*angle); y=radius*sin((t-0.5)*2*angle); label=2;}  //sphere segment
mesh Th = buildmesh(W1r(Nm/4)+W1l(Nm/4)+W4(Nm)+W3l(Nm/4)+W3r(Nm/4)+W2(Nm)+S(Nm));
plot (Th, WindowIndex = 0);


//setting subdomain parameters, i.e. refractive indices of resonator and surrounding air box
fespace Ph(Th, P0);  
Ph reg=region; // function that returns region number of a point XY
int ncav=reg(xl+1, 0);  // cavity
int nair=reg(xr-1, 0);		//air
Ph eorto=no^2*(region==ncav)+nm^2*(region==nair);  //subdomains for epsilon values inside and outside the resonators
Ph epara=ne^2*(region==ncav)+nm^2*(region==nair);  //subdomains for epsilon values inside and outside the resonators

//supplementary variables to store eigenvectors, defined on mesh Th with P2 elements - Largange quadratic.
fespace Supp(Th, P2);
Supp eHsqr;

//3d vector FE space 
fespace Vh(Th,[P2,P2,P2]);
Vh [Hr, Hphi, Hz], [vHr, vHphi, vHz]; //magnetic field components on Vh space and test functions vH

real sigma =(M/(ne*radius))^2+2; // value of the shift (k^2), where the modes will be found

//boundary condition macros
macro EWall(Hr,Hphi, Hz) (dy(Hr)-dx(Hz)+Hr*N.x+Hz*N.y-epara*(Hz*M-dy(Hphi)*x)*N.y+eorto*(Hphi-Hr*M+dx(Hphi)*x)*N.x) // eom
macro MWall(Hr,Hphi, Hz) (Hphi+Hz*N.x-Hr*N.y+epara*(Hz*M-dy(Hphi)*x)*N.x+eorto*(Hphi-Hr*M+dx(Hphi)*x)*N.y ) // eom
//b matrix (see freefem++ manual, section 9 - eigenvalue problems)
varf b([Hr, Hphi, Hz], [vHr, vHphi, vHz]) = int2d(Th)( x*(Hr*vHr+Hphi*vHphi+Hz*vHz) );
// OP = A - sigma B ; // the shifted matrix
varf op([Hr, Hphi, Hz], [vHr, vHphi, vHz])=
int2d(Th)(
(
	(eorto*(vHphi*Hphi-M*(vHphi*Hr+Hphi*vHr)+M^2*vHr*Hr)+epara*M^2*vHz*Hz)/x   //A/r
	+eorto*(dx(vHphi)*(Hphi-M*Hr)+dx(Hphi)*(vHphi-M*vHr))-epara*M*(vHz*dy(Hphi)+Hz*dy(vHphi)) //B
	+x*(eorto*dx(vHphi)*dx(Hphi)+epara*((dx(vHz)-dy(vHr))*(dx(Hz)-dy(Hr))+dy(vHphi)*dy(Hphi))) //C
)/(eorto*epara)
+alpha*(
	(vHr*Hr-M*(vHphi*Hr+Hphi*vHr)+M^2*vHphi*Hphi)/x  //D/r
	+(dx(vHr)+dy(vHz))*(Hr-M*Hphi)+(vHr-M*vHphi)*(dx(Hr)+dy(Hz)) //E
	+x*(dx(vHr)+dy(vHz))*(dx(Hr)+dy(Hz))  //F
)  
-sigma*x*(vHr*Hr+vHphi*Hphi+vHz*Hz)
)
//electric wall boundary condition on the boundary of computation domain
+int1d(Th,1)(EWall(Hr,Hphi,Hz)*EWall(vHr, vHphi, vHz));

//setting sparce matrices and assigning the solver UMFPACK to solve eigenvalue problem
matrix B=b(Vh,Vh,solver="SPARSESOLVER");
matrix OP= op(Vh,Vh,solver="SPARSESOLVER");

real[int] ev(nev); // to store the nev eigenvalue
Vh[int] [eHr,eHphi,eHz](nev); // to store the nev eigenvector
//calling ARPACK on sparce matrices with the assigned solver UMFPACK:
int k=EigenValue(OP,B,sym=true,sigma=sigma,value=ev,vector=eHr,tol=1e-10,maxit=0,ncv=0);

k=min(k,nev); //  sometimes the number of converged eigen values 
              // can be greater than nev;


//file to output mode values
ofstream f("modes.txt");
//setting number of digits in the file output
int nold=f.precision(11);


for (int i=0;i<k;i++)
{
 real lambda=2*pi/sqrt(ev[i]);
 eHsqr=(sqrt(eHr[i]^2+eHphi[i]^2+eHz[i]^2)); // intensity from magnetic field components
 plot(eHsqr,WindowIndex=i, value=1, nbiso=20,LabelColors=1, aspectratio=1, cmm="Mode "+i+", lambda="+lambda+", F="+(299792.458/lambda));
 f << "Mode "<<i<<", ka=" << sqrt(ev[i])*radius << endl;

};
