[data,varargout] =sifPlot(varargin)

%SIFPLOT Plot a sif file with options.
%
% Description
% 
%     Plots and saves a sif file taken with Andor Newton camera. The plot
%     data can be written to an exported figure, which gives a quick way of
%     processing sif files. The main benefit of using this function over
%     sifreadnk is the variety of plotting options.
% 
%     sifPlot can also take a data structure as an input (aka the output of
%     sifreadnk).
%
% Input and Output
% 
%     data = sifPlot(fname,options). fname is the name of a sif file to be
%       imported and plotted. options is a list of desired plot options
%       explained below. data is an output structure giving the sif file's
%       characteristics.
% 
%     The data output structure contains a lot of information. Here are
%     some descriptions of the most important fields
% 
%         'imageData': Contains the actual data. The imageData field is an
%         array of dimension [1024, 1, length]. The first dimension, 1024,
%         corresponds to the number of pixels in the x-direction for the
%         camera. The second dimension is 1 when we acquire in full
%         vertical binning (FVB) mode. The 3rd dimension is the length of
%         the kinetic series if more than one frame were acquired.
%
%         'axisWavelength': The x-axis for the spectrum in units of nm.
%
%         'axisEnergy': The x-axis in units of eV. 'axisGHZ': the x-axis in
%         units of GHz.
%
% Options 
%
%     Here are all the possible options and the default values.
% 
%         'xunits'  : [('eV') 'nm' 'GHz' 'pixels'] change x-units of plot
% 
%         'yunits'  : ['normalized' 'counts' {'rate'}]. 'Normalized' scales
%         the data linearly to the interval [0,1]. 'Counts' displays the
%         total number of counts in the spectrum. 'rate' displays the data
%         in units of counts per second.
% 
%         'scale'   : 1. Scales the image data by the number specified.
%         This can help zoom in on a section of a plot.
% 
%         'offset'  : 0. Offsets the plot by the number specified. This is
%         helpful for easily plotting several sif files on top of each
%         other.
% 
%         'crop'    : 'units', cropMin, cropMax. Crops the x axis in the
%         specified units (either 'nm','eV', 'GHz', or 'pixels'). The
%         program finds the xaxis pixel closest to the cropMin and cropMax
%         specified.
% 
%         'frame'   : 1. Give a number to choose which frame of a kinetic
%         series to plot.
% 
%         'genPDF'  : generate a pdf containing the plotted spectrum and
%         some information about it. This is especially useful for quick
%         processing of sif files for later investigation. The file name is
%         the same as the imported .sif file except with the extension .pdf
% 
%         'noPlot'  : do not plot the spectrum. This speeds up execution
%         time massively.
% 
%         'saveData': Saves the generated data as a .mat file with the same
%         name as the imported file name.
% 
%         'filter'  : [filterMin,filterMax] Draws a box around a certain
%         region in the plot.
% 
%         'plotSpec': plotSpec specifies the style of the plotted line.
% 
%         'fontsize': 12  sets font size of plot axes and labels.
%       
%         'interp': numToInterp. Interpolates numToInterp points in the
%         spectrum.
%
% Examples
%
%     The following are examples of function calls. The various options can
%     be conbined. Coding 'data =' gives the imported information about the
%     sif file as a structure. It is not necessary to write 'data ='
% 
%         data = sifPlot('example.sif') : plots a sif file using default
%         settings. data = sifPlot('spectra.sif','xunits', 'eV') : uses eV
%         for x axis. Can also be 'nm', 'GHz' or 'pixels'
%
%         data = sifPlot('spectra.sif','yunits', 'normalize') : normalizes
%         the y-axis. Options are 'normalize, 'counts', 'rate' (counts per
%         sec).
%
%         data = sifPlot('spectra.sif','crop','nm',800,826) : crops the
%         y-axis to 800nm - 826 nm. Can also specify crop dimensions in
%         'nm' or 'pixels'. If crop dimensions are out of bounds, then
%         program will find the closest crop value that is valid.
%
%         data = sifPlot('spectra.sif','cosmicrays',[200 415]) : removes
%         cosmic rays at pixels 200 and 415.
%
%         data = sifPlot('spectra.sif','nopdf') : do not produce a pdf data
%         file, improves program runtime.
%
%         data = sifPlot('spectra.sif','noplot') : does not plot the
%         spectrum. data = sifPlot('spectra.sif','savedata') : saves the
%
%         data as mat file data = sifPlot('spectra.sif','frame',100) :
%         displays the 100th frame of a kinetic series.
%
%         data = sifPlot('spectra.sif','offset',.3') : adds a .3 offset to
%         the plotted data.
%
%         sifPlot('spectra.sif','plotSpec','--k') : draws the spectrum with
%         a dotted black line.
%
% See also:
%   sifreadnk, convertUnits
%
%
% Todd Karin
% 08/22/2012


% Version History
%
% 3/28/2012
% Added a 'nopdf' option
%
% 3/5/2012
% remade input parameters to be easier to use.
%
% 3/1/2012
% Allowed for variable number of input arguments.
%
% 2/25/2012
% Made into a function.
%
% 06/02/2012
% Added some more options: genpdf, genplot, saveData
%
% 06/26/2012
% Added ability to select which frame of a kinetic series to plot.
%
% 06/29/2012
% Cleaned up option specification.


% PROGRAM
function [data,varargout] =sifPlot(varargin) 

% Read options
fname       = varargin{1};
[ret, xswitch]     = readOption('xunits',1,'nm',varargin);
[ret, yunits]   = readOption('yunits',1,'rate',varargin);
[ret, scale]       = readOption('scale',1,1,varargin);
[ret, offset] = readOption('offset',1,0,varargin);
[ret,cropswitch,crop(1),crop(2)] = readOption('crop',3,{'pixels',0,1e7},varargin);
[removeCosmicRays cosmicrays]  = readOption('cosmicRays',1,[],varargin);
genpdf      = readOption('genPDF',0,0,varargin);
noplot      = readOption('noPlot',0,0,varargin);
saveData    = readOption('saveData',0,0,varargin);
[ret,frame]   = readOption('frame',1,1,varargin);
[axesHandleOn,axesHandle]= readOption('axesHandle',1,0,varargin);
[plotRangeOn,plotRange] = readOption('plotrange',1,[1 1],varargin);
[filterOn, filter] = readOption('filter',1,[0,0],varargin);
[fontSizeOn, fontSize] = readOption('fontSize',1,12,varargin);
[ret,plotSpec] = readOption('plotSpec',1,'k',varargin);
[interpOn,interpPoints] = readOption('interp',1,1,varargin);

[ret,xoffset] = readOption('xoffset',1,0,varargin);

% Make a name 'fname2' name without underscores or .sif at the end.
if ischar(fname)
    fname2 = remove_underscores(fname);
else
    fname2 = 'untitled';
end

% Import the spectrum
if ischar(fname)    % from a file
    data = sifreadnk(fname); 
else                % from a data structure
    data = fname;
end

% Interpolate more data points.
if interpOn
    
    axisWavelengthNew = increaseAxisResolution(data.axisWavelength,interpPoints);

    
    data.imageData = interp1( data.axisWavelength, data.imageData, axisWavelengthNew,'spline')';
    
    
    data.axisWavelength = axisWavelengthNew;
    data.axisEnergy = convertUnits(data.axisWavelength,'nm','eV');
    data.axisGHz = convertUnits(data.axisWavelength,'nm','GHz');
end

% Crop the image to desired wavelength range
if strcmpi(cropswitch,'nm')
    data.cropi(1) = findElement(crop(1),data.axisWavelength);
    data.cropi(2) = findElement(crop(2),data.axisWavelength);
elseif strcmpi(cropswitch,'eV')
    data.cropi(1) = findElement(crop(1),data.axisEnergy);
    data.cropi(2) = findElement(crop(2),data.axisEnergy);
elseif strcmpi(cropswitch,'pixels')
    data.cropi(1) = findElement(crop(1),1:length(data.axisWavelength));
    data.cropi(2) = findElement(crop(2),1:length(data.axisWavelength));
elseif strcmpi(cropswitch,'GHz')
    data.cropi(1) = findElement(crop(1),data.axisGHz);
    data.cropi(2) = findElement(crop(2),data.axisGHz);
else
    data.cropi(1) = 1;
    data.cropi(2) = size(data.imageData,1);
end

% Sort crop indices so that program does not mind if max<min
data.cropi = sort(data.cropi);
cax = data.cropi(1):data.cropi(2);

% Cosmic Ray removal (needs work)
if removeCosmicRays
    c=cosmicrays;
    n=4;
    m=3;
    for k=1:length(c)
        lowerave=(c(k)-n):(c(k)-m);
        upperave=(c(k)+m):(c(k)+n);
        data.imageData(c(k),1,frame) = (mean(data.imageData(lowerave))+ ...
            mean(data.imageData(upperave)))/2;
    end
    data.cosmicRays = cosmicrays;
else
    data.cosmicRays = [];
end
  
% x-axis type
if strcmpi(xswitch,'nm')
    data.xaxis = data.axisWavelength(cax);        
elseif strcmpi(xswitch,'eV')
    data.xaxis = data.axisEnergy(cax);
elseif strcmpi(xswitch,'pixels')
    data.xaxis = cax;
elseif strcmpi(xswitch,'GHz')
    data.xaxis = data.axisGHz(cax);
else 
    error('plot axis type not understood')
end

% Offset xaxis
data.xaxis = data.xaxis + xoffset;

% y-axis type: counts, counts/sec or normalized.
if strcmpi(yunits,'counts')
    imagedata = data.imageData(cax,1,frame);
    ylab = 'Counts';
elseif strcmpi(yunits,'rate')
    imagedata = data.imageData(cax,1,frame)/data.exposureTime;
    ylab = 'Counts/sec';
elseif strcmpi(yunits,'normalized')
    imagemin = min(data.imageData(cax,1,frame)); % sets up normalization, for use in cosmic ray removal.
    imagedata1 = data.imageData(cax,1,frame) - imagemin;
    imagemax = max(imagedata1);
    imagedata = imagedata1/imagemax;
    ylab = 'Normalized counts';
else 
    error('plot style not recognized')
end

% Add scaling and an offset
imagedata = scale*imagedata+offset;

% This is the figure that is exported with data attached.
if genpdf

    figure('Visible','off')
    
    if strcmpi(xswitch,'nm')
        subplot(2,1,1), plot(data.axisWavelength(cax), imagedata);
        xlabel('\lambda (nm)','fontsize',fontSize);
    elseif strcmpi(xswitch,'eV')
        subplot(2,1,1), plot(data.axisEnergy(cax), imagedata);
        xlabel('energy (eV)','fontsize',fontSize);    
    elseif strcmpi(xswitch,'pixels')
        subplot(2,1,1), plot(cax, imagedata);
        xlabel('pixel number');
    elseif strcmpi(xswitch,'GHz','fontsize',fontSize)
        subplot(2,1,1), plot(data.axisGHz(cax), imagedata);
        xlabel('f (GHz)','fontsize',fontSize);
    else 
        error('plot axis type not understood')
    end

    ylabel(ylab,'fontsize',fontSize);
    title(fname2);
    set(gcf,'PaperUnits','inches');
    xSize = 8;
    ySize = 10;
    axis tight;
    set(gcf,'PaperPosition',[0 -1 xSize ySize]); % Set output figure's paper position and size. [left bottom width height]
    %set(gcf,'Position',[0 0 xSize*50 ySize*50]) % set the size and location of the figure on the screen.

    %daspect('manual'),daspect([1 .3 1]); % set the aspect ratio of the final plot.
    box on

    % include the sif file characteristics on the plot.
    subplot(2,1,2);
    axis off;

    nf=55;

    d(1)={['File Name:      ' data.fileName(1:nf)]};
    d(2)={['                ' data.fileName((nf+1):length(data.fileName))]};
    d(3)={['Exposure Time:  ' num2str(data.exposureTime) ' seconds']};
    d(4)={['Grating:        ' num2str(data.grating)]};
    d(5)={['Center Wavelength of original spectrum: ' num2str(data.centerWavelength) ' nm' ]};
    d(6)={['Kinetic Length: ' num2str(data.kineticLength) ]};
    d(7)={['Cosmic rays removed from pixels: ' num2str(data.cosmicRays) ]};
    d(8)={['Spectrum cropped to pixels:  ' num2str(data.cropi(1)) ' - ' num2str(data.cropi(2)) ' ']};
    d(9)={['Date this plot was created: ' date]};

    text(0,.8,d);
    print('-dpdf',[fname(1:(length(fname)-4)) '.pdf']);

end


% FIGURE DISPLAYED IN MATLAB

if ~noplot           
    % Generate a plot in specified axes. Export plot handle.
    if axesHandleOn
        varargout{1} = plot(axesHandle,data.xaxis, imagedata,plotSpec);
    else
        varargout{1} = plot(data.xaxis, imagedata,plotSpec);
    end
    
    % Label y axis.
    ylabel(ylab,'fontsize',fontSize);
    if ischar(fname)
        title(fname2,'fontsize',fontSize);    
    else
        title('','fontsize',fontSize);
    end
    
    % Label x axis
    if xswitch(1)=='n'
        xlabel('\lambda (nm)','fontsize',fontSize);
    elseif xswitch(1)=='e'
        xlabel('Energy (eV)','fontsize',fontSize);
    elseif xswitch(1)=='p'
        xlabel('Pixel Number','fontsize',fontSize);
    elseif xswitch(1)=='G'
        xlabel('f (GHz)','fontsize',fontSize);
    else 
        error('plot axis type not understood')
    end
    
    axis tight;
    
    % Draw a box around the 'filter' limits.
    if filterOn
        ylimits = get(gca,'ylim');
        xlimits = get(gca,'xlim');
        filter = sort(filter);
        if filter(2)>xlimits(2)
            filter(2)=xlimits(2);
        end
        if filter(1)<xlimits(1)
            filter(1)=xlimits(1);
        end
        rectangle('Position',[min(filter),ylimits(1)+.001*abs(ylimits(2)-ylimits(1)),abs(filter(2)-filter(1)),.995*abs(ylimits(2)-ylimits(1))],...
            'LineWidth',1,...
            'LineStyle','--')
        %line([filter(1) filter(1)],ylimits,'color',[0 0 0]);
        %line([filter(2) filter(2)],ylimits,'color',[0 0 0]);
        set(gca,'layer','top')
    end
end

set(gca,'fontsize',fontSize)
axis tight;

if plotRangeOn
    ylim(plotRange)
end

% Output the manipulated data.
data.imageDataScaled = imagedata;

% Save the data as a mat file.
if saveData
    fnameLength = length(fname);
    fname3 = [fname(1:(fnameLength-4)) '.mat'];
    save(fname3,'data');
end

end


% -------------------------------------------------------------------------
% FUNCTION DEFINITIONS
% -------------------------------------------------------------------------


%REMOVE_UNDERSCORES remove underscores from a text string.
%
% Removes the underscores and last 4 letters from a file name. Useful for
% making a plot title from the file name.

function fprime = remove_underscores(fname)
fprime =fname;
for n=1:(length(fname)-4)
    if fname(n) == '_'
        fprime(n) = ' ';
    else
        fprime(n) = fname(n);
    end
end
end


%PUTINBOUDNS squeeze a value between a lower and upper bound
%
% output = putInBounds(input,lowBound,highBound) finds the element in the
% range [lowBound, highBound] that is closest to input, and returns it as
% output.
%
%
function output = putInBounds(input,lowBound,highBound)
if lowBound<=input & input<=highBound
    output = input;
elseif input<=lowBound;
    output = lowBound;
elseif input>=highBound;
    output = highBound;
end
end

%READOPTION : reads the value specified by an option in optionList.
%
% [ret, varargout] = readOption(option,numValues,defaultValue,optionList)
% looks for the char array option in the cell array optionList. If option
% is an element of optionList, then readOption reads the values of the next
% numValues elements of option List and sends them as varargout. Default
% Value is a cell array of the default values to set varargout if the
% option is not found.
%
% This all sounds rather complicated, but it is not. For example,
% calling readOption('fontsize',1,14,{'opt1','val1','fontsize',12}) with
% returns [1, 12], where readOption('founsize',1,14,{'opt1','val1'})
% returns [0, 14].
%
% defaultValue should be a cell array if it is more than one element.
% The program can read defaultValue being a single number as well.
%
% Todd Karin
% 06/28/2102

function [ret varargout] =  readOption(option,numValues,defaultValue,optionList)

ret = 0; % ret returns whether the option was found.

if iscell(defaultValue)
    varargout = defaultValue;
elseif isnumeric(defaultValue)&&length(defaultValue)==1 % If defaultValue is just one number, convert to a cell.
    for i=1:length(defaultValue)
        defaultValueCell{i} = defaultValue(i);
    end
    varargout = defaultValueCell; % return if no option is specified.
else
    varargout{1} = defaultValue;
end


if numValues >=1
    for i=1:length(optionList)
        if strcmpi(optionList{i},option)
            varargout = optionList((i+1):(i+numValues));
            ret=1;
        end
    end
elseif numValues ==0 % Print a 0 or 1 depedning on the presence of option
    for i=1:length(optionList)
        if strcmpi(optionList{i},option)
            ret = 1;
        end
    end
end

end


%FINDELEMENT Find the element of a vector closest to a value.
%
% element = findElement(x,v) returns the element of v that is closest to x.
%
% For example findElement(7.3,[5 6 7.4 7.5]) returns 3. If multiple
% elements of v are equadistant from x, then the first value found is
% returned.
%
% This function threads over lists and matrices, so for example
% findElement([7.3, 7.7],[5 6 7.4 7.5]) returns [3 4]
% 
% Todd Karin
% 06/26/2012

function element = findElement(x,v)

%[~,roundDirection] = readOption('round',1,'down',varargin);

for i=1:size(x,1)
    for j=1:size(x,2)
        [ret, indexes] = min(abs(v-x(i,j)));
        element(i,j) = indexes(1);
    end
end

end