-
Notifications
You must be signed in to change notification settings - Fork 13
/
PhC2D_hex_PWE_f.m
86 lines (60 loc) · 2.81 KB
/
PhC2D_hex_PWE_f.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
function[E,f0]=PhC2D_hex_PWE_f(Xhex,Yhex,Gxhex,Gyhex,k,HHH,nmodes,TE,TM);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c=2.99792458e8;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Interpolation on a grid that have 2^N points %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nx=length(Xhex(1,:));
Ny=length(Yhex(:,1));
NGx=length(Gxhex(1,:));
NGy=length(Gyhex(:,1));
NG=NGx*NGy;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Building Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GXX=Gxhex(:);
GYY=Gyhex(:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if TE==1
GkXX = ( GXX + k(1) )*( GXX + k(1) )'; % Gk(i,j) = (G(i) + k)*(G(j) + k)
GkYY = ( GYY + k(2) )*( GYY + k(2) )';
Gk=GkXX+GkYY;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if TM==1
Gk1 = sqrt( ( GXX + k(1) ) .^2 + ( GYY + k(2) ) .^2 ) ;
Gk2 = sqrt( ( GXX + k(1) )'.^2 + ( GYY + k(2) )'.^2 );
Gk=Gk1*Gk2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
HH=reshape(Gk,[NGy,NGx,NGy,NGx]);
H=HH.*HHH;
H=reshape(H,NG,NG);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Solving Hamiltonian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[psik,k0] = eig(H); %% eigen values are ordered
f0 = sqrt(diag(k0)) ; % actually it is w0
lambda= 2*pi ./ sqrt(diag(k0)) * 1e6 ;
f0=f0(1:nmodes);
psik = psik(:,1:nmodes);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% Reverse Hexagonal Fourier Transform %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:nmodes
PSI = reshape(psik(:,j),[NGy,NGx]);
Ghex=PSI;
for m=1:Nx
for l=1:Ny
wwhex = exp(1i*( (Xhex(l,m)-Xhex(1))*(Nx-1)/Nx *Gxhex + (Yhex(l,m)-Yhex(1))*(Ny-1)/Ny *Gyhex ) );
Fhex(l,m) = (1/(Nx*Ny))*sum(sum(Ghex.*wwhex));
end
end
E(:,:,j)= Fhex /max(Fhex(:));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%