Add files via upload

Functions/BlurGeomPostGrad.m

@@ -0,0 +1,15 @@
+% Apply blurring filters when specificed in OptParm
+function DevicePattern = BlurGeomPostGrad(DevicePattern, iter, OptParm, GridScale)  
+    MaxIterations = OptParm.Optimization.Iterations;
+
+    % Large blur every X iterations
+    if ((mod(iter,OptParm.Optimization.Filter.BlurLargeIter)==0)&&(iter<(MaxIterations - OptParm.Optimization.Filter.BlurLargeIterStop)))
+        FilterLarge = fspecial('disk',0.5*floor(OptParm.Optimization.Filter.BlurRadiusLarge/GridScale));
+        DevicePattern = imfilter(DevicePattern,FilterLarge,'circular');
+
+    % Small blur every Y iterations
+    elseif ((mod(iter,OptParm.Optimization.Filter.BlurSmallIter)==0)&&(iter<(MaxIterations - OptParm.Optimization.Filter.BlurSmallIterStop)))
+        FilterSmall = fspecial('disk',OptParm.Optimization.Filter.BlurRadiusSmall/GridScale);
+        DevicePattern = imfilter(DevicePattern,FilterSmall,'circular');
+    end
+end

Functions/DefineGrid.m

@@ -0,0 +1,27 @@
+% Compute the simulation grid for given geometry
+function [xGrid, yGrid, dr] = DefineGrid(Grid, Period, Wavelength)
+
+    %Number of grid points
+    NGrid = ceil(Grid*Period/Wavelength);
+    Nx = NGrid(1);
+    Ny = NGrid(2);
+
+    %Device period
+    Px = Period(1);
+    Py = Period(2);
+
+    %Compute external grid coordinates
+    xBounds = linspace(0,Px,Nx+1); 
+    yBounds = linspace(0,Py,Ny+1);
+
+    %Compute size of each grid box
+    dx = xBounds(2) - xBounds(1);
+    dy = yBounds(2) - yBounds(1);
+
+    %Compute coordinates of center of each box
+    xGrid = xBounds(2:end)- 0.5*dx;
+    yGrid = yBounds(2:end)- 0.5*dy;
+
+    %Compute average grid size
+    dr = mean([dx dy]);
+end

Functions/DensityFilter2D.m

@@ -0,0 +1,18 @@
+% Distance weighted averaging filter of radius R
+function PatternOut = DensityFilter2D(PatternIn,Radius)
+    % If the radius is less than 1 pixel, no filter can be applied
+    if Radius<1
+        PatternOut = PatternIn;
+    else
+        % Define grid
+        [X1,X2] = meshgrid(-ceil(Radius):ceil(Radius),-ceil(Radius):ceil(Radius));
+
+        % Compute weights
+        Weights = Radius - (X1.^2+X2.^2).^(1/2);
+        Weights(Weights<0) = 0;
+        B = sum(Weights(:));
+
+        %Apply filter
+        PatternOut = imfilter(PatternIn, Weights/B, 'circular');
+    end
+end

Functions/EnforceSymmetry.m

@@ -0,0 +1,14 @@
+% Enforces required symmetries in the pattern
+% by folding over midline and averaging
+function PatternOut = EnforceSymmetry(PatternIn, SymX, SymY)
+    PatternOut = PatternIn;
+    % Enforce X symmetry
+    if SymX
+        PatternOut = 0.5*(PatternIn+flipud(PatternIn));
+    end
+
+    % Enforce Y symmetry
+    if SymY
+        PatternOut = 0.5*(PatternIn+fliplr(PatternIn));
+    end
+end

Functions/FilteredGrad2D.m

