# MATLAB optical generation grid conversion script for lumsen project
# (c) 2010 Michael Kelzenberg
# California Institute of Technology

#setdep @previous@
#setdep @node|sde@
#setdep @node|lumcad@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB script for importing arbitrary profiles onto DF-ISE simulation
% grids for use with Sentaurus TCAD.
% 
% Michael Kelzenberg, 2010
%
% This script reads a DF-ISE .dat file to determine the spatial coordinates of 
% of each grid point.  It then generates a new .dat file containing an
% OpticalGeneration profile that was calculated by FDTD.  The key input/output 
% settings and the mapping function are indicated by comments in the code.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%DAT file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This should be a valid DF-ISE .dat file (i.e. generated by mesh or
%noffset3d.  The meshing program must be scripted to store the x- and y-
%position of each vertex of the grid as "PMIUserField0" and
%"PMIUserField1", respectively.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

datFile = 'n@node|sde@_msh.dat';
grdFile = 'n@node|sde@_msh.grd';

%FDTD MAT file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This should be the Matlab MAT file generated by the Lumerical CAD script
%including:
% Pabs_x,Pabs_y		X and Y specification of grid (m)
% freq			Freq. of simulation (Hz)
% Pabs *		Matrix of power absorption (W/m3)
% Ngen *		Matrix of optical generation rate (per cm3 per s)
% IntgPwr *		Total power absorption (W/m)
% Current *		Total photocurrent (A per um device length)
% Absfrac *		Fraction of absorbed light, i.e. Absorption Quantum Efficiency
%   *these variables followed by '_pavg' corespond to partial spectral averaging
%
% Note: presently, only Pabs_x, Pabs_y, and Ngen_pavg are used by this script.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FDTDFile = 'n@node|lumcad1@_lummat.mat';

%Regions to process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%These are the regions to perform the optical generation mesh conversion.
%This must be a cell array of region names, including double-quotes (")
%around each region name!!!
%  Example syntax: regionsToProcess = {'"Base_region"', "Emitter_region"' };
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
regionsToProcess = {'"InnerContact_region"' '"Emitter_region"' '"Base_region"'};

%The output dat and grd files are used for monochromatic-illumination device 
% physics simulations (the next step in this experiment).
outputFile = ['n@node@_optgen.dat'];
outputGrid = ['n@node@_optgen.grd'];

%The export data file (.mat) is saved so that a MATLAB script can sum together 
% multiple single-wavelength simulations to approximate solar illumination.
% We chose a file name that is unique to the wavelength, polarization, and 
% the device physics simulation grid:
exportFile = 'n@node|sde@_OptGen@Wavelen@@Polarization@.mat';

%Number of data values to output per line in output DAT file
numperline = 10;

try

disp('');
disp('-----------------------------------------------------------------');
disp(['OptGenConverter Version 2']);
disp(['(c) 2010 Michael Kelzenberg']);
disp(['California Institute of Technology']);
disp('-----------------------------------------------------------------');
disp(' ');

% A wavelength value of '15' is the signal to assemble 1-sun AM1.5G illumination
if (@Wavelen@ == 15)
	nodenum = @node@;  nodenum_sde = @node|sde@;
	disp(['Invoking am15proc.m script to generate combined-wavlength OptGen' ...
        ' profile...']);
	disp(' ');
	disp('This will fail if the prerequisite single-wavelength FDTD profiles');
	disp('  have not been generated for this device physics simulation grid,');
    disp('  (n@node|sde@)');
	disp(' ');
	am15proc;
	exit(0);
end

disp(['Opening DAT file ' datFile ]);

grd = fopen(datFile);
if (grd < 1)
    error(['Error opening file ' datFile ' for reading.']);
    %exit
end

if ( ~isequal( fgetl(grd), 'DF-ISE text'))
    disp('Error with grid file format: It might not be a DF_ISE text file.');
    disp('Please double-check input file.  The first line should read:');
    disp('     DF-ISE text');
    error('File parse error');
    
end
fln = 1;

verts = [];
regions = {};

nl = fgetl(grd); fln=fln+1;
while( isempty( regexp(nl, 'nb_vertices *= *[0-9]+') ) && ~feof(grd) )
    nl = fgetl(grd); fln=fln+1;
end
tmp=regexp(nl, '[0-9]+','match');
numverts = str2num(tmp{1});
disp(['  File reports ' num2str(numverts) ' vertices']);

nl = fgetl(grd); fln=fln+1;
while( isempty( regexp(nl, 'nb_edges *= *[0-9]+') ) && ~feof(grd) )
    nl = fgetl(grd); fln=fln+1;
end
tmp=regexp(nl, '[0-9]+','match');
numedges = str2num(tmp{1});
disp(['  File reports ' num2str(numedges) ' edges']);

nl = fgetl(grd); fln=fln+1;
while( isempty( regexp(nl, 'nb_elements *= *[0-9]+') ) && ~feof(grd) )
    nl = fgetl(grd); fln=fln+1;
end
tmp=regexp(nl, '[0-9]+','match');
numelems = str2num(tmp{1});
disp(['  File reports ' num2str(numelems) ' elements']);

nl = fgetl(grd); fln=fln+1;
while( isempty( regexp(nl, 'nb_regions *= *[0-9]+') ) && ~feof(grd) )
    nl = fgetl(grd); fln=fln+1;
end
tmp = regexp(nl, '[0-9]+','match');
numregions = str2num(tmp{1});
disp(['  File reports ' num2str(numregions) ' regions']);

%Advance to data section of file...
nl = fgetl(grd); fln=fln+1;
while( isempty( regexp(nl, 'Data.*\{', 'once') ) && ~feof(grd) )
    nl = fgetl(grd); fln=fln+1;
end

if ( feof(grd) )
    disp('Unexpected end-of-file, no data processed.');
    disp(['Line: ' num2str(fln)]);
    error('File parse error.');
end


regionArray = [];
disp(' ');
disp('Reading data points...');
%Main reading loop.  Look for PMIUserField 0 or 1 data sets...
while ~feof(grd)
    
    nl = fgetl(grd); fln=fln+1;
    while ( isempty( regexpi(nl,'\s*function\s*=\s*PMIUserField[01]','once'))...
            && ~feof(grd) )
        nl = fgetl(grd); fln=fln+1;
    end
    if (feof(grd))
        break
    end
    
    tmp = regexp(nl, '[01]','match');
    axisNumber = str2num(tmp{1});
    
    nl = fgetl(grd); fln=fln+1;
    while ( isempty( regexpi(nl, '\s*validity\s*=\s*\[\s*".*"\s*\]', 'once'))...
            && ~feof(grd) )
        nl = fgetl(grd); fln=fln+1;
    end
    
    if (feof(grd))
        error(['File Parse Error near line ' num2str(fln)]);
        break
    end
    tmp = regexp(nl, '".*"','match');
    regionName = tmp{1};
    
    nl = fgetl(grd); fln=fln+1;
    while ( isempty( regexpi(nl, '\s*Values\s*\(\s*[0-9]+\s*\)', 'once') )...
            && ~feof(grd) )
        nl = fgetl(grd); fln=fln+1;
    end
    
    if (feof(grd))
        disp(['File Parse Error near line ' num2str(fln)]);
        break
    end
    tmp = regexp(nl, '[0-9]+','match');
    numElems = str2num(tmp{1});
    
    dataPoints = [];
    while (1)
        nl = fgetl(grd); fln = fln+1;
        if(isempty(regexp(nl,'[0-9]+') ) )
            break
        else
            thisline = regexp(nl,'[\.\-\e\E\+0-9][\s\.\-\e\E\+0-9]*','match');
            thisline = thisline{1};
            dataPoints = [dataPoints str2num(thisline)];
        end
        if ( ~isempty(regexp(nl,'}','once') ))
            break
        end
    end
    
    disp(['  Region ' regionName ' read ' num2str(length(dataPoints)) '/' ...
        num2str(numElems) ' elements for axis ' num2str(axisNumber) ]);
    
    %Error if data points disagree with number stated in header
    if ( numElems ~= length(dataPoints) )
        disp(['Error: number of data points does not match file header']);
        disp(['Parse error near line ' num2str(fln)]);
        error(['File structure error in region ' regionName]);
    end
    
    existingRegion = 0;
    for n=1:length(regionArray)
        canRegion = regionArray{n};
        if (isequal(regionName,canRegion.name))
            existingRegion = n;
        end
    end
    
    if (existingRegion)
        if (axisNumber == 0)
            regionArray{existingRegion}.xdata = dataPoints;
        else
            regionArray{existingRegion}.ydata = dataPoints;
        end
        
        if ~isequal( length(regionArray{existingRegion}.xdata), ...
                length(regionArray{existingRegion}.ydata) )
            disp(['Error: number of x data points does not match number of ' ...
                     'y data points']);
            error(['File structure error in region ' regionName ]);
        end
        
    else
        newRegion.name = regionName;
        if (axisNumber == 0)
            newRegion.xdata = dataPoints;
            newRegion.ydata = [];
        else
            newRegion.ydata = dataPoints;
            newRegion.xdata = [];
        end
	newRegion.gdata = zeros(size(dataPoints));
        regionArray{end+1} = newRegion;
    end
    
end

for n=1:length(regionArray)
    if ~isequal( length(regionArray{n}.xdata), length(regionArray{n}.ydata) )
            disp(['Error: number of x data points does not match number of '...
                'y data points']);
            error(['File structure error in region ' regionArray{n}.name ]);
    end
end

disp(' ');
disp('Completed reading DAT file');
disp(['  Read ' num2str(length(regionArray)) ' region(s)']);
disp(' ');
fclose(grd);

%Now ensure that data was successfully read for all requested regions
regionsToProcess = unique(regionsToProcess);
for n=1:length(regionsToProcess)
    reqName = regionsToProcess{n};
    hasRegion = 0;
    for m=1:length(regionArray)
        if isequal( reqName, regionArray{m}.name )
             hasRegion=1;
        break;
        end
    end
    if ~hasRegion
        disp(['Error: Vertex position information for requested region ' ...
            reqName ' not contained within this grid.']);
        error(['Unable to process region: ' reqName ]);
    end
end

% We're done parsing grid file -- Now we load FDTD results and define the 
% mapping function.  Note that the spatial translation applied to the x- and
% y-coordinates in the mapping function is specific to the project geometry:
disp(' ');
disp(['Loading MAT file ' FDTDFile ]);
load(FDTDFile);
optGenMatrix = OptGen';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Mapping function:
newoptgen = @(xi, yi) interp2(Pabs_x,Pabs_y, optGenMatrix, xi*1e-6,yi*1e-6 );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Now ready to write the output data file...
disp(['Opening output file ' outputFile ]);

ogo = fopen(outputFile,'w');
if (ogo < 1)
    error(['Error opening file ' outputFile ' for writing.']);
end

fprintf(ogo, 'DF-ISE text\n\n');
fprintf(ogo, ...
   'Info {\n  version     = 1.0\n  type        = dataset\n  dimension   = 2\n');
fprintf(ogo, '  nb_vertices = %d\n  nb_edges    = %d\n  nb_faces    = 0\n',...
    numverts, numedges);
fprintf(ogo, '  nb_elements = %d\n  nb_regions  = %d\n  datasets    = [ ',...
    numelems, numregions);
for n=1:length(regionsToProcess)
    fprintf(ogo,'"OpticalGeneration" ');
end
fprintf(ogo, ']\n  functions   = [ ');
    for n=1:length(regionsToProcess)
        fprintf(ogo,'OpticalGeneration ');
    end
fprintf(ogo, ']\n}\n\nData {\n\n');

for n=1:length(regionsToProcess)
    reqName = regionsToProcess{n};
    hasRegion = 0;
    for m=1:length(regionArray)
        if isequal( reqName, regionArray{m}.name )
            hasRegion=m;
            break;
        end
    end
    if (hasRegion)
        reg = regionArray{hasRegion};
         
        disp( ['Proessing Optical Generation for region ' reg.name '...'] );
        
        fprintf(ogo,['Dataset ("OpticalGeneration") {\n  function  = '...
            'OpticalGeneration\n  type      = scalar\n  dimension = 1\n'...
            'location  = vertex\n  validity  = [ ' reg.name ' ]\n' ] );
        fprintf(ogo, '  Values (%d) {\n', length(reg.xdata) );
        
	gdata = zeros(size(reg.xdata));
        nl = 1;
        for nv=1:length(reg.xdata)
            
            ogi = newoptgen(reg.xdata(nv), reg.ydata(nv));
            fprintf(ogo, ' %22e', ogi);
	    gdata(nv) = ogi;
            nl = nl + 1;
            if (nl > 10)
                fprintf(ogo, '\n');
                nl = 1;
            end
        end
        if (nl > 1)
            fprintf(ogo,'\n');
        end
        fprintf(ogo, '  }\n}\n\n');
        
        disp( ['   ' num2str(length(reg.xdata))  ' processed'] );
        regionArray{hasRegion}.gdata = gdata;
      end    
 end    

fprintf(ogo,'\n\n}');
fclose(ogo);
disp(['Finished writing output file ' outputFile ]);

disp(' ');
disp(['Copying from grid file: ' grdFile]);
copyfile(grdFile,outputGrid);
disp(['To grid file: ' outputGrid]);

disp(' ');
disp(['Exporting generation profile: ' exportFile ]);
save( exportFile ,'regionArray','numverts','numedges','numelems','numregions');

disp(' ');
disp('Processing complete!');

exit(0);

catch ME
    disp(ME);
    exit(1); %This will mark the node as 'failed' in SWB.
end
