% MATLAB script for flattening 3D optical generation profiles for
% 2D-cylindrical device physics simulations 
% (c) 2010 Michael Kelzenberg
% California Institute of Technology

% This script anticipates the output of my FDTD processing script
% (lumerical_getOptGen) to be loaded as variables in the MATLAB workspace.

% Lumerical automatically "unfolds" the data for the simulation, so that
% data is available throughout the entire unit cell (even though only 1/4
% of the cell was simulated.  Thus, the center of the wire is located at
% the midpoint of the x and y indices.  Due to symmetry, we need only
% process 1/4 of the unit cell.
startx = (1+length(Pabs_x))/2;
% Note: This script assumes a square unit cell with identical x and y grids.

% Buffer for output data:
OptGenPlane = zeros( startx, length(Pabs_z))+1e17;
PlanePixelCt = zeros( startx, length(Pabs_z));

for ny=1:startx;
    
    for nz = 1:length(Pabs_z)
        
        for nx = 1:startx
            
            % "pixel radius" for mapping:
            myrad = round(((nx-1)^2 + (ny-1)^2)^0.5)+1;
            
            if (myrad > startx)
                break;
            end
            
            genHere = Ngen_pavg(startx+nx-1,startx+ny-1,nz);
            OptGenPlane(myrad, nz) = OptGenPlane(myrad,nz) + genHere;
            PlanePixelCt(myrad,nz) = PlanePixelCt(myrad,nz) + 1;
        end
    end
end

OptGen = OptGenPlane ./ PlanePixelCt;
OptGen_x = Pabs_x(startx:end);
OptGen_z = Pabs_z;

% This "touches up" the start and end grid locations for the optical generation 
% data set.  Sometimes the endpoints are slightly off, which can be problematic 
% for the mesh conversion step.
% Note that the geometry parameters are hard-coded here, these must be
% changed if the simulation structure is modified.
OptGen_x(end) = 8e-7;
OptGen_z(end) = 3.75e-5;
OptGen_x(1) = 0;
OptGen_z(1) = -3.75e-5;

%Plot the profile:
imagesc(OptGen_x*1e6,flipud(OptGen_z*1e6),log(OptGen')./log(10))
hold on
imagesc(fliplr(-OptGen_x*1e6),flipud(OptGen_z*1e6),log(OptGen')./log(10))
set(gca,'XLim',[-max(OptGen_x) max(OptGen_x)]);
set(gca,'CLim',[18 21]);
axis image
set(gcf,'Position',[1 1 100 1000])
set(gca,'YTick',[])
set(gca,'XTick',[])
set(gca,'Position',[0.01 0.01 0.98 0.98])
set(gcf,'PaperPosition',[0.25 0.25 8 10.50])



    