% PCF_Structure.m - Photonic Crystal Fiber Structure Visualization
clear all; close all;

% Parameters (converted to meters)
d = 2.55e-6;         % Diameter of air holes
r = d / 2;           % Radius of air holes
pitch = 3e-6;        % Pitch (can be changed to 2.7, 2.8, 2.9, or 3.0 μm)
r_clad = 15e-6;      % Cladding radius (changed from thickness to radius)
PML = 2e-6;          % Perfectly Matched Layer thickness (changed to 2μm as requested)
r_total = r_clad + PML; % Total radius including PML

% For annotation
pitch1e6 = pitch * 1e6;
d1e6 = d * 1e6;
r1e6 = r * 1e6;
r_clad1e6 = r_clad * 1e6;
PML1e6 = PML * 1e6;
r_total1e6 = r_total * 1e6;

% Create figure
figure('Position', [100, 100, 800, 800]);
hold on; axis equal;

% Plot PML layer (outermost, dark gray)
theta_pml = linspace(0, 2*pi, 200);
x_pml = r_total1e6 * cos(theta_pml);
y_pml = r_total1e6 * sin(theta_pml);
fill(x_pml, y_pml, [0.3 0.3 0.3], 'EdgeColor', 'none');

% Plot cladding (silica, light gray)
theta_clad = linspace(0, 2*pi, 200);
x_clad = r_clad1e6 * cos(theta_clad);
y_clad = r_clad1e6 * sin(theta_clad);
fill(x_clad, y_clad, [0.85 0.85 0.85], 'EdgeColor', 'k', 'LineWidth', 1);

% Get lattice positions
[x, y] = hex_lattice(pitch, 3);

% Sort holes by distance from center for proper layering
distances = sqrt(x.^2 + y.^2);
[~, idx] = sort(distances);
x = x(idx); y = y(idx);

% Plot holes (only those within the cladding)
for i = 1:length(x)
    dist = sqrt(x(i)^2 + y(i)^2);
    
    % Check if hole center is within cladding
    if dist + r < r_clad
        % First layer (semi-circular holes)
        if abs(dist - pitch) < pitch * 0.1 % Tolerance for first ring
            % Calculate angle to center
            theta = atan2(y(i), x(i));
            % Semi-circle with flat side towards center
            plot_semicircle(x(i), y(i), r, theta + pi/2, theta + 3*pi/2);
        else
            % Other layers (circular holes)
            plot_circle(x(i), y(i), r);
        end
    end
end

% Mark the core region
plot(0, 0, 'r+', 'MarkerSize', 15, 'LineWidth', 3);
text(0.3, 0.3, 'Core', 'FontSize', 14, 'Color', 'r', 'FontWeight', 'bold');

% Add scale bar
plot([-5, -3], [-16, -16], 'k-', 'LineWidth', 3);
text(-4, -16.5, '2 μm', 'FontSize', 12, 'HorizontalAlignment', 'center');

% Add annotations for different regions
% First ring annotation
theta_ann = pi/6;
x_ann = pitch * cos(theta_ann) * 1e6;
y_ann = pitch * sin(theta_ann) * 1e6;
plot([0, x_ann], [0, y_ann], 'b--', 'LineWidth', 1);
text(x_ann/2, y_ann/2 + 0.5, 'First ring', 'FontSize', 10, 'Color', 'b');

% Add pitch annotation
plot([0, pitch1e6], [-12, -12], 'k-', 'LineWidth', 1);
plot([0, 0], [-11.8, -12.2], 'k-', 'LineWidth', 1);
plot([pitch1e6, pitch1e6], [-11.8, -12.2], 'k-', 'LineWidth', 1);
text(pitch1e6/2, -13, sprintf('Pitch = %.1f μm', pitch1e6), ...
    'FontSize', 12, 'HorizontalAlignment', 'center');

% Add cladding radius annotation
plot([0, r_clad1e6], [0, 0], 'k--', 'LineWidth', 1);
text(r_clad1e6/2, 1, sprintf('r_{clad} = %.0f μm', r_clad1e6), ...
    'FontSize', 10, 'HorizontalAlignment', 'center');

