% Data plotting script for MATLAB
% (c) 2010 Michael Kelzenberg

% Prompt user for file location
global mksearchpath
if (~isempty(mksearchpath))
   [FileName,PathName] = uigetfile({'*.asc;*.dat;*.txt;*.csv' '*'},...
     'Select norm and data file - norm file must include ''norm'' in name.',...
     'MultiSelect','on',mksearchpath);
else
   [FileName,PathName] = uigetfile({'*.asc;*.dat;*.txt;*.csv' '*'},...
     'Select norm and data file - norm file must include ''norm'' in name.',...
     'MultiSelect','on');
end

% Figure out which file is which, and determine sweep settings
if (~isempty(FileName))
  if( iscell(FileName))
    
    normfilename = [PathName FileName{1}];
    [pathstr, name, ext, versn] = fileparts(normfilename);
    one_data = dlmread( normfilename );

    samplefilename = [PathName FileName{2}];
    [pathstr, name, ext, versn] = fileparts(samplefilename);
    two_data = dlmread( samplefilename );

    if isempty( strfind( [name ext], 'norm' ) )
        norm_data = one_data; sample_data = two_data;
    else
        sample_data = one_data; norm_data = two_data;
    end   

    probablename = name;
    nameend = strfind(name,'_');
    if (~isempty(nameend))
        probablename = name(1:nameend(1)-1);
    end

    probabletype = 'T';
    if (~isempty(strfind(name,'refl')) || ~isempty(strfind(name,'br')))
        probabletype = 'R';
    end

    sampleWL = sample_data(:,2);
    wl_range = [ min(sampleWL) (sampleWL(2)-sampleWL(1)) max(sampleWL)]; 
                                                           % (min step max)
    numpersweep = length(wl_range(1):wl_range(2):wl_range(3));
    numsweeps = length(sampleWL)/numpersweep;
    sampleTheta = sample_data(:,1);
    
    normWL = norm_data(:,2);
    WL_max = min(max(normWL),wl_range(3)); % truncate data range to that
    WL_min = max(min(normWL),wl_range(1)); % available in normalization file
    a1name = '\lambda (nm)';
    if (WL_max < 400) a1name = '\theta_y (\circ)'; end
    if (WL_max < 10 ) a1name = 'x (mm)'; end
        
    if (numsweeps > 1)
        a2name = '\theta_x (\circ)';
        a2range = [num2str(min(sampleTheta)) ' ' num2str(max(sampleTheta))];
        dTheta = sampleTheta(numpersweep+1)-sampleTheta(1);
        if (max(sampleTheta) < 10) a2name = 'y (mm)'; end
    else
        a2name = 'N/A'; a2range = 'N/A'; dTheta = 0;
    end
    
    global mkrunavg mkdeglitch;
    if (isempty(mkrunavg) || isempty(mkdeglitch))
        mkdeglitch = '1054 1074'; mkrunavg = '0';
    end

    %Now prompt user for user input
    defAns = {probablename '' probabletype '' '\lambda (nm)' ...
       [num2str(WL_min) ' ' num2str(WL_max)] a2name a2range ...
       mkrunavg mkdeglitch };            
    sampleinfo = inputdlg({ ...
      'Variable name for MALTAB workspace (cannot be empty)'...
      'Title for graph (if different from above)' ...
      'Normalization type (R, T, N, Rabs, Tabs, or user-defined code)'...
      'Append to previous graph? (write anything here to do so)' ...
      'Axis 1 title' 'Axis 1 range' 'Axis 2 title' 'Axis 2 range' ...
      'Axis 1 running average (num. data points or ''0'' for no averaging)'...
      'Axis 1 deglitch range (min max) [leave empty for none]'}, ...
      'Settings', 1, defAns);        
  end
  mksearchpath = samplefilename;  % remember folder for next time.
end

