import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');
comp1 = model.component.create('comp1',true);
geom1 = comp1.geom.create('geom1', 3);
mphlaunch;

% --- Define all parameters as COMSOL parameters ---
model.param.set('wg_height', '300[nm]');
model.param.set('wg_width', '600[nm]');
model.param.set('wg_length', '20[um]');
model.param.set('a0_lattice', '350[nm]');
model.param.set('at_lattice', '315[nm]');
model.param.set('r0_parallel', '95[nm]');
model.param.set('r0_perpendicular', '105[nm]');
model.param.set('rt_parallel', '20[nm]');
model.param.set('rt_perpendicular', '20[nm]');
model.param.set('taper_period', '15');

num_taper = model.param.evaluate('taper_period');

% --- Parameterize taper holes ---
for i = 1:num_taper
    tfac = (i-1)/(num_taper-1);
    model.param.set(['taper_lattice_' num2str(i)], ...
        ['at_lattice + (a0_lattice-at_lattice)*' num2str(tfac)]);
    model.param.set(['taper_r_par_' num2str(i)], ...
        ['rt_parallel + (r0_parallel-rt_parallel)*' num2str(tfac)]);
    model.param.set(['taper_r_perp_' num2str(i)], ...
        ['rt_perpendicular + (r0_perpendicular-rt_perpendicular)*' num2str(tfac)]);
end

% --- Calculate taper length and mirror count ---
taper_length_expr = '0';
for i = 1:num_taper
    taper_length_expr = [taper_length_expr '+taper_lattice_' num2str(i)];
end
model.param.set('taper_length_one_side', taper_length_expr);
model.param.set('available_length_for_mirrors', 'wg_length/2 - taper_length_one_side');
model.param.set('num_mirror', 'ceil(available_length_for_mirrors/a0_lattice)');
num_mirror = model.param.evaluate('num_mirror');

% --- Geometry construction using parameters ---
blk = geom1.feature.create('blk', 'Block');
blk.set('size', {'wg_length', 'wg_width', 'wg_height'});
blk.set('base', 'center');
blk.set('pos', {'0', '0', '0'});

% SiO2 layer
sio2 = geom1.feature.create('sio2', 'Block');
sio2.set('size', {'wg_length', 'wg_width', 'wg_height/2'});
sio2.set('base', 'center');
sio2.set('pos', {'0', '0', '-3*wg_height/4'});

% Substrate
sub = geom1.feature.create('sub', 'Block');
sub.set('size', {'wg_length', '20*wg_width', '5*wg_height'});
sub.set('base', 'center');
sub.set('pos', {'0', '0', '-3.5*wg_height'});

% --- Create holes using parameterized values ---
ell_names = {};

% LEFT taper holes (center to -x)
x_pos = '0';
for i = 1:num_taper
    ell = geom1.feature.create(['ell_LT_' num2str(i)], 'Ellipsoid');
    ell.set('semiaxes', {['taper_r_par_' num2str(i)], ['taper_r_perp_' num2str(i)], 'wg_height/2*1.1'});
    ell.set('pos', {x_pos, '0', '0'});
    ell_names{end+1} = ['ell_LT_' num2str(i)];
    if i < num_taper
        x_pos = [x_pos ' - taper_lattice_' num2str(i)];
    end
end

% LEFT mirror holes
for i = 1:num_mirror
    x_m = [x_pos ' - ' num2str(i) '*a0_lattice'];
    ell = geom1.feature.create(['ell_LM_' num2str(i)], 'Ellipsoid');
    ell.set('semiaxes', {'r0_parallel', 'r0_perpendicular', 'wg_height/2*1.1'});
    ell.set('pos', {x_m, '0', '0'});
    ell_names{end+1} = ['ell_LM_' num2str(i)];
end

% RIGHT taper holes (center to +x)
x_pos = '0';
for i = 2:num_taper % skip center (already created)
    x_pos = [x_pos ' + taper_lattice_' num2str(i-1)];
    ell = geom1.feature.create(['ell_RT_' num2str(i)], 'Ellipsoid');
    ell.set('semiaxes', {['taper_r_par_' num2str(i)], ['taper_r_perp_' num2str(i)], 'wg_height/2*1.1'});
    ell.set('pos', {x_pos, '0', '0'});
    ell_names{end+1} = ['ell_RT_' num2str(i)];
end

% RIGHT mirror holes
for i = 1:num_mirror
    x_m = [x_pos ' + ' num2str(i) '*a0_lattice'];
    ell = geom1.feature.create(['ell_RM_' num2str(i)], 'Ellipsoid');
    ell.set('semiaxes', {'r0_parallel', 'r0_perpendicular', 'wg_height/2*1.1'});
    ell.set('pos', {x_m, '0', '0'});
    ell_names{end+1} = ['ell_RM_' num2str(i)];
end

% Subtract all ellipsoids from the waveguide block
dif = geom1.feature.create('dif', 'Difference');
dif.selection('input').set({'blk'});
dif.selection('input2').set(ell_names);
% 
% geom1.run();
% % mphlaunch();