%
% MirrorEigenModes.m
%
% Model exported on Mar 25 2011, 12:41 by COMSOL 4.1.0.88.

%% Import
import com.comsol.model.*
import com.comsol.model.util.*

%% Create Model

model = ModelUtil.create('Model');

%model.modelPath('/home/aso/exp/LCGT/MIF/IFO/OptConf/PI/COMSOL');

model.name('TM_Eigen_Modes_Sapphire.mph');

% === Define parameters ===
model.param.set('r', '12.5[cm]', 'Radius of the mirror');
model.param.set('d', '15[cm]', 'Thickness of the mirror');
model.param.set('wl', '1064[nm]', 'Wavelength');
model.param.set('k0', '2*pi/wl', 'Wave number');
model.param.set('zr', '3464.1[m]', 'Rayleigh Range');
model.param.set('q', 'zr*i+3000', 'q-parameter at the mirror surface');

model.modelNode.create('mod1');

% === Define Hermite Polynoimnals ===
model.func.create('an1', 'Analytic');
model.func.create('an2', 'Analytic');
model.func.create('an3', 'Analytic');
model.func.create('an4', 'Analytic');
model.func.create('an5', 'Analytic');
model.func.create('an6', 'Analytic');
model.func.create('an7', 'Analytic');
model.func.create('an8', 'Analytic');
model.func.create('an9', 'Analytic');

model.func('an1').name('Hermite 0');
model.func('an1').set('funcname', 'H0');
model.func('an1').set('expr', '1');
model.func('an1').set('plotargs', {'x' '0' '10'});
model.func('an2').name('Hermite 1');
model.func('an2').set('funcname', 'H1');
model.func('an2').set('expr', '2*x');
model.func('an2').set('plotargs', {'x' '0' '10'});
model.func('an3').name('Hermite 2');
model.func('an3').set('funcname', 'H2');
model.func('an3').set('expr', '4*x^2-2');
model.func('an4').name('Hermite 3');
model.func('an4').set('funcname', 'H3');
model.func('an4').set('expr', '8*x^3-12*x');
model.func('an5').name('Hermite 4');
model.func('an5').set('funcname', 'H4');
model.func('an5').set('expr', '16*x^4-48*x^2+12');
model.func('an6').name('Hermite 5');
model.func('an6').set('funcname', 'H5');
model.func('an6').set('expr', '32*x^5-160*x^3+120*x');
model.func('an7').name('Hermite 6');
model.func('an7').set('funcname', 'H6');
model.func('an7').set('expr', '64*x^6-480*x^4+720*x^2-120');
model.func('an8').name('Hermite 7');
model.func('an8').set('funcname', 'H7');
model.func('an8').set('expr', '128*x^7-1344*x^5+3360*x^3-1680*x');
model.func('an9').name('Hermite 8');
model.func('an9').set('funcname', 'H8');
model.func('an9').set('expr', '256*x^8-3584*x^6+13440*x^4-13440*x^2+1680');


% === Variables ===
model.variable.create('var1');
model.variable('var1').name('Variables 1a');
% Common Gaussian part
model.variable('var1').set('K', 'exp(-i*k0/q*(X^2+Y^2)/2)');

% Coordinates conversion: X -> xi, Y -> zeta
model.variable('var1').set('xi', 'sqrt(k0*zr)/abs(q)*X');
model.variable('var1').set('zeta', 'sqrt(k0*zr)/abs(q)*Y');
% TEM00 Mode
model.variable('var1').set('E00', 'H0(xi)*H0(zeta)*K');

% === Geometry ===
model.geom.create('geom1', 3);
model.geom('geom1').feature.create('cyl1', 'Cylinder');
model.geom('geom1').feature('cyl1').set('r', 'r');
model.geom('geom1').feature('cyl1').set('h', 'd');
model.geom('geom1').run;

% === Material ===
model.material.create('mat1');
model.material('mat1').materialModel.create('RefractiveIndex', 'Refractive index');
model.material('mat1').materialModel.create('Anisotropic', 'Anisotropic');

