%Copyright (c) 2010, Luuk J.G.W. van Wilderen
%
%plot_figs.m

function plot_figs
%Main plotting script.
%The figures that are generated contain:
%-The raw data (in 3D), with the fit in continous color through it (patches
%represent sections where the fit lies above or below the data), with the 
%concentration curves below it. Figure name: Raw data Concentration profiles.
%-A maximum of 200 time traces, in blocks of 40 (if the data set contains 
%more, these are not shown, although they are fitted).  Figure names:
%trace1, trace2, etc... 
%-A surface plot of the fitted data, and slices through the data can be 
%visualised by clicking two times in the surface plot (the first one 
%activates the plot, the second one selects a data point).  Figure name:
%Explore fit. 
%-The spectra corresponding to the components.  Figure name: spectra.
%-The residuals if selected in the gui.  Figure name: residuals

%Defined in globfit_bio_gui
global nr_waves nr_svd_comps xdata wave y_fit ydata start_pars
global RESNORM RESIDUAL opt_pars obj opt_specs opt_times time_fix_counter

%Defined in fit_data
global SVDcontrib_switch axis_switch S V
global time_limits wave_limits start_par_time_zero
global filename comps analysis IRF Par_fix_time_zero
global Error_switch_lin_spectra plot_residuals Par_fix_IRF
%Previously defined, if doesn't exist yet (f.i. if 'Replot fit') is
%pressed, then it is defined here.
global lin_log

%Defined in glob_fun_bio
global decays spectra_std Par_IRF IRF_fun

%Defined here
global svd_comp_perc

% %Color order for graphs
Color_matrix=[  0   0   1.0;     
    1.0 0   0;
    0   1.0 0;
    0   1.0 1.0;
    1.0 0   1.0;
    0   0   0
    1   1   0];

%Plotting
%Time axis toggle. If not switched on, all time axes are logarithmic.
if exist('lin_log')==0 || isempty('lin_log')==1
    lin_log=get(findobj('Tag','lin_log_toggle'),'Value');
end

%Disable 'negative timepoints' warning if logarithmic time axis used. Is
%switched back on at end of this file
warning('off','MATLAB:Axes:NegativeDataInLogAxis')

%Showing raw data
plot_raw_data
set(findobj('Name','Loaded_raw_data'),'WindowStyle','docked'); %Create single figure window that contains all following plots as tabs inside the 'Figures Window'
%Raw data explorer figure is deleted because same data present in fit
%explorer plot
if isempty(findobj('Name','Explore raw data'))==0
    close ('Explore raw data');
end

%----Plotting all traces in one plot and concentration profiles----
raw_conc_prof=figure;
set(raw_conc_prof,'name','Raw data Concentration profiles','Units','normalized','position',[0.02 0.05 0.48 0.8],...
    'NumberTitle','off','WindowStyle','docked');
