Skip to content

Commit d3ae630

Browse files
twhughestwhughes
committed
initial commit
0 parents  commit d3ae630

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

68 files changed

+8762
-0
lines changed

‎.DS_Store

6 KB
Binary file not shown.

‎LICENSE

+21
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
MIT License
2+
3+
Copyright (c) 2017 Tyler Hughes
4+
5+
Permission is hereby granted, free of charge, to any person obtaining a copy
6+
of this software and associated documentation files (the "Software"), to deal
7+
in the Software without restriction, including without limitation the rights
8+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9+
copies of the Software, and to permit persons to whom the Software is
10+
furnished to do so, subject to the following conditions:
11+
12+
The above copyright notice and this permission notice shall be included in all
13+
copies or substantial portions of the Software.
14+
15+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21+
SOFTWARE.

‎README.md

+37
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
# DLA Optimization Software
2+
3+
## Contents
4+
This matlab program is used to generate optimized dielectric laser accelerator structures.
5+
The code produces two separate structures:
6+
- optimized for maximum acceleration gradient
7+
- optimized for maximum acceleration gradient divided by the maximum electric field amplitude (either in material or design region).
8+
9+
## Paper information
10+
This code was used to for work related to this paper:
11+
12+
[Method for Computationally Efficient Design of Dielectric Laser Accelerator Structures](https://www.osapublishing.org/oe/abstract.cfm?uri=oe-25-13-15414 "Method for Computationally Efficient Design of Dielectric Laser Accelerator Structures")
13+
14+
## How to run simulation
15+
Just run the script "optimization.m". The FDFD code and sparse matrix LU-factoring code is located in "dependencies/".
16+
The entire directory needs to be added to the path to work correctly
17+
18+
## Citing our work
19+
20+
If you would like to use this code, please cite the paper:
21+
22+
23+
@article{Hughes:17,
24+
author = {Tyler Hughes and Georgios Veronis and Kent P. Wootton and R. Joel England and Shanhui Fan},
25+
journal = {Opt. Express},
26+
keywords = {Gratings; Optical devices; Subwavelength structures},
27+
number = {13},
28+
pages = {15414--15427},
29+
publisher = {OSA},
30+
title = {Method for computationally efficient design of dielectric laser accelerator structures},
31+
volume = {25},
32+
month = {Jun},
33+
year = {2017},
34+
url = {http://www.opticsexpress.org/abstract.cfm?URI=oe-25-13-15414},
35+
doi = {10.1364/OE.25.015414},
36+
}
37+

‎dependencies/.DS_Store

8 KB
Binary file not shown.

‎dependencies/FDFD/.DS_Store

6 KB
Binary file not shown.

‎dependencies/FDFD/FDFD.m

+78
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
function [fields,extra] = FDFD(ER,MuR,RES,NPML,BC,lambda0,Pol,b,AF)
2+
3+
extra = {};
4+
fields = {};
5+
k0 = 2*pi/lambda0;
6+
kinc = [1,0];
7+
dx = RES(1);
8+
dy = RES(2);
9+
10+
NGRID = size(ER);
11+
NGRID2 = size(MuR);
12+
%assert(NGRID==NGRID2,'mu and epsilon matrices are not the same size')
13+
[Nx,Ny] = size(ER);
14+
%compute PML
15+
[sx,sy] = calcpml2d(NGRID,NPML);
16+
%
17+
ERxx = ER./sx.*sy;
18+
ERyy = ER.*sx./sy;
19+
ERzz = ER.*sx.*sy;
20+
MuRxx = MuR./sx.*sy;
21+
MuRyy = MuR.*sx./sy;
22+
MuRzz = MuR.*sx.*sy;
23+
24+
%compute material matrices
25+
N = Nx*Ny;
26+
ERxx_inv = spdiags(1./ERxx(:),0,N,N);
27+
ERyy_inv = spdiags(1./ERyy(:),0,N,N);
28+
ERzz = spdiags(ERzz(:),0,N,N);
29+
MuRxx_inv = spdiags(1./MuRxx(:),0,N,N);
30+
MuRyy_inv = spdiags(1./MuRyy(:),0,N,N);
31+
MuRzz = spdiags(MuRzz(:),0,N,N);
32+
%compute derivative matrices
33+
[DEX,DEY,DHX,DHY] = yeeder(NGRID,k0*RES,BC,kinc/k0);
34+
%construct system matrix A
35+
36+
if (nargin<9)
37+
if (strcmp(Pol,'Ez'))
38+
A = DHX*MuRyy_inv*DEX + DHY*MuRxx_inv*DEY + ERzz;
39+
elseif (strcmp(Pol,'Hz'))
40+
A = DEX*ERyy_inv*DHX + DEY*ERxx_inv*DHY + MuRzz;
41+
else
42+
exit(1)
43+
end
44+
AF = factorize(A);
45+
end
46+
b = b(:);
47+
f = AF\b;
48+
49+
%f = mldivide(A,b);
50+
51+
SX = spdiags(sx(:),0,Nx*Ny, Nx*Ny);
52+
SY = spdiags(sy(:),0,Nx*Ny, Nx*Ny);
53+
Hz = reshape(f,[Nx,Ny]);
54+
Ex = reshape(1i*ERxx_inv*DHY*f,[Nx,Ny]);
55+
Ey = reshape(-1i*ERyy_inv*DHX*f,[Nx,Ny]);
56+
Ez = Hz;
57+
Hx = Ex;
58+
Hy = Ey;
59+
60+
x = [1i*ERxx_inv*DHY*f; -1i*ERyy_inv*DHX*f];
61+
62+
fields.Ex = Ex;
63+
fields.Ey = Ey;
64+
fields.Ez = Ez;
65+
fields.Hx = Hx;
66+
fields.Hy = Hy;
67+
fields.Hz = Hz;
68+
fields.x = x;
69+
70+
extra.eps = ER;
71+
extra.derivatives = {};
72+
extra.derivatives.DEX = DEX;
73+
extra.derivatives.DEY = DEY;
74+
extra.derivatives.DHX = DHX;
75+
extra.derivatives.DHY = DHY;
76+
extra.AF = AF;
77+
78+
end

‎dependencies/FDFD/FDFD_TFSF.m

+90
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
function [fields,extra] = FDFD_TFSF(ER,MuR,RES,NPML,BC,lambda0,Pol,TFSF,k,AF)
2+
3+
extra = {};
4+
fields = {};
5+
k0 = 2*pi/lambda0;
6+
kinc = k/norm(k);
7+
dx = RES(1);
8+
dy = RES(2);
9+
10+
NGRID = size(ER);
11+
NGRID2 = size(MuR);
12+
%assert(NGRID==NGRID2,'mu and epsilon matrices are not the same size')
13+
[Nx,Ny] = size(ER);
14+
%compute PML
15+
[sx,sy] = calcpml2d(NGRID,NPML);
16+
%
17+
ERxx = ER./sx.*sy;
18+
ERyy = ER.*sx./sy;
19+
ERzz = ER.*sx.*sy;
20+
MuRxx = MuR./sx.*sy;
21+
MuRyy = MuR.*sx./sy;
22+
MuRzz = MuR.*sx.*sy;
23+
24+
%compute material matrices
25+
N = Nx*Ny;
26+
ERxx_inv = spdiags(1./ERxx(:),0,N,N);
27+
ERyy_inv = spdiags(1./ERyy(:),0,N,N);
28+
ERzz = spdiags(ERzz(:),0,N,N);
29+
MuRxx_inv = spdiags(1./MuRxx(:),0,N,N);
30+
MuRyy_inv = spdiags(1./MuRyy(:),0,N,N);
31+
MuRzz = spdiags(MuRzz(:),0,N,N);
32+
%compute derivative matrices
33+
[DEX,DEY,DHX,DHY] = yeeder(NGRID,k0*RES,BC,kinc/k0);
34+
%construct system matrix A
35+
36+
xpoints = (0:Nx-1)*dx;
37+
ypoints = (0:Ny-1)*dy;
38+
39+
Q = diag(sparse(TFSF(:)));
40+
41+
fsrc = transpose(exp(1i*2*pi/lambda0*kinc(1)*xpoints))*exp(1i*2*pi/lambda0*kinc(2)*ypoints);
42+
fsrc = fsrc(:);
43+
44+
if (nargin<10)
45+
if (strcmp(Pol,'Hz')==0)
46+
A = DHX*MuRyy_inv*DEX + DHY*MuRxx_inv*DEY + ERzz;
47+
elseif (strcmp(Pol,'Ez')==0)
48+
A = DEX*ERyy_inv*DHX + DEY*ERxx_inv*DHY + MuRzz;
49+
else
50+
exit(1)
51+
end
52+
AF = factorize(A);
53+
end
54+
55+
% size(fsrc)
56+
% size(Q)
57+
% size(A)
58+
59+
b = (Q*A-A*Q)*fsrc;
60+
61+
f = AF\b;
62+
63+
%f = mldivide(A,b);
64+
65+
Hz = reshape(f,[Nx,Ny]);
66+
Ex = reshape(1i*ERxx_inv*DHY*f,[Nx,Ny]);
67+
Ey = reshape(-1i*ERyy_inv*DHX*f,[Nx,Ny]);
68+
Ez = Hz;
69+
Hx = Ex;
70+
Hy = Ey;
71+
72+
x = [1i*ERxx_inv*DHY*f; -1i*ERyy_inv*DHX*f];
73+
74+
fields.Ex = Ex;
75+
fields.Ey = Ey;
76+
fields.Ez = Ez;
77+
fields.Hx = Hx;
78+
fields.Hy = Hy;
79+
fields.Hz = Hz;
80+
fields.x = x;
81+
82+
extra.eps = ER;
83+
extra.derivatives = {};
84+
extra.derivatives.DEX = DEX;
85+
extra.derivatives.DEY = DEY;
86+
extra.derivatives.DHX = DHX;
87+
extra.derivatives.DHY = DHY;
88+
extra.AF = AF;
89+
90+
end

‎dependencies/FDFD/FDFD_fast.m

+59
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
function [fields,extra] = FDFD(ER,MuR,RES,NPML,BC,lambda0,Pol,b,AF)
2+
3+
extra = {};
4+
fields = {};
5+
k0 = 2*pi/lambda0;
6+
kinc = [1,0];
7+
8+
NGRID = size(ER);
9+
%assert(NGRID==NGRID2,'mu and epsilon matrices are not the same size')
10+
[Nx,Ny] = size(ER);
11+
%compute PML
12+
[sx,sy] = calcpml2d(NGRID,NPML);
13+
%
14+
ERxx = ER./sx.*sy;
15+
ERyy = ER.*sx./sy;
16+
ERzz = ER.*sx.*sy;
17+
MuRxx = MuR./sx.*sy;
18+
MuRyy = MuR.*sx./sy;
19+
MuRzz = MuR.*sx.*sy;
20+
21+
%compute material matrices
22+
N = Nx*Ny;
23+
ERxx_inv = spdiags(1./ERxx(:),0,N,N);
24+
ERyy_inv = spdiags(1./ERyy(:),0,N,N);
25+
ERzz = spdiags(ERzz(:),0,N,N);
26+
MuRxx_inv = spdiags(1./MuRxx(:),0,N,N);
27+
MuRyy_inv = spdiags(1./MuRyy(:),0,N,N);
28+
MuRzz = spdiags(MuRzz(:),0,N,N);
29+
%compute derivative matrices
30+
[DEX,DEY,DHX,DHY] = yeeder(NGRID,k0*RES,BC,kinc/k0);
31+
%construct system matrix A
32+
33+
if (nargin<9)
34+
if (strcmp(Pol,'Hz')==0)
35+
A = DHX*MuRyy_inv*DEX + DHY*MuRxx_inv*DEY + ERzz;
36+
elseif (strcmp(Pol,'Ez')==0)
37+
A = DEX*ERyy_inv*DHX + DEY*ERxx_inv*DHY + MuRzz;
38+
else
39+
exit(1)
40+
end
41+
AF = factorize(A);
42+
end
43+
b = b(:);
44+
f = AF\b;
45+
46+
x = [1i*ERxx_inv*DHY*f; -1i*ERyy_inv*DHX*f];
47+
y = f;
48+
49+
fields.x = x;
50+
fields.y = y;
51+
52+
extra.derivatives = {};
53+
extra.derivatives.DEX = DEX;
54+
extra.derivatives.DEY = DEY;
55+
extra.derivatives.DHX = DHX;
56+
extra.derivatives.DHY = DHY;
57+
extra.AF = AF;
58+
59+
end

‎dependencies/FDFD/calcpml2d.m

+40
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
function [sx,sy] = calcpml2d(NGRID,NPML)
2+
3+
eta0 = 376.73; % free space impedence
4+
a_max = 2.0; % 0 <= a_max <= 5
5+
p = 4.0; % 3 <= p <= 5
6+
sigma_p_max = 0.008; % sigma_prime ~ 1
7+
8+
Nx = NGRID(1); Ny = NGRID(2);
9+
10+
Nxlo = NPML(1); Nxhi = NPML(2); Nylo = NPML(3); Nyhi = NPML(4);
11+
12+
13+
sx = ones(Nx,Ny);
14+
sy = ones(Nx,Ny);
15+
16+
for nx = (1:Nxlo)
17+
xp = nx/Nxlo;
18+
sig_p_x = sigma_p_max*(sin(pi*xp/2)^2);
19+
a = 1 + a_max*(xp)^p;
20+
sx(Nxlo-nx+1,:) = a*(1+1i*eta0*sig_p_x);
21+
end
22+
for nx = (1:Nxhi)
23+
xp = nx/Nxhi;
24+
sig_p_x = sigma_p_max*(sin(pi*xp/2)^2);
25+
a = 1 + a_max*(xp)^p;
26+
sx(Nx-Nxhi+nx,:) = a*(1+1i*eta0*sig_p_x);
27+
end
28+
for ny = (1:Nylo)
29+
yp = ny/Nylo;
30+
sig_p_y = sigma_p_max*(sin(pi*yp/2)^2);
31+
a = 1 + a_max*(yp)^p;
32+
sy(:,Nylo-ny+1) = a*(1+1i*eta0*sig_p_y);
33+
end
34+
for ny = (1:Nyhi)
35+
yp = ny/Nyhi;
36+
sig_p_y = sigma_p_max*(sin(pi*yp/2)^2);
37+
a = 1 + a_max*(yp)^p;
38+
sy(:,Ny-Nyhi+ny) = a*(1+1i*eta0*sig_p_y);
39+
end
40+
end

‎dependencies/FDFD/calcpml_2d.m

+40
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
function [sx,sy] = calcpml2d(NGRID,NPML)
2+
3+
eta0 = 376.73; % free space impedence
4+
a_max = 3.0; % 0 <= a_max <= 5
5+
p = 3.0; % 3 <= p <= 5
6+
sigma_p_max = 1.0; % sigma_prime ~ 1
7+
8+
Nx = NGRID(1); Ny = NGRID(2);
9+
10+
Nxlo = NPML(1); Nxhi = NPML(2); Nylo = NPML(3); Nyhi = NPML(4);
11+
12+
13+
sx = ones(Nx,Ny);
14+
sy = ones(Nx,Ny);
15+
16+
for nx = (1:Nxlo)
17+
xp = nx/Nxlo;
18+
sig_p_x = sigma_p_max*(sin(pi*xp/2)^2);
19+
a = 1 + a_max*(xp)^p;
20+
sx(Nxlo-nx+1,:) = a*(1+1i*eta0*sig_p_x);
21+
end
22+
for nx = (1:Nxhi)
23+
xp = nx/Nxhi;
24+
sig_p_x = sigma_p_max*(sin(pi*xp/2)^2);
25+
a = 1 + a_max*(xp)^p;
26+
sx(Nx-Nxhi+nx,:) = a*(1+1i*eta0*sig_p_x);
27+
end
28+
for ny = (1:Nylo)
29+
yp = ny/Nylo;
30+
sig_p_y = sigma_p_max*(sin(pi*yp/2)^2);
31+
a = 1 + a_max*(yp)^p;
32+
sy(:,Nylo-ny+1) = a*(1+1i*eta0*sig_p_y);
33+
end
34+
for ny = (1:Nyhi)
35+
yp = ny/Nyhi;
36+
sig_p_y = sigma_p_max*(sin(pi*yp/2)^2);
37+
a = 1 + a_max*(yp)^p;
38+
sy(:,Ny-Nyhi+ny) = a*(1+1i*eta0*sig_p_y);
39+
end
40+
end

‎dependencies/FDFD/little_yeeder.p

508 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)