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

model = ModelUtil.create('Model');
comp1 = model.component.create('comp1',true);
geom3 = comp1.geom.create('geom3', 3); % 3D geometry

% --- Parameters ---
model.param.set('wg_length', '20[um]');
model.param.set('wg_width', '600[nm]');
model.param.set('wg_height', '300[nm]');
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');

% Parameters for taper cavity holes

model.param.set('cavity_lattice','range(at_lattice,taper_period,a0_lattice)');
model.param.set('cavity_r_par','range(rt_parallel,taper_period,r0_parallel)');
model.param.set('cavity_r_perp','range(rt_perpendicular,taper_period,r0_perpendicular)');

% Evaluate the vector parameters as MATLAB arrays
cavity_lattice = model.param.evaluate('cavity_lattice');
cavity_r_par = model.param.evaluate('cavity_r_par');
cavity_r_perp = model.param.evaluate('cavity_r_perp');

% --- Work Plane in 3D geometry ---
wp1 = geom3.feature.create('wp1', 'WorkPlane');
wp1.set('quickplane', 'xy');
wp_geom = wp1.geom;

% Nanobeam rectangle in work plane
beam = wp_geom.feature.create('beam', 'Rectangle');
beam.set('size', {'wg_length', 'wg_width'});
% Center the nanobeam at (0,0)
beam.set('pos', {'-wg_length/2', '-wg_width/2'});

% Calculate the center offset for symmetry
center_offset = sum(cavity_lattice)/2;

x_pos = 0;
for i = 1:numel(cavity_lattice)
    r_par_val = cavity_r_par(i);
    r_perp_val = cavity_r_perp(i);

    % Position along x, centered
    x_pos = sum(cavity_lattice(1:i)) - center_offset;

    % Create ellipse at +x
    ellipse_p = wp_geom.feature.create(['ellipse_px' num2str(i)], 'Ellipse');
    ellipse_p.set('pos', {num2str(x_pos), '0'});
    ellipse_p.set('semiaxes', {num2str(r_par_val), num2str(r_perp_val)});

    % Optionally, if you want symmetry, you can also create at -x (mirror)
    ellipse_m = wp_geom.feature.create(['ellipse_mx' num2str(i)], 'Ellipse');
    ellipse_m.set('pos', {num2str(-x_pos), '0'});
    ellipse_m.set('semiaxes', {num2str(r_par_val), num2str(r_perp_val)});
end
% Build the work plane geometry
wp_geom.run;

% Build the 3D geometry (extrude the work plane)
extr1 = geom3.feature.create('extr1', 'Extrude');
extr1.set('distance', 'wg_height');
extr1.selection('input').set({'wp1'});

geom3.run;
% mphgeom(model,'geom3');
mphlaunch(model);
% Launch COMSOL's built-in geometry viewer
% ModelUtil.show(model);