% Add PML annotation
theta_pml_ann = 3*pi/4;
x_pml_start = r_clad1e6 * cos(theta_pml_ann);
y_pml_start = r_clad1e6 * sin(theta_pml_ann);
x_pml_end = r_total1e6 * cos(theta_pml_ann);
y_pml_end = r_total1e6 * sin(theta_pml_ann);
plot([x_pml_start, x_pml_end], [y_pml_start, y_pml_end], 'm-', 'LineWidth', 2);
text(x_pml_end - 2, y_pml_end + 1, sprintf('PML = %.0f μm', PML1e6), ...
    'FontSize', 10, 'Color', 'm', 'HorizontalAlignment', 'center');

% Add hole diameter annotation (find a visible hole)
for i = 1:length(x)
    if x(i) > 0 && y(i) > 0 && sqrt(x(i)^2 + y(i)^2) < r_clad - r
        x_hole = x(i) * 1e6;
        y_hole = y(i) * 1e6;
        plot([x_hole - r1e6, x_hole + r1e6], [y_hole, y_hole], 'g-', 'LineWidth', 2);
        text(x_hole, y_hole - 1, sprintf('d = %.2f μm', d1e6), ...
            'FontSize', 10, 'Color', 'g', 'HorizontalAlignment', 'center');
        break;
    end
end

% Labels and title
xlabel('x (μm)', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('y (μm)', 'FontSize', 14, 'FontWeight', 'bold');
title(sprintf('Photonic Crystal Fiber Cross-Section\nPitch = %.1f μm, Hole diameter = %.2f μm, Cladding radius = %.0f μm', ...
    pitch1e6, d1e6, r_clad1e6), 'FontSize', 16, 'FontWeight', 'bold');

% Set axis properties
xlim([-18, 18]);
ylim([-18, 18]);
grid on; grid minor;
set(gca, 'FontSize', 12);
box on;

% Add legend
legend_items = [
    fill(nan, nan, [0.85 0.85 0.85], 'EdgeColor', 'k', 'LineWidth', 1);
    fill(nan, nan, [0.3 0.3 0.3], 'EdgeColor', 'none');
    fill(nan, nan, 'w', 'EdgeColor', 'k', 'LineWidth', 1.5);
];
legend(legend_items, {'Silica Cladding', 'PML', 'Air Holes'}, ...
    'Location', 'northeast', 'FontSize', 10);

hold off;

% Print structure parameters
fprintf('\n=== PCF Structure Parameters ===\n');
fprintf('Pitch: %.2f μm\n', pitch * 1e6);
fprintf('Hole diameter: %.2f μm\n', d * 1e6);
fprintf('Hole radius: %.2f μm\n', r * 1e6);
fprintf('Cladding radius: %.0f μm\n', r_clad * 1e6);
fprintf('PML thickness: %.0f μm\n', PML * 1e6);
fprintf('Total radius: %.0f μm\n', r_total * 1e6);
fprintf('Number of rings: 3\n');
fprintf('First ring: Semi-circular holes\n');
fprintf('Other rings: Circular holes\n');
fprintf('Core: Missing central hole\n');
fprintf('================================\n');

% --- Functions ---

function [x, y] = hex_lattice(p, layers)
    x = [];
    y = [];
    for m = -layers:layers
        for n = -layers:layers
            if abs(m + n) <= layers % Ensure hexagonal shape
                x_pos = p * (m + n/2);
                y_pos = p * sqrt(3)/2 * n;
                % Exclude the central point (core)
                if ~(m == 0 && n == 0)
                    x = [x; x_pos];
                    y = [y; y_pos];
                end
            end
        end
    end
end

function plot_semicircle(cx, cy, r, theta_start, theta_end)
    theta = linspace(theta_start, theta_end, 100);
    x = cx + r * cos(theta);
    y = cy + r * sin(theta);
    % Add the flat side
    x = [x, cx + r * cos(theta_end), cx + r * cos(theta_start)];
    y = [y, cy + r * sin(theta_end), cy + r * sin(theta_start)];
    fill(x * 1e6, y * 1e6, 'w', 'EdgeColor', 'k', 'LineWidth', 1.5); % Convert to μm
end

function plot_circle(cx, cy, r)
    theta = linspace(0, 2*pi, 100);
    x = cx + r * cos(theta);
    y = cy + r * sin(theta);
    fill(x * 1e6, y * 1e6, 'w', 'EdgeColor', 'k', 'LineWidth', 1.5); % Convert to μm
end