function spectra(varargin)
% spectra
% function spectra(spectra,legend titles,units)
%
% Function that plots spectra and their cummulative spectra. The maximum
% number of spectra is three. The legend title and the legend titles and
% units are optional
%
% Examples of use:
%
% [x1,fx1]=pwelch(y1,2^12,2^11,2^12,25)
% [x2,fx2]=pwelch(y2,2^12,2^11,2^12,25)
% [x3,fx3]=pwelch(y3,2^12,2^11,2^12,25)
% spectra(fx1,x1,fx2,x2,fx3,x3);
% spectra(fx1,x1,fx2,x2,fx3,x3,'Case a','Case b','Case c');
% spectra(fx1,x1,fx2,x2,fx3,x3,'Case a','Case b','Case c','Nm');

%   Project: IMPACT
%   Title of Code: setSignalSelection
%   Relevant additional documentation:
%   Developer: Carlos Gonzalez
%   External Dependencies: NA
%
%   Disclaimer
%   CAUTION: It is strongly recommended that changes to this code
%   should be avoided without first consulting with authors.

spect = 0;
units = '';

for i=1:nargin
    if ~isnumeric(varargin{i})
        break
    end
    spect=spect+1;
end

if mod(spect,2)~=0
    error('Wrong number of input arguments')
end

% if (spect/2>3)
%     error('The maximum number of spectra is three')
% end

for i=spect+1:nargin
    if ~ischar(varargin{i})
        error('Last arguments must be char type');
    end
end

if spect==nargin
    for i=1:spect/2
        eval(['string',num2str(i),'=[''Spectrum '',num2str(i)];']);
    end
else
    if (nargin==spect*3/2)
        for i=1:spect/2
            eval(['string',num2str(i),'=varargin{spect+i};']);
        end
    else
        if (nargin==spect*3/2+1)
            for i=1:spect/2
                eval(['string',num2str(i),'=varargin{spect+i};']);
            end
            units=char(varargin(nargin));
        else
            error('Wrong number of input arguments');
        end
    end
end

spect=spect/2;


l_style1={'-','-','-','-','-','-','-'};
%color={'0 0 0','0.4 0.4 0.4','0.7 0.7 0.7','1 0 0','0 1 0','0 0 1','0 1 1'};
color={'0 0 1','0 0.5 0','1 0 0',' 0 0.75 0.75','0.75 0 0.75 ','0.75 0.75 0','0.25 0.25 0.25'};