subplot(2,1,1)
if nr_svd_comps==0
    surf(xdata,wave,y_fit')
    shading interp;
    hold on;
    surf(xdata,wave,ydata')
    title(['SSE: ',num2str(RESNORM),'. Grid is data; continuous colors are fit']);
    axis tight;
    
    if lin_log==0
        set(gca,'XScale','log');
    end
    if axis_switch==0 && nr_svd_comps==0 && SVDcontrib_switch==0
        set(gca,'YDir','reverse');
    end
elseif nr_svd_comps>0
    title(['Check the seperate spectra figure']);
end

%----Plotting traces individually----
trace1=figure;
set(trace1,'name','trace1','Units','normalized','position',[0.02 0.05 0.97 0.8],...
    'NumberTitle','off','WindowStyle','docked');
wave_sel=1:nr_waves;
nr_plots_row=8;         %Number of plots per row
if nr_waves<nr_plots_row
    nr_plots_row=nr_waves+1;
end
if nr_svd_comps>0   %For SVD analysis also show singular values
    for t=1:nr_waves
        subplot(ceil(nr_waves/nr_plots_row)+1,nr_plots_row,t)
        if lin_log==0
            semilogx(xdata,y_fit(:,wave_sel(t)))
            hold on;
            semilogx(xdata,ydata(:,wave_sel(t)),'k')
        elseif lin_log==1
            plot(xdata,y_fit(:,wave_sel(t)))
            hold on;
            plot(xdata,ydata(:,wave_sel(t)),'k')
        end
        axis tight;
        if t==1
            legend('fit','LSV')
        end
        title(['SV = ',num2str(t)])
        subplot(ceil(nr_waves/nr_plots_row)+1,nr_plots_row,t+nr_plots_row)
        if SVDcontrib_switch==0
            plot(wave,V(:,t),'k')
        elseif SVDcontrib_switch==1
            for i=1:nr_svd_comps
                SV(:,i)=S(i,i)*V(:,i);
            end
            plot(wave,SV(:,t),'k')
        end
        axis tight;
        title(['SV = ',num2str(t)])
    end
    clear t
    subplot(ceil(nr_waves/nr_plots_row)+1,nr_plots_row,1)
    ylabel('LSV');
    subplot(ceil(nr_waves/nr_plots_row)+1,nr_plots_row,1+nr_plots_row)
    ylabel('RSV');
    for t=nr_waves+1
        subplot(ceil(nr_waves/nr_plots_row)+1,nr_plots_row,t)
        plot(diag(S),'r*-')
        title('Singular values (SV)');
        axis tight;
    end
    clear t
    if SVDcontrib_switch==0
        annotation('textbox',[0.4 0.9 .1 .1],'string','SVD analysis of selected number of components',...
            'linestyle','none');
    elseif SVDcontrib_switch==1
        annotation('textbox',[0.4 0.9 .1 .1],'string','Trace (LSV)*Singular value (S)','linestyle','none');
        annotation('textbox',[0.4 0.05 .1 .01],'string','Singular value (S)*Spectra (RSV)','linestyle','none');
    end
elseif nr_svd_comps==0
    wave_digits=get(findobj('Tag','wave_digits_box'),'String'); %Only for displaying purposes
    plot_wave_divider=str2num(get(findobj('Tag','plot_wave_divider_box'),'String')); %Only for displaying purposes
    if plot_wave_divider==0;
        plot_wave_divider=1;
    end
    precision_wavelengths=['%10.',wave_digits,'f'];
    clear wave_digits
    %Depending on the amount of traces to be plotted, this section
    %determines where to plot them. One page (one big figure) can only 
    %contain 40 traces (8 x 5 small figures), and a maximum of 6 pages is
    %shown. The plot_wave_divider parameter is read from the gui and is
    %only for plotting. All traces are fit, but not all are shown.
    for t=1:nr_waves
        if t<41  && t*plot_wave_divider<nr_waves
            if nr_waves<41
                subplot(floor((nr_waves-1)/nr_plots_row)+1,nr_plots_row,t)
            elseif nr_waves>40
                subplot(5,nr_plots_row,t)
            end
            if lin_log==0
                semilogx(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                semilogx(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            elseif lin_log==1
                plot(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                plot(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            end
            axis tight;
            %grid on;
        elseif nr_waves>40 && t>40 && t<81 && t*plot_wave_divider<nr_waves
            if t==41
                trace2=figure;
                set(trace2,'name','trace2','Units','normalized','position',[0.02 0.05 0.97 0.8],...
                    'NumberTitle','off','WindowStyle','docked');
            end
            if nr_waves<81
                subplot(floor((nr_waves-41)/nr_plots_row)+1,nr_plots_row,t-40)
            elseif nr_waves>80
                subplot(5,nr_plots_row,t-40)
            end
            if lin_log==0
                semilogx(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                semilogx(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            elseif lin_log==1
                plot(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                plot(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            end
            axis tight;
            %grid on;
        elseif nr_waves>80 && t>80 && t<121 && t*plot_wave_divider<nr_waves
            if t==81
                trace3=figure;
                set(trace3,'name','trace3','Units','normalized','position',[0.02 0.05 0.97 0.8],...
                    'NumberTitle','off','WindowStyle','docked');
            end
            if nr_waves<121
                subplot(floor((nr_waves-81)/nr_plots_row)+1,nr_plots_row,t-80)
            elseif nr_waves>120
                subplot(5,nr_plots_row,t-80)
            end
            if lin_log==0
                semilogx(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                semilogx(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            elseif lin_log==1
                plot(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                plot(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            end
            axis tight;
        elseif nr_waves>120 && t>120 && t<161 && t*plot_wave_divider<nr_waves
            if t==121
                trace4=figure;
                set(trace4,'name','trace4','Units','normalized','position',[0.02 0.05 0.97 0.8],...
                    'NumberTitle','off','WindowStyle','docked');
            end
            if nr_waves<161
                subplot(floor((nr_waves-121)/nr_plots_row)+1,nr_plots_row,t-120)
            elseif nr_waves>160
                subplot(5,nr_plots_row,t-120)
            end
            if lin_log==0
                semilogx(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                semilogx(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            elseif lin_log==1
                plot(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                plot(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            end
            axis tight;
            %grid on;
        elseif nr_waves>160 && t>160 && t<201 && t*plot_wave_divider<nr_waves
            if t==161
                trace5=figure;
                set(trace5,'name','trace5','Units','normalized','position',[0.02 0.05 0.97 0.8],...
                    'NumberTitle','off','WindowStyle','docked');
            end
            if nr_waves<201
                subplot(floor((nr_waves-161)/nr_plots_row)+1,nr_plots_row,t-160)
            elseif nr_waves>200
                subplot(5,nr_plots_row,t-160)
            end
            if lin_log==0
                semilogx(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                semilogx(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            elseif lin_log==1
                plot(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                plot(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            end
            axis tight;
            %grid on;
        elseif nr_waves>200 && t>200 && t<241 && t*plot_wave_divider<nr_waves
            if t==201
                trace6=figure;
                set(trace6,'name','trace6','Units','normalized','position',[0.02 0.05 0.97 0.8],...
                    'NumberTitle','off','WindowStyle','docked');
            end
            if nr_waves<241
                subplot(floor((nr_waves-201)/nr_plots_row)+1,nr_plots_row,t-200)
            elseif nr_waves>240
                subplot(5,nr_plots_row,t-200)
            end
            if lin_log==0
                semilogx(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                semilogx(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            elseif lin_log==1
                plot(xdata,y_fit(:,wave_sel(t*plot_wave_divider)))
                hold on;
                title(num2str(wave(wave_sel(t*plot_wave_divider)),precision_wavelengths));
                plot(xdata,ydata(:,wave_sel(t*plot_wave_divider)),'k')
            end
            axis tight;
            %grid on;
        end
    end
    clear t wave_sel precision_wavelengths nr_plots_row
end

%------Interactive plot that shows raw data and its fits----
if nr_svd_comps==0 %Plot is not created during SVD analysis
    fit_explorer=figure;
    set(fit_explorer,'name','Explore fit','Units','normalized','position',[0.02 0.05 0.97 0.8],...
        'NumberTitle','off','Toolbar','figure','WindowStyle','docked',...
        'WindowButtonDownFcn','fit_explorer_plot');
    %Renderer setting below needed to prevent cross-hair from drawing over figure
    if lin_log==1
        set(fit_explorer,'renderer','opengl')
    end
    %Raw data surface plot
    subplot(3,3,[1:2,4:5])
    if time_limits==0
        time_limits(1)=1;
        time_limits(2)=size(xdata,1);
    end
    if wave_limits==0
        wave_limits(1)=1;
        wave_limits(2)=size(wave,2);
    end
    pcolor(xdata,wave,ydata');%(time_limits(1):time_limits(2),wave_limits(1):wave_limits(2))')
    shading interp;
    colorbar
    title({'Click in surface plot to activate cross-hair, click again for slice selection',['Datafile: ' filename]},'Interpreter','None');    %The TEX Interpreter is switched off to prevent error message due to possible interpretation of forward slashes
    xlabel('Time');
    ylabel('Wavelength');
    hold on;
    if lin_log==0
        set(gca,'XScale','log');
    end
    set(gca,'YDir','reverse');
    if axis_switch==1 && nr_svd_comps==0 && SVDcontrib_switch==0
        set(gca,'YDir','normal');
    end
    
    %Time trace plot
    subplot(3,3,7:8)
    plot(xdata,ydata(:,1),'k');%time_limits(1):time_limits(2),1))
    hold on;
    plot(xdata,y_fit(:,1),'b');%time_limits(1):time_limits(2),1),'r')
    if lin_log==0
        set(gca,'XScale','log');
    end
    xlabel(['Time trace at wavelength ',num2str(wave(1))]);
    axis tight;
    grid on;
    
    %Spectrum plot
    subplot(3,3,[3,6])
    plot(wave,ydata(1,:),'k');%wave_limits(1):wave_limits(2)))
    hold on;
    plot(wave,y_fit(1,:),'b');%wave_limits(1):wave_limits(2)),'r')
    
    xlabel(['Spectrum at timepoint ',num2str(xdata(1))]);
    axis tight;
    if axis_switch==1 && nr_svd_comps==0 && SVDcontrib_switch==0
        set(gca,'XDir','reverse');
    end
    legend('Raw data','Fit');
    set(findobj('Tag','legend'),'Position',[0.75 0.15 0.0929 0.0942]);
    grid on;
end

%----Plotting concentration profiles----
%Could be >1 at t=0 due to fitted time zero
figure(raw_conc_prof);
subplot(2,1,2)
clr='brgcmkybrgcmky';
trg_tst=1;  %Checks if target model is loaded that has equal amount of species and reactions
if size(obj.species,1)-1<size(obj.reactions,1)
    trg_tst=0;
end
for c=1:size(obj.Species,1)-trg_tst
    if lin_log==0
        semilogx(xdata,decays(:,c),clr(c));
    elseif lin_log==1
        plot(xdata,decays(:,c),clr(c));
    end
        hold on;
end
clear c
axis tight;
if length(IRF)==2  %Plot also instrument response function
    if lin_log==0
        semilogx(xdata,IRF_fun,'k--');
    elseif lin_log==1
        plot(xdata,IRF_fun,'k--');
    end
    hold on
end
set(gca,'XLim',[start_par_time_zero max(xdata)],'YLim',[0  1]);

clear legend_text_old legend_text_new;
rev_counter=0;      %Count number of reversible reactions
for c=1:size(obj.Reactions,1)
    if c==1
        if obj.Reactions(c).Reversible==0
            legend_text_old=num2str(opt_times(c));
            xlabels_old=num2str(round(opt_times(c)));
        elseif obj.Reactions(c).Reversible==1
            legend_text_old=strcat(num2str(opt_times(c)),'; ',num2str(opt_times(c+1)));
            xlabels_old=strcat(num2str(round(opt_times(c))),'; ',num2str(round(opt_times(c+1))));
            rev_counter=rev_counter+1;
        end
    elseif c>1
        if obj.Reactions(c).Reversible==0
            legend_text_new=strvcat(legend_text_old,num2str(opt_times(c+rev_counter)));
            xlabels_new=strvcat(xlabels_old,num2str(round(opt_times(c+rev_counter))));
            legend_text_old=legend_text_new;
            xlabels_old=xlabels_new;
        elseif obj.Reactions(c).Reversible==1
            legend_text_new=strvcat(legend_text_old,strcat(num2str(opt_times(c+rev_counter)),'; ',...
                num2str(opt_times(c+1+rev_counter))));
            xlabels_new=strvcat(xlabels_old,strcat(num2str(round(opt_times(c+rev_counter))),'; ',...
                num2str(round(opt_times(c+1+rev_counter)))));
            legend_text_old=legend_text_new;
            xlabels_old=xlabels_new;
            rev_counter=rev_counter+1;
        end
    end
end
clear legend_text legend_text_new c
if size(obj.species,1)-1>size(obj.reactions,1)
    legend_text=strvcat(legend_text_old,'End product');
    xlabels=strvcat(xlabels_old,'End product');
elseif size(obj.species,1)-1==size(obj.reactions,1)
    legend_text=legend_text_old;
    xlabels=xlabels_old;
elseif size(obj.species,1)-1<size(obj.reactions,1)
    legend_text=legend_text_old;
    xlabels=xlabels_old;
end
clear legend_text_old xlabels_old 
if length(IRF)==2 && Par_fix_IRF==0 
    legend_text=strvcat(legend_text,['IRF (\Delta=',num2str(opt_pars(comps+2-time_fix_counter-Par_fix_time_zero)),')']);
elseif (Par_fix_IRF==1 && length(IRF)==2)
    legend_text=strvcat(legend_text,['IRF (\Delta=',num2str(start_pars(comps+2)),')']);
end
legend (legend_text)
title('Concentration profiles')

%----Plotting difference spectra of sequential, parallel and target models---
spectra_fig=figure;
set(spectra_fig,'Toolbar','figure','name','spectra','Units','normalized','position',[0.515 0.05 0.475 0.8],...
    'NumberTitle','off','WindowStyle','docked');
set(gca,'colororder',Color_matrix,'NextPlot','replacechildren');
if nr_svd_comps>0
    if SVDcontrib_switch==0
        if Error_switch_lin_spectra==0
            if Par_IRF==0
                plot(1:nr_svd_comps,opt_specs);
            elseif Par_IRF==1
                 plot(1:nr_svd_comps,opt_specs(:,1:comps));
                 hold on
                 plot(1:nr_svd_comps,opt_specs(:,end),'k--');
            end
        elseif Error_switch_lin_spectra==1
            for i=1:comps
                errorbar(1:nr_svd_comps,opt_specs(:,i),spectra_std(i,:),clr(i));
                hold on;
            end
            clear i
            if Par_IRF==1
                errorbar(1:nr_svd_comps,opt_specs(:,end),spectra_std(end,:),'k--');
            end
        end
        set(gca,'xtick',[1:nr_svd_comps],'Xticklabel',num2str([1:nr_svd_comps]'));
        xlabel('SVD component number')
    elseif SVDcontrib_switch==1
        subplot(1,2,1)
        set(gca,'colororder',Color_matrix,'NextPlot','replacechildren');
        if Error_switch_lin_spectra==0
            if Par_IRF==0
                plot(1:size(opt_specs,1),opt_specs);
            elseif Par_IRF==1
                plot(1:size(opt_specs,1),opt_specs(:,1:comps));
                hold on
                plot(1:size(opt_specs,1),opt_specs(:,end),'k--');
            end
        elseif Error_switch_lin_spectra==1
            for i=1:size(obj.Species,1)-1
                errorbar(1:nr_svd_comps,opt_specs(:,i),spectra_std(i,:),clr(i));
                hold on;
            end
            clear i
            if Par_IRF==1
                errorbar(1:nr_svd_comps,opt_specs(:,end),spectra_std(end,:),'k--');
            end
        end
        set(gca,'xtick',[1:nr_svd_comps],'Xticklabel',num2str([1:nr_svd_comps]'));
        xlabel('SVD component number')
        grid on;
        axis tight;
        subplot(1,2,2)
        svd_comp_perc=100*sum(abs(opt_specs))/sum(sum(abs(opt_specs)));
        plot(1:size(obj.Species,1)-1+Par_IRF-rev_counter,svd_comp_perc);
        if Par_IRF==0
            if size(obj.species,1)-1>size(obj.reactions,1)
                set(gca,'xtick',1:size(obj.Species,1)-1,'Xticklabel',strvcat(num2str(round(opt_times')),'End product'));
            elseif size(obj.species,1)-1==size(obj.reactions,1)
                set(gca,'xtick',1:size(obj.Species,1)-1,'Xticklabel',xlabels);
            end
        elseif Par_IRF==1
            if size(obj.species,1)-1>size(obj.reactions,1)
                set(gca,'xtick',1:size(obj.Species,1)-1+Par_IRF-rev_counter,'Xticklabel',strvcat(num2str(round(opt_times')),'End product','IRF'));
            elseif size(obj.species,1)-1==size(obj.reactions,1)
                set(gca,'xtick',1:size(obj.Species,1)-1+Par_IRF-rev_counter,'Xticklabel',strvcat(xlabels,'IRF'));
            end
        end            
        clear xlabels
    end
elseif nr_svd_comps==0
    if Error_switch_lin_spectra==0
        if Par_IRF==0
            plot(wave,opt_specs);
        elseif Par_IRF==1
            plot(wave,opt_specs(:,1:comps));
            hold on
            plot(wave,opt_specs(:,end),'k--')
        end
        
    elseif Error_switch_lin_spectra==1
        for i=1:size(obj.Species,1)-trg_tst
            errorbar(wave,opt_specs(:,i),spectra_std(i,:),clr(i));
            hold on
        end
        if Par_IRF==1
            errorbar(wave,opt_specs(:,end),spectra_std(end,:),'k--');
        end
        clear i
    end
end
clear clr
if axis_switch==1 && nr_svd_comps==0 && SVDcontrib_switch==0
    set(gca,'XDir','reverse');
end
axis tight;
grid on;
if SVDcontrib_switch==0 && nr_svd_comps==0
    if analysis==1
        title('Parallel model');
    elseif analysis==2
        title('Sequential model');
    elseif analysis==3
        title('Target model');
    end
    legend (legend_text)
elseif SVDcontrib_switch==0 && nr_svd_comps>0
    if analysis==1
        title('Parallel model: Amplitude per component');
    elseif analysis==2
        title('Sequential model: Amplitude per component');
    elseif analysis==3
        title('Target model: Amplitude per component');
    end
    legend (legend_text)
elseif SVDcontrib_switch==1 && nr_svd_comps>0
    subplot(1,2,1)
    if analysis==1
        title('Parallel model (LSV*S): Amplitude per component');
    elseif analysis==2
        title('Sequential model (LSV*S): Amplitude per component');
    elseif analysis==3
        title('Target model (LSV*S): Amplitude per component');
    end
    legend (legend_text)
    subplot(1,2,2)
    title('Percentual contribution per species')
    xlabel('Decay time (ps)')
end
clear Color_matrix S V legend_text
if nr_svd_comps>0
    figure(trace1); %Bring traces to front
end

%----Plotting residuals----
if plot_residuals==1
    residuals=figure;
    set(residuals,'Toolbar','figure','name','residuals','Units','normalized','position',[0.515 0.05 0.475 0.8],...
        'NumberTitle','off','WindowStyle','docked');
    subplot(2,1,1)
    if lin_log==0
        semilogx(xdata,sum(RESIDUAL,2),'r');
    elseif  lin_log==1
        plot(xdata,sum(RESIDUAL,2),'r');
    end
    xlabel('Time')
    axis tight;
    grid on;
    title(['Sum of residual analysis in time and spectrally. SSE: ',num2str(RESNORM),'.']);
    subplot(2,1,2)
    if nr_svd_comps==0
        plot(wave,sum(RESIDUAL,1),'r');
        xlabel('Wavelength')
    elseif nr_svd_comps>0
        plot(1:nr_svd_comps,sum(RESIDUAL,1),'r');
        xlabel('SVD component number')
        set(gca,'xtick',[1:nr_svd_comps],'Xticklabel',num2str([1:nr_svd_comps]'));
    end
    if axis_switch==1 && nr_svd_comps==0 && SVDcontrib_switch==0
        set(gca,'XDir','reverse');
    end
    axis tight;
    grid on;
end

%Turn warning about negative timepoints back on again (see above)
warning('on','MATLAB:Axes:NegativeDataInLogAxis')