model.material('mat1').name('Sapphire');
model.material('mat1').materialModel('def').set('relpermeability', {'1' '0' '0' '0' '1' '0' '0' '0' '1'});
model.material('mat1').materialModel('def').set('electricconductivity', {'1e-14[S/m]' '0' '0' '0' '1e-14[S/m]' '0' '0' '0' '1e-14[S/m]'});
model.material('mat1').materialModel('def').set('thermalexpansioncoefficient', {'0.55e-6[1/K]' '0' '0' '0' '0.55e-6[1/K]' '0' '0' '0' '0.55e-6[1/K]'});
model.material('mat1').materialModel('def').set('heatcapacity', '703[J/(kg*K)]');
model.material('mat1').materialModel('def').set('relpermittivity', {'2.09' '0' '0' '0' '2.09' '0' '0' '0' '2.09'});
model.material('mat1').materialModel('def').set('density', '3970[kg/m^3]');
model.material('mat1').materialModel('def').set('thermalconductivity', {'1.38[W/(m*K)]' '0' '0' '0' '1.38[W/(m*K)]' '0' '0' '0' '1.38[W/(m*K)]'});
model.material('mat1').materialModel('RefractiveIndex').set('n', {'1.45' '0' '0' '0' '1.45' '0' '0' '0' '1.45'});
model.material('mat1').materialModel('Anisotropic').set('D', {'4.97300000e+11' '1.62800000e+11' '1.16000000e+11' '0.00000000e+00' '-2.19000000e+10' '0.00000000' '1.62800000e+11' '4.97300000e+11' '1.16000000e+11' '0.00000000e+00' '2.19000000e+10' '0.00000000e+00' '1.16000000e+11' '1.16000000e+11' '5.00900000e+11' '0.00000000e+00' '0.00000000e+00' '0.00000000e+00' '0.00000000e+00' '0.00000000e+00' '0.00000000e+00' '1.67250000e+11' '0.00000000e+00' '-2.19000000e+10' '-2.19000000e+10' '2.19000000e+10' '0.00000000e+00' '0.00000000e+00' '1.46800000e+11' '0.00000000e+00' '0.00000000e+00' '0.00000000e+00' '0.00000000e+00' '-2.19000000e+10' '0.00000000e+00' '1.46800000e+11'});

% === Physics ===
model.physics.create('solid', 'SolidMechanics', 'geom1');
model.physics('solid').prop('AddMixedFormPressure').set('AddMixedFormPressure', '1');
model.physics('solid').feature('lemm1').set('SolidModel', 'Anisotropic');

% === Mesh ===
model.mesh.create('mesh1', 'geom1');
model.mesh('mesh1').feature('size').set('hauto', 2);
%Mapped mesh on the HR surface
model.mesh('mesh1').feature.create('map1', 'Map');
model.mesh('mesh1').feature('map1').selection.set([4]);
model.mesh('mesh1').feature('map1').feature.create('size1', 'Size');
model.mesh('mesh1').feature('map1').feature('size1').set('hauto', 2);

%Sweep the mapped mesh
model.mesh('mesh1').feature.create('swe1', 'Sweep');
model.mesh('mesh1').feature('swe1').feature.create('size1', 'Size');
model.mesh('mesh1').feature('swe1').feature('size1').set('hauto', 2);

% Generate mesh
model.mesh('mesh1').run;

%Get number of elements
nMeshElems = model.mesh('mesh1').stat.getNumElem;

%% Study

model.study.create('std1');
model.study('std1').feature.create('eig', 'Eigenfrequency');

model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').attach('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature.create('e1', 'Eigenvalue');
model.sol('sol1').feature('e1').feature.create('adDef', 'Adaption');

model.study('std1').feature('eig').set('neigs', '200');
model.study('std1').feature('eig').set('shift', '50e3');

model.sol('sol1').feature('st1').name('Compile Equations: Eigenfrequency');
model.sol('sol1').feature('st1').set('studystep', 'eig');
model.sol('sol1').feature('e1').set('neigs', '80');
model.sol('sol1').feature('e1').set('transform', 'eigenfrequency');
model.sol('sol1').feature('e1').set('shift', '50e3');
model.sol('sol1').feature('e1').set('control', 'eig');
model.sol('sol1').feature('e1').set('eigfunscale', 'mass');
model.sol('sol1').feature('e1').feature('adDef').set('adapsol', 'none');

%% Run
tic
model.sol('sol1').runAll;
toc

model.result.dataset.create('surf1', 'Surface');
model.result.dataset('surf1').selection.set([4]);