%Parse user input
if (isempty(sampleinfo{2})) titlename = sampleinfo{1};    %graph title name
else titlename = sampleinfo{2}; end  
shouldhold = ~isempty(sampleinfo{4});           %whether to open new window
userwlrange = str2num(sampleinfo{6});
userthetarange = str2num(sampleinfo{8});
WL_max = min(userwlrange(2),WL_max);            %truncate wavelength range?
WL_min = max(userwlrange(1),WL_min);
if (numsweeps > 1)                                   %truncate theta range?
   Theta_min = max(min(sampleTheta),userthetarange(1)); 
   Theta_max = min(max(sampleTheta),userthetarange(2));
end
deglitchrange = 0;   %whether or not to deglitch the noisy peak at ~1064 nm
if ~isempty(sampleinfo{10}) 
    deglitchrange = str2num(sampleinfo{10}); 
    if (deglitchrange) 
        deglitchrange = min(WL_max, deglitchrange);
        deglitchrange = max(WL_min, deglitchrange);
    end
end
runavg = floor(abs(str2num(sampleinfo{9})));              % running average
exptype = sampleinfo{3};      %What type of normalization (T, R, or other?)
samplename = sampleinfo{1};               %base name for workspace variable

%Save deglitch/avg settings for next time program is run
mkrunavg = num2str(runavg); mkdeglitch = mat2str(deglitchrange); 

%Initialize arrays for normalized data
WL = WL_min:wl_range(2):WL_max;
nWL = length(WL);

if ( numsweeps > 1)
    Theta = Theta_min:dTheta:Theta_max;
    nT = length(Theta);
else
    Theta = 0; nT = 1; Theta_min = 0; Theta_max = 0;
end

normalizedData = zeros(nT,nWL);
sampleReading = sample_data(:,3)./sample_data(:,4);

%Re-sample the normalization data
normReadingOrig = norm_data(:,3)./norm_data(:,4);
normWLOrig = normWL;
normReading = interp1(normWLOrig, normReadingOrig, WL);

%Iterate through the file data to construct the 2D array of normalized data
ThetaIndex = 1;
WLIndex = 1;
recognizedNormalization = 1;
for n=1:length(sampleReading)
    thisWL = sampleWL(n);
    thisTheta = sampleTheta(n);
    
    if (thisTheta < Theta_min || thisTheta > Theta_max) continue
    elseif (thisWL < WL_min || thisWL > WL_max) continue
    end
    
    thisRelative = sampleReading(n)/normReading(WLIndex);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  Normalization - this is where we determine what the "reference"
    %  value is.  We can add additional experiment types here, e.g. for a 
    %  reference photodiode (which must be interpolated by wavelength):
    %
    %  if (strcmp(exptype,'EQE'))
    %     refValue = interp1(diodeWL,diodeEQE,thisWL);  
    %  end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (strcmp(exptype,'T'))
        refValue = 0.92; %approx transmission through quartz slide
    elseif (strcmp(exptype,'R'))
        refValue = 0.98; %approx reflectivity of BaSO4 coating on refl. stage
    elseif (strcmp(exptype,'N'))
        refValue = 1; % (Data is divided by max value later)
    elseif (strcmp(exptype,'Tabs') || strcmp(exptype,'Rabs'))
        refValue = 1;  %Data is normalized to unity
    else
        refValue = 1;  %default case
        recognizedNormalization = 0; %This will issue a warning after the loop
    end  
    
    normalizedData(ThetaIndex, WLIndex) = thisRelative * refValue;
    
    WLIndex = WLIndex + 1;
    if (WLIndex > nWL)
        WLIndex = 1;
        ThetaIndex = ThetaIndex + 1;
        if (ThetaIndex > nT)
            break;
        end
    end
end