@@ -0,0 +1,25 @@
+% Computes the base pattern gradient from the density and threshold
+% filtered gradient
+function GradientOut = FilteredGrad2D(GradientIn,PatternIn,Bin,Midpoint,Radius)
+
+% Compute the derivative of the threshold filter
+% Threshold filter is a piecewise function centered around the Midpoint
+if Bin~=0
+    PatternLow = Bin*(exp(-Bin*(1-PatternIn/Midpoint)))+exp(-Bin);
+    PatternHigh = exp(-Bin) + Bin*exp(-Bin*(PatternIn-Midpoint)/(1-Midpoint));
+else
+    PatternLow = ones(size(PatternIn))/(2*Midpoint);
+    PatternHigh = ones(size(PatternIn))/(2*(1-Midpoint));
+end
+
+% Combine the two pieces of the threshold derivative
+PatternLow(PatternIn>Midpoint) = 0;
+PatternHigh(PatternIn<=Midpoint) = 0;
+PatternDeriv = PatternLow + PatternHigh;
+
+% Apply chain rule
+Gradient = PatternDeriv.*GradientIn;
+
+% Chain rule of the DensityFilter is another DensityFilter
+GradientOut = DensityFilter2D(Gradient,Radius);
+end

Functions/FineGrid.m

@@ -0,0 +1,52 @@
+% Transfers the pattern onto a finer/coarser grid specifid by the scaling
+% factor Scale
+function [PatternOut,XOut,YOut] = FineGrid(PatternIn,Period,Scale,Binarize,SmoothGeom)
+
+Pattern = PatternIn;
+[Nx, Ny] = size(PatternIn); % Size of input pattern
+
+% Compute size of output pattern
+NScaled = round(Scale.*[Nx Ny]);
+NxScaled = NScaled(1);
+NyScaled = NScaled(2);
+
+% Blur pattern if option is specfied
+if SmoothGeom == 1
+    Pattern = imgaussfilt(Pattern,Ny/50,'Padding','circular');
+end
+
+% Add extra pixels in order to allow for periodic interpolation
+Pattern = [Pattern(:,end),Pattern];
+Pattern = [Pattern(end,:);Pattern];
+
+% Input grid
+[X,Y] = meshgrid(0:Nx,0:Ny);
+
+% Scaled grid
+[XScaled,YScaled] = meshgrid(linspace(Nx/NxScaled,Nx,NxScaled),...
+    linspace(Ny/NyScaled,Ny,NyScaled));
+
+% Output coordinates
+[XOut,YOut] = meshgrid(linspace(0,Period(1),NxScaled),...
+    linspace(0,Period(2),NyScaled));
+
+if SmoothGeom == 1
+    % Use interpolation if smooth geometry is desired
+    PatternOut = interp2(X,Y,Pattern',XScaled,YScaled,'cubic')';
+else
+    % Place pixels on a square grid otherwise
+    PatternOut = zeros([NxScaled NyScaled]);
+    for ii = 1:NxScaled
+        for jj = 1:NyScaled
+            PatternOut(ii,jj) = Pattern(ceil(ii*Nx/NxScaled),...
+                ceil(jj*Ny/NyScaled));
+        end
+    end
+end
+
+%Binarize if desired
+if Binarize==1
+    PatternOut(PatternOut<0.5) = 0;
+    PatternOut(PatternOut>=0.5) = 1;
+end
+end

Functions/FractureGeom.m

@@ -0,0 +1,78 @@
+% Divides the given geometry into rectangles to be used in Reticolo
+function GeometryOut = FractureGeom(PatternIn,nLow,nHigh,XGrid,YGrid)
+
+% Acceptable refractive index tolerance in fracturing
+Tolerance = .01;
+
+% Extract grid parameters
+dX = XGrid(2)-XGrid(1);
+dY = YGrid(2)-YGrid(1);
+[Nx, Ny] = size(PatternIn);
+
+Geometry = {nLow}; %Define background index
+
+% Define the areas with max refractive index (within tolerance)
+% For complex values, tests the real part
+PatternHigh = (PatternIn >= (nHigh - Tolerance));
+
+% Define areas with higher than background refractive index
+PatternElse = (PatternIn > (nLow + Tolerance));
+
+MaxFractures = 2000;
+RectCount = 1;
+
+% Divide given PatternHigh geometry into rectangles
+while (sum(PatternHigh(:)) > 0) && (RectCount <= MaxFractures)
+    [X1,Y1] = find(PatternHigh==1,1);
+    Rect = [X1,X1,Y1,Y1];
+
+    % Find largest silicon rectangle from inital point
+    Fracturing = 1;
+    while Fracturing
+        % Attempt to expand rectangle in x direction
+        if (Rect(2)<Nx)&&(sum(sum(PatternHigh(Rect(2)+1,Rect(3):Rect(4))))==(Rect(4)-Rect(3)+1))
+            Rect(2) = Rect(2) + 1;
+
+        % Attempt to expand in y direction
+        elseif (Rect(4)<Ny)&&(sum(sum(PatternHigh(Rect(1):Rect(2),Rect(4)+1)))==(Rect(2)-Rect(1)+1))
+            Rect(4) = Rect(4) + 1;
+        else
+            Fracturing = 0;
+        end
+    end
+
+    % Remove pixels of fractured areas
+    PatternHigh(Rect(1):Rect(2),Rect(3):Rect(4)) = 0;
+    PatternElse(Rect(1):Rect(2),Rect(3):Rect(4)) = 0;
+
+    %Record size of the rectangle
+    dXGrid = 0.5*(XGrid(Rect(1)) + XGrid(Rect(2))); % Grid point size
+    dYGrid = 0.5*(YGrid(Rect(3)) + YGrid(Rect(4))); 
+    dXLength= dX*(Rect(2)-Rect(1)+1); % Physical Size
+    dYLength = dY*(Rect(4)-Rect(3)+1); 
+
+    if isequal(Geometry{end},[dXGrid,dYGrid,dXLength,dYLength,nHigh,1])
+        disp('Warning: duplicate structure');
+        break;
+    end
+
+    % Add rectangle to list
+    Geometry = [Geometry,{[dXGrid,dYGrid,dXLength,dYLength,nHigh,1]}];
+
+    RectCount = RectCount + 1;
+    if RectCount > (MaxFractures * 0.95)
+        disp('Warning: High fracture count');
+    end
+end
+
+% Fracture non binarized pixels
+for i = 1:Nx % Defining texture for patterned layer. Probably could have vectorized this.
+    for j = 1:Ny
+        if PatternElse(i,j)==1
+            % Add each reamaining pixel and its refractive index independently
+            Geometry = [Geometry,{[XGrid(i),YGrid(j),dX,dY,PatternIn(i,j),1]}];
+        end
+    end
+end
+GeometryOut = Geometry;
+end

Functions/GaussFilter2D.m

@@ -0,0 +1,4 @@
+%Gaussian weighted averaging within radius Sigma
+function PatternOut = GaussFilter2D(PatternIn,Sigma)
+    PatternOut = imgaussfilt(PatternIn, Sigma, 'padding', 'circular');
+end

Functions/GaussGrad2D.m

@@ -0,0 +1,25 @@
+% Computes the base pattern gradient from the Gaussian and threshold
+% filtered gradient
+function GradientOut = GaussGrad2D(GradientIn,PatternIn,Bin,Midpoint,Sigma)
+
+% Compute the derivative of the threshold filter
+% Threshold filter is a piecewise function centered around the Midpoint
+if Bin~=0
+    PatternLow = Bin*(exp(-Bin*(1-PatternIn/Midpoint)))+exp(-Bin);
+    PatternHigh = exp(-Bin) + Bin*exp(-Bin*(PatternIn-Midpoint)/(1-Midpoint));
+else
+    PatternLow = ones(size(PatternIn))/(2*Midpoint);
+    PatternHigh = ones(size(PatternIn))/(2*(1-Midpoint));
+end
+
+% Combine the two pieces of the threshold derivative
+PatternLow(PatternIn>Midpoint) = 0;
+PatternHigh(PatternIn<=Midpoint) = 0;
+PatternDeriv = PatternLow + PatternHigh;
+
+% Apply chain rule
+Gradient = PatternDeriv.*GradientIn;
+
+% Chain rule of Gaussian filter is another Gaussian
+GradientOut = DensityFilter2D(Gradient ,Sigma);
+end

Functions/GenerateBVector.m

@@ -0,0 +1,20 @@
+% Defines the values of B, the binarization rate, at each iteration
+% Plot BVector to see the progression of binarization
+function BVector = GenerateBVector(MaxIterations, BinParm)
+    % Extract parameters
+    BMin = BinParm.Min;
+    BMax = BinParm.Max;
+    BStart = BinParm.IterationStart;
+    BHold = BinParm.IterationHold;
+
+    BMid = BMax/20;
+
+    BVector = zeros(1,MaxIterations);
+
+    Bmult1 = (BMid/BMin)^(1/floor((round(MaxIterations/2)-BStart)/BHold));
+    Bmult2 = (BMax/BMid)^(1/floor((round(MaxIterations/2))/BHold));
+
+    % The binarization speed is a piecewise function
+    BVector((BStart+1):round(MaxIterations/2)) = BMin*Bmult1.^(floor((1:(round(MaxIterations/2)-BStart))/BHold));
+    BVector((round(MaxIterations/2)+1):end) = BMid*Bmult2.^(floor((1:(MaxIterations-round(MaxIterations/2)))/BHold));
+end

Functions/GenerateThreshVectors.m

@@ -0,0 +1,25 @@
+% Interpolate between the start and end edge deviations to create
+% and edge deviation value for each iteration
+function [ThresholdVectors, NRobustness] = GenerateThreshVectors(OptParm)
+    % Extract necessary parameters
+    MaxIterations = OptParm.Optimization.Iterations;
+    Start = OptParm.Optimization.Robustness.StartDeviation;
+    End = OptParm.Optimization.Robustness.EndDeviation;
+    Ramp = OptParm.Optimization.Robustness.Ramp;
+    NRobustness = length(Start);
+
+    % Error if the dimensions don't match
+    if NRobustness ~= length(End)
+       error('Robustness vectors are not the same length!');
+    end
+
+    % Generate edge deviation vector
+    DeviationVectors = zeros(NRobustness, MaxIterations);
+    for ii = 1:NRobustness
+       DeviationVectors(ii, 1:Ramp) = linspace(Start(ii), End(ii), Ramp);
+       DeviationVectors(ii, (Ramp+1):end) = End(ii);
+    end
+
+    % Generate corresponding thresholding values
+    ThresholdVectors = 0.5*(1-erf(DeviationVectors/OptParm.Optimization.Filter.BlurRadius));
+end

Functions/Initialize.m

@@ -0,0 +1,72 @@
+% Initialize default values for necessary optimization parameters
+function [OptParm] = Initialize()
+addpath('reticolo_allege')
+OptParm.Input.Wavelength = 1000; % Wavelength of incident light
+OptParm.Input.Polarization = 'TM'; % Polarization: can be 'TM', 'TE' or 'Both'
+OptParm.Input.Theta = 0; % Incident angle (angle given in air)
+
+OptParm.Geometry.Period = [1500, 500]; % Period of the device
+OptParm.Geometry.Thickness = 350; % Device layer thickness
+OptParm.Geometry.Substrate = 1.45; % Refractive index of substrate
+OptParm.Geometry.Device = 3.45; % Refractive index of device
+OptParm.Geometry.Top = 1; % Refractive index of upper layer
+
+% Symmetry enables faster computation
+OptParm.Geometry.SymmetryY = 1; % Enforces symmetry in the Y direction
+OptParm.Geometry.SymmetryX = 0; % Enforces symmetry in the X direction
+
+OptParm.Simulation.Grid = 90; % Grid points per wavelength 
+OptParm.Simulation.Fourier = [14 14]; % Number of Fourier orders per wavelength in x and y respectively
+OptParm.Simulation.ZGrid = 50; % Number of slices in device layer at which fields are computed
+
+OptParm.Optimization.Target = 1; % Target diffraction order in the x-direction for each polarization i.e. [1 -1]
+
+OptParm.Optimization.Iterations = 350; % Max iterations 
+OptParm.Optimization.Gradient.StepSize = 0.1; % Initial gradient step size
+OptParm.Optimization.Gradient.StepDecline = 0.992; % Multiplying factor that decreases step size each iteration
+
+OptParm.Optimization.Filter.BlurRadiusSmall = 5; % Blur radius for small blurring
+OptParm.Optimization.Filter.BlurRadius = 20; % Blurring radius used for second blurring
+OptParm.Optimization.Filter.BlurRadiusLarge = 60; % Blurring radius for initial blurring
+OptParm.Optimization.Filter.BlurSmallIter = 5; % Frequency of small blurring
+OptParm.Optimization.Filter.BlurSmallIterStop = 25; % Steps at the end of optimization which have no blurring
+OptParm.Optimization.Filter.BlurLargeIter = 60; % Frequency of large blurring
+OptParm.Optimization.Filter.BlurLargeIterStop = 50; % Steps at the end of optimization which have no blurring
+
+% Parameters for defining the binarization rate B. Smaller B is a
+% less-sharp function, a large B will create a step function that cuts off 
+% at the different values of N.
+OptParm.Optimization.Binarize.Min = 1; % Minimum value of B that is used after it is initialized. 
+OptParm.Optimization.Binarize.Max = 50; % Maximum value of B
+OptParm.Optimization.Binarize.IterationStart = 1; % Iteration at which B is allowed to be nonzero.
+OptParm.Optimization.Binarize.IterationHold = 20; % Number of iterations during which B is kept constant before increasing (after it is initialized)
+
+% Start and End should be the same length. Each element represents a
+% separate thread of the optimization at a different flat edge deviation in nm. 
+OptParm.Optimization.Robustness.StartDeviation = [-5 0 5]; % Starting edge deviation values
+OptParm.Optimization.Robustness.EndDeviation = [-5 0 5]; % Ending edge deviation values
+OptParm.Optimization.Robustness.Ramp = 30; % Iterations over which the thresholding parameter changes
+OptParm.Optimization.Robustness.Weights = [.5 1 .5]; % Gradient weight for each robustness value
+
+OptParm.Optimization.Start = []; % Input pattern matrix as starting point. Leave empty to generate random starting point
+
+% Parameters for random starting point. See RandomStart for details. 
+% Ignored if starting point is already provided
+OptParm.Optimization.RandomStart.Pitch = 100; % Grid size on which to create random Gaussians
+OptParm.Optimization.RandomStart.Average = 0.6; % Average refractive index
+OptParm.Optimization.RandomStart.Sigma = 0.5; % Standard deviation of the refractive index
+
+OptParm.Display.ShowText = 1; % Display progress as console output
+OptParm.Display.PlotGeometry = 1; % Display current device geometry as image
+OptParm.Display.PlotEfficiency = 0; % Plot device efficiency vs iteration
+
+% Checkpointing saves the optimization state to a temporary file 
+% and allows restarting from the previous checkpoint if the optimization is
+% interrupted
+OptParm.Checkpoint.Enable = false; % Enable checkpointing behavior
+OptParm.Checkpoint.File = ''; % Checkpoint file name
+OptParm.Checkpoint.Frequency = 10; % Number of iterations between checkpoints
+
+% Disable parfor warning messages
+warning('off','MATLAB:mir_warning_maybe_uninitialized_temporary')
+end