% spect is the number of spectra to plot
for i=1:spect
    fx=varargin{2*(i-1)+1};     % Frequency
    sx=varargin{2*i};           % Spectra
    % [axn,h1n,h2n]             % n = spect
    % axn -> axes =[ax1,ax2]
    % h1n -> handle to the spectra plot
    % h2n -> handle to the cumulative spectra plot
    eval(['[ax',num2str(i),',h1',num2str(i),',h2',num2str(i),']=plotyy(2*pi*fx,sx,2*pi*fx,cumsum(sx),''loglog'',''semilogx'');']);
    %     eval(['[ax',num2str(i),',h1',num2str(i),',h2',num2str(i),']=plotyy(2*pi*fx,sx,2*pi*fx,cumsum(2*pi*fx.*sx),''loglog'',''semilogx'');']);
    eval(['set(h1',num2str(i),',''linestyle'',''',l_style1{mod(i-1,7)+1},''',''color'',[',color{mod(i-1,7)+1},'],''linewidth'',1);']);
    eval(['set(h2',num2str(i),',''linestyle'',''--'',''color'',[',color{mod(i-1,7)+1},'],''linewidth'',1);']);
    hold on
end

% Legend
% h will contain all the handles to the spectra plots
h=[];
st=[];


if spect~=1
    for i=1:spect
        if i==1
            h='[h11';
        elseif i==spect
            h=[h,[',h1',num2str(i),']']];
        else
            h=[h,[',h1',num2str(i)]];
        end
    end

    leg=['legend(',h,''];
    for i=1:spect
        leg=[leg,',string',num2str(i)];
    end
    leg=[leg,',''Location'',''NorthWest'');'];
    [LEGH,OBJH,OUTH,OUTM] =eval(leg);
else
    [LEGH,OBJH,OUTH,OUTM]=legend(h11,string1,'Location','NorthWest');
end

% Change the axes lines color to black
for i=1:spect
    eval(['set(ax',num2str(i),'(1),''ycolor'',[0 0 0]);']);
    eval(['set(ax',num2str(i),'(2),''ycolor'',[0 0 0]);']);
end

% Put the units on the axes
eval(['set(get(ax',num2str(spect),'(1),''Xlabel''),''String'',''Frequency (rad/s)'');']);
un1=['PSD (',units,'^2/rad)'];
eval(['set(get(ax',num2str(spect),'(1),''Ylabel''),''String'',un1);']);
un2=['Cumulative PSD (',units,'^2)'];
eval(['set(get(ax',num2str(spect),'(2),''Ylabel''),''String'',un2,''VerticalAlignment'',''top'');']);

% Change the limits on the cumulative spectra y-axis
xM=[];
xm=[];
% Get the max/min ticks
for i=1:spect
    x=eval(['get(ax',num2str(i),'(2),''YLim'');']);
    xM=max([xM,x]);
    xm=min([xm,x]);
end

% Change the number of values of the ticks
n=length(eval(['get(ax',num2str(i),'(1),''Ytick'');']))-1;
for i=1:spect
    eval(['set(ax',num2str(i),'(2),''Ylim'',[xm,xM]);']);
    %     eval(['set(ax',num2str(i),'(2),''Ylim'',[0.9999*xm,1.0001*xM]);']);
    %     eval(['set(ax',num2str(i),'(2),''YTick'',[0.9999*xm:(1.0001*xM-xm*0.9999)/n:1.0001*xM]);']);
    %    eval(['set(ax',num2str(i),'(2),''YTickLabel'',two_decimals([0.999*xm:(xM-xm)/n:1.001*xM]));']);
end

% Save the limits of the cumulative spectra axis in AX0
% AX0 is a 3D matrix:
%   1st row: x-axis, 2nd row: y-axis
%   1st thrid dimension : spectra axis, 2nd thrid dimension : cummulative
%   spectra axis
AX0(2,:,2)=[xm,xM];


% Change the limits on the cumulative spectra x-axis
xM=[];
xm=[];
for i=1:spect
    x=eval(['get(ax',num2str(i),'(2),''XLim'');']);
    % % %     xM=max([xM,x]);
    xm=min([xm,x]);
end
% % % xM = 31.6;
% % % xm = 0.0316;
xM = 2*pi*fx(end);
xM = 10^(log10(xM)*1.04);


for i=1:spect
    eval(['set(ax',num2str(i),'(1),''Xlim'',[xm,xM]);']);
    %     eval(['set(ax',num2str(i),'(1),''XTick'',10.^[log10(xm):log10(xM)]);']);
    eval(['set(ax',num2str(i),'(2),''Xlim'',[xm,xM]);']);
    %     eval(['set(ax',num2str(i),'(2),''XTick'',10.^[log10(xm):log10(xM)]);']);
end

% Save the limits of the cumulative spectra axis in AX0
AX0(1,:,1)=[xm,xM];
AX0(1,:,2)=[xm,xM];


% Change the limits on the spectra y-axis
xM=[];
xm=[];
for i=1:spect
    x=eval(['get(h1',num2str(i),',''YData'');']);
    xM=max([xM,x]);
    xm=min([xm,x]);
end


for i=1:spect
    if (xm==xM)
        if(xm==0)
            eval(['set(ax',num2str(i),'(1),''Ylim'',[0,1]);']);
        else
            eval(['set(ax',num2str(i),'(1),''Ylim'',[0.9*xm,1.1*xM]);']);
        end
    else
        eval(['set(ax',num2str(i),'(1),''Ylim'',[xm,xM]);']);
    end
end

AX0(2,:,1)=get(ax1(1),'Ylim');


% Save the axes
AX1=[];         % Spectra axis
AX2=[];         % Cummulative axis
for i=1:spect
    axx1=eval(['ax',num2str(i),'(1)']);
    axx2=eval(['ax',num2str(i),'(2)']);
    AX1=[AX1,axx1];
    AX2=[AX2,axx2];
end
% Link the axes
linkaxes(AX1,'xy');
linkaxes(AX2,'xy');

% Set the zoom functions
h = zoom(gca);
set(h,'ActionPreCallback',@myprecallback);
set(h,'ActionPostCallback',@mypostcallback)

% Save axes, etc as UserData on the figure
Data.h = h;
Data.AX1 =AX1;
Data.AX2 =AX2;
Data.yl = [];
Data.AX0 = AX0;
set(gcf,'UserData',Data)


hold off

function myprecallback(obj,evd)
UD = get(obj,'UserData');
yl = get(UD.AX2(end),'ylim');
UD.yl = yl;
set(obj,'UserData',UD);


function mypostcallback(obj,evd)
UD = get(obj,'UserData');
AX1 = UD.AX1;
AX0 = UD.AX0;
AX2 = UD.AX2;
yl = UD.yl;
h = UD.h;

for i=1:length(AX1)
    if(get(obj,'CurrentAxes')==AX1(i))
        cc=1;
    end
end

for i=1:length(AX2)
    if(get(obj,'CurrentAxes')==AX2(i))
        cc=1;
    end
end

% if cc==1
% CC0001: sort out zoom in, zoom out warnings
if exist('cc', 'var') && cc==1
% CC0001
    if strcmp(get(h,'Direction'),'in')
        ylim=get(AX1(1),'ylim');
        newXlim = get(evd.Axes,'XLim');
        newYlim = get(evd.Axes,'YLim');
        set(AX1(1),'Xlim',newXlim);

        xp=ylim(1)*10^((log10(ylim(2))-log10(ylim(1)))/(yl(2)-yl(1))*(newYlim(1)-yl(1)));
        xpp=ylim(1)*10^((log10(ylim(2))-log10(ylim(1)))/(yl(2)-yl(1))*(newYlim(2)-yl(1)));

        set(AX1(1),'Ylim',[xp,xpp]);
        UD.yl = newYlim;

    else
        set(AX1(1),'Ylim',AX0(2,:,1));
        set(AX2(1),'Ylim',AX0(2,:,2));
        set(AX1(1),'Xlim',AX0(1,:,1));
        set(AX2(1),'Xlim',AX0(1,:,2));
    end
    % for i=length(AX1)
    %     set(AX1(i),'Ytickmode','auto')
    %     set(AX2(i),'Ytickmode','auto')
    % end

    UD.AX1 = AX1;
    UD.AX2 = AX2;

    set(obj,'UserData',UD);
    cc=0;
end