if (~recognizedNormalization)
    warndlg(['The experiment type ''' exptype ''' does not have an '...
        'associated normalization spectrum in this program.  Data have '...
        'been normalized to unity.']);
end

%apply deglitching
if(deglitchrange)
  numdeglitches = size(deglitchrange); numdeglitches = numdeglitches(1);
  for nd = 1:numdeglitches
    temp_ND = normalizedData;
    deglitchFromIndex=max(floor((deglitchrange(nd,1)-WL_min)/wl_range(2)+1),1);
    deglitchToIndex=min(ceil((deglitchrange(nd,2)-WL_min)/wl_range(2)+1),nWL);
    for nTheta = 1:nT
      for nWavelen = deglitchFromIndex:deglitchToIndex;
        temp_ND(nTheta, nWavelen) = ...
          interp1([deglitchFromIndex deglitchToIndex], ...
            [normalizedData(nTheta,deglitchFromIndex) ...
             normalizedData(nTheta,deglitchToIndex)], nWavelen);                
      end
    end
  normalizedData = temp_ND;
  end
end

%apply running average
if (runavg)
  temp_ND = normalizedData;
  for nTheta = 1:nT
    for nWavelen = 1:length(WL)
      temp_ND(nTheta, nWavelen) = ...
        mean(normalizedData(nTheta, ...
          (max(nWavelen-runavg,1)):(min(nWavelen+runavg,nWL))));
    end
  end
  normalizedData = temp_ND;
end

%Normalize data to unity for experiment type 'N':
if (strcmp(exptype,'N'))
    normalizedData = normalizedData./max(max(normalizedData));
end
    
%Prepare figure for plot
if(shouldhold)
    gcf; hold on; legend(gca,'off')
    prettycolors = {'b' 'r' 'k' 'g' 'c' 'y' 'm'}; 
    linespec = prettycolors{mod(length(get(gca,'Children')),7)+1};
else
    figure; linespec = 'b'; %Note: linespec only affects 2D (line) plots.
end

%Determine human-readable experiment type for plots
if (exptype(1) == 'T')
    measuredQuantity = 'Transmission';
elseif (exptype(1) == 'R')
    measuredQuantity = 'Reflection';
elseif (strcmp(exptype,'EQE'))
    measuredQuantity = 'External Quantum Efficiency';
elseif (strcmp(exptype,'N'))
    measuredQuantity = '(Normalized)';
else
    measuredQuantity = [exptype ' (rel. to reference reading'];
end

%Make 3-D plot for 2-dimensional data sets
if ( nT > 1 )
    hmesh = mesh(Theta, WL, normalizedData','DisplayName',titlename);
    xlabel(sampleinfo{7},'FontName','Helvetica','FontSize',18)
    ylabel(sampleinfo{5},'FontName','Helvetica','FontSize',18)
    zlabel(measuredQuantity,'FontName','Helvetica','FontSize',18)
    view(65,30)
    set(hmesh,'FaceColor','interp','EdgeColor', [0 0 0],'MeshStyle','column')
    set(gca,'FontSize',14,'FontName','Helvetica','LineWidth',2)
    set(gca,'YLim',[WL_min WL_max], 'XLim',[Theta_min Theta_max])
    title(titlename,'FontName','Helvetica','FontSize',18);
    if (exptype(1) == 'R') 
        set(gca,'XDir','reverse'); %This usually looks better
    end 

%2D line plot for 1-dimensional data set    
else 
    plot(WL,normalizedData',linespec,'DisplayName',titlename);
    xlabel(a1name,'FontName','Helvetica','FontSize',18);
    title(titlename,'FontName','Helvetica','FontSize',18);
    ylabel(measuredQuantity,'FontName','Helvetica','FontSize',18);
    set(gca,'LineWidth',2,'FontName','Helvetica','FontSize',14);
end

hold on;
if (shouldhold) legend('show'); title(''); end

% Save data to workspace:
assignin('base',[samplename '_' exptype],normalizedData');
assignin('base',[samplename '_WL'],WL);
assignin('base',[samplename '_Theta'],Theta);
