function [] = PlotMetabolics(data,varargin)
    % PLOTMETABOLICS    Plot gas exchange data   
    %   Written by Grant Elliott, June 2011
    %
    %   PlotMetabolics(data) Generates a figure for a metabolic data structure produced by
    %   ProcessMetabolics.
    %
    %   PlotMetabolics(data, paramName, paramValue,...) takes additional optional parameters.
    %
    %   Accepts the following optional parameters:
    %   LineWidth                   - Width of VO2 and VCO2 plots in pixels (Defaults to 1)
    %   WindowLineWidth             - Width of VO2 and VCO2 plots within the analyzed Window in pixels (Defaults to 2 if available)
    %   DisplayGasExchangeRateMeans - Add numeric values of gas exchange rates to legend (Defaults to true if available)
    %   DisplayPower                - Add numeric value of metabolic power to legend (Defaults to true if available)
    %   DisplayNormalizedPower      - Add numeric value of normalized metabolic power to legend (Defaults to true if available)
    %   PlotGasExchangeRateMeans    - Plot horizontal lines for mean gas exchange rates (Defaults to true if available)
    %   HighlightForbidden          - Highlight forbidden regions (Defaults to false)
    %   O2Color                     - Color for plotting O2 data (Defaults to [1,0,0])
    %   CO2Color                    - Color for plotting CO2 data (Defaults to [0,0,1])
    %   ForbiddenColor              - Color for highlighting forbidden region (Defaults to [1,.7,.7])
    %   ErrorBarAlpha               - Alpha value to apply to highlighted error bars (Defaults to .3)
    %   YAxisMax                    - Max Scale of Gas Exchange Rate (Defaults to next 10mL/s)
    %   Axes                        - Axes in which to create plot (Defaults to axes in new figure)
    %   FigureSize                  - Size of new figure, if created. Also disables resize (Defaults to MatLab default; No effect if Axes is specified)
    %   AxesFontName                - Font to use for axes (Defaults to Helvetica)
    %   AxesFontSize                - Font size to use for axes (Defaults to 10)
    %   LabelFontName               - Font to use for axis labels (Defaults to Helvetica)
    %   LabelFontSize               - Font size to use for axis labels (Defaults to 10)   
    %   LegendFontName              - Font to use for legend labels (Defaults to Lucida Console)
    %   LegendFontSize              - Font size to use for legend labels (Defaults to 8)
    %   TitleText                   - Text to use in title (Defaults to description of raw data filename))
    %   TitleFontName               - Font to use for legend labels (Defaults to Lucida Console)
    %   TitleFontSize               - Font size to use for legend labels (Defaults to 10)
    %
    %   See also ProcessMetabolics, Measurement.
        
    if(~(isstruct(data) && isfield(data,'Data') && isstruct(data.Data) && isfield(data.Data,'VO2') && isfield(data.Data,'VCO2')))
        error('Metabolics:PlotMetabolics:BadFormat','Data structure does not contain VO2 and VCO2 data');
    end
    
    fileNameDefined = isfield(data,'FileName') && ~isempty(data.FileName);
    windowDefined = isfield(data,'Window') && isstruct(data.Window) && isfield(data.Window,'WindowIndices');
    meansDefined = isfield(data,'Results') && isstruct(data.Results) && isfield(data.Results,'VO2') && isfield(data.Results,'VCO2');
    powerDefined = isfield(data,'Results') && isstruct(data.Results) && isfield(data.Results,'Power');
    normalizedPowerDefined = isfield(data,'Results') && isstruct(data.Results) && isfield(data.Results,'NormalizedPower');
    forbiddenDefined = isfield(data.Data,'GlitchIndices') && isfield(data.Data,'HoldIndices') && isfield(data.Data,'LockoutStart') && isfield(data.Data,'LockoutStop') && isfield(data.Data,'LockoutHold');

    defaultTitle = '';
    if (fileNameDefined)
        [~,f,e]=fileparts(data.FileName);
        defaultTitle=['Metabolic Analysis of File: ',f,e];
    end
    
    p=inputParser;
    p.addParamValue('LineWidth', 1, @(x) isnumeric(x) && numel(x)==1 && x>0);
    p.addParamValue('WindowLineWidth', 2, @(x) isnumeric(x) && numel(x)==1 && x>0);
    p.addParamValue('DisplayGasExchangeRateMeans', meansDefined , @(x) islogical(x) && numel(x)==1 && (~x || meansDefined));
    p.addParamValue('DisplayPower', powerDefined , @(x) islogical(x) && numel(x)==1 && (~x || powerDefined));
    p.addParamValue('DisplayNormalizedPower', normalizedPowerDefined , @(x) islogical(x) && numel(x)==1 && (~x || powerDefined));
    p.addParamValue('PlotGasExchangeRateMeans', meansDefined , @(x) islogical(x) && numel(x)==1 && (~x || meansDefined)); 
    p.addParamValue('HighlightForbidden', false , @(x) islogical(x) && numel(x)==1 && (~x || forbiddenDefined));
    p.addParamValue('O2Color',[1,0,0], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
    p.addParamValue('CO2Color',[0,0,1], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
    p.addParamValue('ForbiddenColor',[1,.7,.7], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
    p.addParamValue('ErrorBarAlpha',.3, @(x) isnumeric(x) && numel(x)==1 && x<=1 && x>=0);
    p.addParamValue('Axes', NaN, @(x) ishandle(x) && strcmp(get(x,'type'),'axes'));
    p.addParamValue('YAxisMax', NaN, @(x) isnumeric(x) && numel(x)==1 && x>0);
    p.addParamValue('FigureSize', NaN, @(x) isnumeric(x) && numel(x)==2 && min(x)>0);
    p.addParamValue('AxesFontName', 'Helvetica', @(x) ischar(x) && any(strcmp(x,listfonts())));
    p.addParamValue('AxesFontSize', 10, @(x) isnumeric(x) && numel(x)==1 && x>0);
    p.addParamValue('LabelFontName', 'Helvetica', @(x) ischar(x) && any(strcmp(x,listfonts())));
    p.addParamValue('LabelFontSize', 10, @(x) isnumeric(x) && numel(x)==1 && x>0);
    p.addParamValue('LegendFontName', 'Lucida Console', @(x) ischar(x) && any(strcmp(x,listfonts())));
    p.addParamValue('LegendFontSize', 8, @(x) isnumeric(x) && numel(x)==1 && x>0);
    p.addParamValue('TitleText', defaultTitle, @(x) ischar(x));
    p.addParamValue('TitleFontName', 'Helvetica', @(x) ischar(x) && any(strcmp(x,listfonts())));
    p.addParamValue('TitleFontSize', 12, @(x) isnumeric(x) && numel(x)==1 && x>0);
    p.parse(varargin{:});
    settings = p.Results;

    if (isnan(settings.Axes))
        figure();
        set(gcf,'Color','White');
        if (~isnan(settings.FigureSize))
            fs1=get(gcf,'outerposition');
            fs2=get(gcf,'position');
            ss=get(0,'screensize');
            fp=min([fs2(1:2);ss(3:4)-settings.FigureSize-(fs1(1:2)+fs1(3:4)-fs2(1:2)-fs2(3:4))]);
            set(gcf,'position',[fp,settings.FigureSize],'resize','off');
        end
        settings.Axes=gca;
    end
    cla(settings.Axes);
    
    plotTime = data.Data.TimeResampled;
    plotVO2 = data.Data.VO2*1000;
    plotVCO2 = data.Data.VCO2*1000;
    plotTime(data.Data.GlitchIndices) = nan;
    plotVO2(data.Data.GlitchIndices) = nan;
    plotVCO2(data.Data.GlitchIndices) = nan;
    if (windowDefined)
        windowIndices = data.Window.WindowIndices;
        plotWindowTime=nan(size(plotTime));
        plotWindowVO2=nan(size(plotVO2));
        plotWindowVCO2=nan(size(plotVCO2));
        plotWindowTime(windowIndices)=plotTime(windowIndices);
        plotWindowVO2(windowIndices)=plotVO2(windowIndices);
        plotWindowVCO2(windowIndices)=plotVCO2(windowIndices);
    end
    for i=fliplr(data.Data.HoldIndices)
        plotTime = [plotTime(1:i-1),NaN,plotTime(i:end)];
        plotVO2 = [plotVO2(1:i-1),NaN,plotVO2(i:end)];
        plotVCO2 = [plotVCO2(1:i-1),NaN,plotVCO2(i:end)];
        if (windowDefined)
            plotWindowTime = [plotWindowTime(1:i-1),NaN,plotWindowTime(i:end)];
            plotWindowVO2 = [plotWindowVO2(1:i-1),NaN,plotWindowVO2(i:end)];
            plotWindowVCO2 = [plotWindowVCO2(1:i-1),NaN,plotWindowVCO2(i:end)];
        end
    end
    
    hold on;
    if (isnan(settings.YAxisMax))
        settings.YAxisMax=ceil(max([max(plotVO2),max(plotVCO2)])/10)*10;
    end
    axis([0,plotTime(find(isnan(plotTime)==0,1,'last')),0,settings.YAxisMax]);

    if (settings.HighlightForbidden)
        for i=data.Data.GlitchIndices
            rect(data.Data.TimeResampled(i-1),0,data.Data.TimeResampled(i+1),settings.YAxisMax,'none',settings.ForbiddenColor,1);
        end
        for i=data.Data.HoldIndices
            rect(data.Data.TimeResampled(i-1),0,data.Data.TimeResampled(i)+data.Data.LockoutHold,settings.YAxisMax,'none',settings.ForbiddenColor,1);
        end
        if (data.LockoutStart~=0)
            rect(0,0,data.Data.LockoutStart,settings.YAxisMax,'none',settings.ForbiddenColor,1);
        end
        if (data.LockoutStop~=0)
            rect(data.Data.TimeResampled(end)-data.Data.LockoutStop,0,data.Data.TimeResampled(end),settings.YAxisMax,'none',settings.ForbiddenColor,1);
        end
    end
    
    if (settings.PlotGasExchangeRateMeans)
        Measurement.plotFunctionDistribution([0,plotTime(end)],1000*data.Results.VO2*[1,1],'lineColor',settings.O2Color,'borderLineColor',settings.O2Color,'fillColor',settings.O2Color,'fillAlpha',settings.ErrorBarAlpha);
        Measurement.plotFunctionDistribution([0,plotTime(end)],1000*data.Results.VCO2*[1,1],'lineColor',settings.CO2Color,'borderLineColor',settings.CO2Color,'fillColor',settings.CO2Color,'fillAlpha',settings.ErrorBarAlpha);
    end
    
    O2Plot = plot(plotTime,plotVO2,'color',settings.O2Color,'linewidth',settings.LineWidth);
    CO2Plot = plot(plotTime,plotVCO2,'color',settings.CO2Color,'linewidth',settings.LineWidth);
    if (windowDefined && settings.WindowLineWidth~=settings.LineWidth)
        O2Plot = plot(plotWindowTime,plotWindowVO2,'color',settings.O2Color,'linewidth',settings.WindowLineWidth);
        CO2Plot = plot(plotWindowTime,plotWindowVCO2,'color',settings.CO2Color,'linewidth',settings.WindowLineWidth);
    end
    legp = [O2Plot, CO2Plot];
    legt = {'O_2','CO_2'};
    if (settings.DisplayPower)
        legp(length(legp)+1)=gca();
        legt{length(legt)+1}=sprintf('Power     % 5.f    %s% 3.f    W',data.Results.Power.Value,'\pm',data.Results.Power.Uncertainty);
    end
    if (settings.DisplayNormalizedPower)
        legp(length(legp)+1)=gca();
        legt{length(legt)+1}=sprintf('Power/Mass  % 5.2f %s % 5.2f W/kg',data.Results.NormalizedPower.Value,'\pm',data.Results.NormalizedPower.Uncertainty);
    end
    if (settings.DisplayGasExchangeRateMeans)
        legt{1}=sprintf('O_2         % 7.2f %s % 4.2f mL/s',data.Results.VO2.Value*1000,'\pm',data.Results.VO2.Uncertainty*1000);
        legt{2}=sprintf('CO_2        % 7.2f %s % 4.2f mL/s',data.Results.VCO2.Value*1000,'\pm',data.Results.VCO2.Uncertainty*1000);
    end
    set(legend(legp,legt),'fontname',settings.LegendFontName,'fontsize',settings.LegendFontSize,'location','southeast');
    ylabel('Gas Exchange Rate (mL/s)','fontname',settings.LabelFontName,'fontsize',settings.LabelFontSize);
    xlabel('Time (s)','fontname',settings.LabelFontName,'fontsize',settings.LabelFontSize);
    set(gca,'fontname',settings.AxesFontName,'fontsize',settings.AxesFontSize,'xtick',0:60:max(plotTime));
    if (~isempty(settings.TitleText))
        title(settings.TitleText,'fontname',settings.TitleFontName,'fontsize',settings.TitleFontSize);
    end
    set(gca,'layer','top')
    hold off;
end

function [h]=rect(x1,y1,x2,y2,e,f,a)
    h=patch([x1,x1,x2,x2],[y1,y2,y2,y1],[1,1,1]*(1-a)+f*a,'edgecolor',e);
end