%% Save model
model.save('TM_Elastic_Modes');

%% Load saved model
model = mphload('TM_Elastic_Modes');

%% Eigenfrequency list
model.result.numerical.create('egfreqs', 'EvalGlobal');
model.result.numerical('egfreqs').set('expr', 'freq');
model.result.numerical('egfreqs').set('unit', '1/s');
model.result.numerical('egfreqs').set('descr', 'Frequency');
efreqs=model.result.numerical('egfreqs').getReal();
%model.result.numerical('egfreqs').run()
model.result.numerical.remove('egfreqs');

%% Volume plot of elastic modes

model.result.create('pg1', 'PlotGroup3D');
model.result('pg1').set('solnum',[113]);
model.result('pg1').feature.create('surf1', 'Surface');
model.result('pg1').feature.create('arws1', 'ArrowSurface');
model.result('pg1').feature('surf1').set('expr', 'w');
model.result('pg1').feature('surf1').set('unit', 'm');
model.result('pg1').feature('surf1').set('descr', 'z displacement');
mphplot(model, 'pg1')
model.result.remove('pg1')

%% Surface plot of optical modes

m=0;
n=1;
HomExpr = sprintf('abs(H%d(xi)*H%d(zeta)*K)',m,n);

model.result.create('pg2', 'PlotGroup2D');
model.result('pg2').feature.create('surf1', 'Surface');
model.result('pg2').set('solnum', '1');
model.result('pg2').feature('surf1').set('expr', HomExpr);
model.result('pg2').feature('surf1').set('descr', HomExpr);
mphplot(model, 'pg2')
model.result.remove('pg2')

%% Calculate overlap factors

tic
%Volume of the mirror
V = mphint(model, '1', 'Edim',3, 'Solnum', 'end', 'Intvolume','on');
%Normalizing factor for E00
NormE00 = mphint(model, 'abs(E00)^2', 'Edim',2, 'Solnum', 'end','Intsurface','on', 'Selection',[4]);
%Volume integration of displacement for normalizing the elastic modes
Vdisp = mphint(model, 'abs(solid.disp)^2', 'Edim',3, 'Intvolume','on');

%Maximum order of the optical HOMs.
idxMax = 8;

%Number of elastic modes
sz = mphgetp(model);
nElasticModes = sz(2);

%Array to hold overlap factor
%Lambda(m+1,n+1,p) returns the overlap factor for TEM(m,n) and p-th elastic
%mode.
Lambda = zeros(idxMax+1, idxMax+1, nElasticModes);


for m=0:idxMax
    for n=0:idxMax
        OvlpExpr = sprintf('E00*H%d(xi)*H%d(zeta)*K*w',m,n);
        HomExpr = sprintf('abs(H%d(xi)*H%d(zeta)*K)^2',m,n);
        
        %Surface integral of overlap
        Sovlp = mphint(model, OvlpExpr, 'Edim',2, 'Intsurface','on', 'Selection',[4]);
        %Normalization factor for HOM
        NormHOM = mphint(model, HomExpr, 'Edim',2, 'Solnum', 'end','Intsurface','on', 'Selection',[4]);
        
        Lambda(m+1,n+1,:) = V*abs(Sovlp).^2./(NormE00*NormHOM*Vdisp);
    end
end
toc


%% Save overlap factors
save('OverlapFactorsETM.mat','Lambda','efreqs');


%% Integration
OvlpExpr = 'E00*H1(xi)*H2(zeta)*K*w';
HomExpr = 'abs(H1(xi)*H2(zeta)*K)^2';
%Volume of the mirror
V = mphint(model, '1', 'Edim',3, 'Solnum', 'end', 'Intvolume','on');

Sovlp = mphint(model, OvlpExpr, 'Edim',2, 'Intsurface','on', 'Selection',[4]);
SE00 = mphint(model, 'abs(E00)^2', 'Edim',2, 'Solnum', 'end','Intsurface','on', 'Selection',[4]);
Shom = mphint(model, HomExpr, 'Edim',2, 'Solnum', 'end','Intsurface','on', 'Selection',[4]);
Vdisp = mphint(model, 'abs(solid.disp)^2', 'Edim',3, 'Intvolume','on');

OvlpFact = V*abs(Sovlp).^2./(SE00*Shom*Vdisp);

%