forked from lydia1895/PMM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
main_ellipse_Huygens.m
executable file
·123 lines (94 loc) · 2.59 KB
/
main_ellipse_Huygens.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
clc
clear all
% Find all windows of type figure, which have an empty FileName attribute.
allPlots = findall(0, 'Type', 'figure', 'FileName', []);
% Close.
delete(allPlots);
verbose = 6;
format long
eta=0;
f1=0;
%%%%%%%%%conditions
figure_shape = 'ellipse';
dispersion = 'no';
N_intervals_x = 3;
N_intervals_y = 3;
N_b = 6;
n_points = 150;
N_basis_x = N_b*ones(N_intervals_x,1);
N_basis_y = N_b*ones(N_intervals_y,1);
lambda = linspace(1000,1700,100);
theta = 0*pi/180;
phi = 0*pi/180;
R1 = 242;
R2 = 242;
P1 = 666;
P2 = 666;
Q2 = R2/sqrt(2);
Q1 = (R1/R2)*sqrt(R2^2-Q2^2);
ellipse_parameters = [R1 R2 P1 P2 Q1 Q2];
b_x1 = zeros(N_intervals_x+1);
b_x2 = zeros(N_intervals_y+1);
x1_t_plus = P1/2+Q1;
x1_t_minus = P1/2-Q1;
x2_t_plus = P2/2+Q2;
x2_t_minus = P2/2-Q2;
b_x1 = [0 x1_t_minus x1_t_plus P1];
b_x2 = [0 x2_t_minus x2_t_plus P2];
[Nx1, NNxx1] = size(b_x1);
[Nx2, NNxx2] = size(b_x2);
periodx = b_x1(NNxx1)-b_x1(1);
periody = b_x2(NNxx2)-b_x2(1);
%delta is the angle between E and the incidence plane
%delta = 0 TM, delta = pi/2 TE
delta = 0;
nSi = 3.5;
nSiO2 = 1.45;
epsSi = nSi^2;
epsSiO2 = nSiO2^2;
epsAir = 1.0;
refIndices = [1.0 nSi];
L=4;
epsilon = zeros(L,2);
epsilon(4,:) = epsAir*ones(1,2);
epsilon(3,:) = epsAir*ones(1,2);
epsilon(2,:) = epsSiO2*ones(1,2);
epsilon(1,:) = epsSi*ones(1,2);
epsilon(3,2) = epsSi; %ellipses
h = zeros(L,1);
h(4) = 0.0;
h(3) = 220;
h(2) = 2000;
h(1) = 0.0;
alpha_ref = -sin(pi/6)/periodx;
beta_ref = -sin(pi/6)/periody;
tau_x = exp(1j*alpha_ref*periodx);
tau_y = exp(1j*beta_ref*periody);
%La is lambda in Gegenbauer polynomials; La>-1/2
La = 0.5;
N_FMM = 1;
%calculate reflection and transmission
[Rsum_ellipses,Tsum_ellipses] = ...
PMM_main_function(figure_shape, dispersion, lambda, theta, phi, delta,...
h, L, N_FMM, epsilon, refIndices, La, tau_x, tau_y, alpha_ref, beta_ref,...
b_x1, b_x2, N_basis_x, N_basis_y, N_intervals_x, N_intervals_y,ellipse_parameters,...
n_points, eta, f1, verbose);
L=3;
epsilon = zeros(L,2);
epsilon(3,:) = epsAir*ones(1,2);
epsilon(2,:) = epsSiO2*ones(1,2);
epsilon(1,:) = epsSi*ones(1,2);
h = zeros(L,1);
h(3) = 0.0;
h(2) = 2000;
h(1) = 0.0;
[Rsum_etched,Tsum_etched] = ...
PMM_main_function(figure_shape, dispersion, lambda, theta, phi, delta,...
h, L, N_FMM, epsilon, refIndices, La, tau_x, tau_y, alpha_ref, beta_ref,...
b_x1, b_x2, N_basis_x, N_basis_y, N_intervals_x, N_intervals_y,ellipse_parameters,...
n_points, eta, f1, verbose);
T_ref_T0 = Tsum_ellipses./Tsum_etched;
figure(1)
plot(lambda, transpose(T_ref_T0), 'r', 'Linewidth', 2);
%plot(lambda, transpose(Rsum), 'r', 'Linewidth', 2);
hold off