function wspecplot(S,varargin)
%WSPECPLOT Plot a spectral density 
%
% CALL:  wspecplot(S,plotflag,linetype) 
%
%    S       = an array of spectral density structs:
%              1D  (see dat2spec)
%              2D  (see createspec)
%   1D:
%   plotflag = 1 plots the density, S, (default)
%              2 plot 10log10(S)
%	       3 plots both the above plots 
%   2D:
%   Directional spectra:             
%   plotflag = 1 polarplot S (default)
%              2 plots spectral density and the directional 
%                spreading
%              3 mesh of S
%              4 mesh of S in polar coordinates  
%   Wavenumber spectra:
%   plotflag = 1 contourplot  S (default)
%
%   lintype : specify color and lintype, see PLOT for possibilities.
%
% NOTE: - lintype may be given anywhere after S.
%       wspecplot adds  point-and-click editing of all the text 
%       objects (title, xlabel, ylabel, etc ) of the current figure
%
% See also: dat2spec, createspec, simpson

% tested on Matlab 5.3
% History:
% revised by gl 12.07.2000
%  - added third argument 0 on line 125 to prevent warning from findrfc
% revised pab 14.06.2000
% - added the possibility that S.theta is given in degrees
% - fixed bug for plotflag 2 S.type = dir. (simpson integration)
% - cleaned up some code and put it into separate functions:
%        num2pistr, findzlevel cltext1 and fixthetalabels
% revised pab 08.05.2000
%  - fixed a bug for S an array of structs
% revised jr 10.04.2000
%  - made a comment of line 127
% revised pab 28.01.2000
%  - added plotflag 4 for dir. spec. 2D
%  - added point-and-click editing of all the text objects (title, xlabel, ylabel) of the current figure
% revised es 21.01.2000 lintype also for dir.spec.   
% revised pab 20-01.2000
%  - added the possibility that S is a array of structs
% revised pab 12-01.2000
%  - added plotflag 3 for dir. spec. i.e., meshplot
%  - made the xticks and xticklabels print better 
%    for dir. spec. when plotflag==2 (line 244-257)
% revised pab 06.01.2000
%  - added checking for empty c from contourc line 190
% revised pab 22.10.1999
%  - fixed a bug: legend(txt) -> legend(txt{:}) line 147  
% revised pab 12.10.1999
%
% by pab,es 17.09.1999


P   = varargin;

ih=ishold;
Ns=length(S);
if Ns>1
  cfig=gcf;
  for ix=1:Ns,
    if ih
      newplot
    else
      figure(cfig-1+ix)
    end
    wspecplot(S(ix),P)
  end
  return
end

% Default values: 
%~~~~~~~~~~~~~~~
plotflag = 1;
lintype  = 'b-';
Nlevel   = 7;

Np    = length(P);
if Np>0
  strix = zeros(1,Np);
  for ix=1:Np, % finding symbol strings 
    strix(ix)=ischar(P{ix});
  end
  k = find(strix);
  if any(k)
    Nk = length(k);
    if Nk>1
      warning('More than 1 strings are not allowed')
    end
    lintype = P{k(1)};
    Np = Np-Nk;
    P  = {P{find(~strix)}}; % remove strings from input
  end
  if Np>0 & ~isempty(P{1}), plotflag=P{1};end
end


ftype = freqtype(S); %options are 'f' and 'w' and 'k'
switch ftype
case {'w'},
  freq      = S.w; 
  xlbl_txt  = 'Frequency (rad/s)';  
  ylbl1_txt = 'S(w) (m^2 s / rad)';
  ylbl3_txt = 'Directional Spectrum';
  zlbl_txt  = 'S(w,\theta) (m^2 s / rad^2)';
  funit     = ' (rad/s)';  
  Sunit     = ' (m^2 s / rad)'; 
case {'f'}, 
  freq      = S.f; 
  xlbl_txt  = 'Frequency (Hz)';     
  ylbl1_txt = 'S(f) (m^2 s)';
  ylbl3_txt = 'Directional Spectrum';
  zlbl_txt  = 'S(f,\theta) (m^2 s / rad)';
  funit     = ' (Hz)';  Sunit=' (m^2 s )'; 
case {'k'},
  freq      = S.k; 
  xlbl_txt  = 'Wave number (rad/m)';
  ylbl1_txt = 'S(k) (m^3/ rad)';
  funit     = ' (rad/m)'; 
  Sunit     = ' (m^3 / rad)';
  if isfield(S,'k2')
    ylbl4_txt='Wave Number Spectrum';
  end
otherwise,  error('Frequency type unknown')
end

if isfield(S,'norm') & S.norm 
  Sunit=[];funit=[];
  ylbl1_txt = 'Normalized Spectral density';
  ylbl3_txt = 'Normalized Directional Spectrum';
  ylbl4_txt = 'Normalized Wave Number Spectrum';
  if strcmp(ftype,'k')
    xlbl_txt = 'Normalized Wave number';
  else
    xlbl_txt = 'Normalized Frequency'; 
  end
end	

ylbl2_txt = 'Power spectrum (dB)';

switch lower(S.type(end-2:end))
  case {'enc','req','k1d'} % 1D plot
    S.S  = S.S(:);
    Fn   = freq(end); % Nyquist frequency
    indm = findpeaks(S.S,4,0);    % not tested!!!
    if isempty(indm)
      %disp('indm empty')  % Commented, 00.04.10 
      [maxS,indm] = max(S.S);  %to be removed when findpeaks works
    end
    maxS = max(S.S(indm));
    if isfield(S,'CI')& ~isempty(S.CI),
      maxS  = maxS*S.CI(2);
      txtCI = [num2str(100*S.p), '% CI'];
    end
    Fp = freq(indm); %peak frequency/wave number
    
    if length(indm)==1
      txt = {['fp = ' num2str(Fp,2), funit];};
    else
      txt = {['fp1 = ' num2str(Fp(1),2), funit]};
      for j=2:length(indm)
	txt(end+1) = {['fp',num2str(j),' = ', num2str(Fp(j),2), funit]};
      end
    end
    if (plotflag==3), subplot(211),end 
    if (plotflag == 1) | (plotflag ==3),% Plot in normal scale
      plot([Fp(:) Fp(:)]',[zeros(length(indm),1) S.S(indm)]',':')
      hold on
      if isfield(S,'CI'),
	plot(freq,S.S*S.CI(1), 'r:' )
	plot(freq,S.S*S.CI(2), 'r:' )
      end
      plot(freq,S.S,lintype)  ,if ~ih, hold off,end
      if ih, a=axis; else a=zeros(1,4); end
      axis([0 max(min(Fn,10*max(Fp)),a(2)) 0 max(1.01*maxS,a(4))]) 
      title('Spectral density')
      xlabel(xlbl_txt)
      ylabel(ylbl1_txt )
    end
    
    if (plotflag==3), subplot(212),end    
    
    if (plotflag == 2) | (plotflag ==3), % Plot in logaritmic scale
      ind=find(S.S>0);
      plot([1 1]*Fp,[min(10*log10(S.S(ind)/maxS)) 10*log10(S.S(indm)/maxS)],':')
      hold on 
      if isfield(S,'CI'),
	plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' )
	plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' )
      end
      plot(freq(ind),10*log10(S.S(ind)/maxS),lintype)
      if ~ih, hold off, end
      if ih, a=axis; else a=[0 0 0 0]; end
      axis([0 max(min(Fn,10*Fp),a(2)) -20 max(1.01*10*log10(1),a(4))]) % log10(maxS)
      title('Spectral density')
      xlabel(xlbl_txt)
      ylabel(ylbl2_txt )
    end      
   
    if isfield(S,'CI'),
      legend(txt{:},txtCI,0)
    else
      legend(txt{:},0)
    end
  case {'k2d'}
    c       = contour(freq,S.k2,S.S,'b');
    z_level = findzlevel(c);
    
    % label the contour levels
    textstart_x=0.05; textstart_y=0.94;
    cltext1(z_level,textstart_x,textstart_y);
       
    xlabel(xlbl_txt)
    ylabel(xlbl_txt)
    title(ylbl4_txt)
    %return
    km=max([-freq(1) freq(end) S.k2(1) -S.k2(end)]);
    axis([-km km -km km])
    hold on
    plot([0 0],[ -km km],':')
    plot([-km km],[0 0],':')
    axis('square')
    %cltext(z_level);
    %axis('square')
    if ~ih, hold off,end
  case {'dir'}
    thmin = S.theta(1);thmax=S.theta(end);
    if plotflag==1 % polar plot
      if 0, % alternative but then z_level must be chosen beforehand
	h = polar([0 2*pi],[0 freq(end)]);
	delete(h);hold on
	[X,Y]=meshgrid(S.theta,freq);
	[X,Y]=pol2cart(X,Y);
	contour(X,Y,S.S',lintype)
      else
	if (abs(thmax-thmin)<3*pi), % angle given in radians
	  theta = S.theta;
	else
	  theta = S.theta*pi/180; % convert to radians
	end
	c = contours(theta,freq,S.S');%,Nlevel); % calculate levels
	if isempty(c)
	  c = contours(theta,freq,S.S);%,Nlevel); % calculate levels
	end
	[z_level c] = findzlevel(c); % find contour levels
	polar(c(1,:),c(2,:),lintype);
      end
      title(ylbl3_txt)
      % label the contour levels
      textstart_x = -0.1; textstart_y=1.00;
      cltext1(z_level,textstart_x,textstart_y) % a local variant of cltext
      %cltext(z_level)
    elseif plotflag==2 
      subplot(211)
      Sf = spec2spec(S,'freq'); % frequency spectrum
      wspecplot(Sf,1,lintype)

      subplot(212)
      Dtheta  = simpson(freq,S.S,2); %Directional spreading
      [y,ind] = max(Dtheta);
      Wdir    = S.theta(ind); % main wave direction
      txtwdir = ['\thetap=' num2pistr(Wdir,3)]; % convert to text string
           
      plot([1 1]*S.theta(ind),[0 Dtheta(ind)],':'), hold on
      lh=legend(txtwdir,0);
      plot(S.theta,Dtheta,lintype)
      
      if ~ih, hold off, end

      fixthetalabels(thmin,thmax,'x',2)  % fix xticklabel and xlabel for theta
      ylabel('D(\theta)')
      title('Spreading function')
      %legend(lh) % refresh current legend
    elseif plotflag==3 % mesh
      mesh(freq,S.theta,S.S)
      xlabel(xlbl_txt);
      fixthetalabels(thmin,thmax,'y',3) % fix yticklabel and ylabel for theta
      zlabel(zlbl_txt)
      title(ylbl3_txt)
    elseif plotflag==4 % mesh
      %h=polar([0 2*pi],[0 freq(end)]);
      %delete(h);hold on
      [X,Y]=meshgrid(S.theta,freq);
      [X,Y]=pol2cart(X,Y);
      mesh(X,Y,S.S')
      % display the unit circle beneath the surface
      hold on, mesh(X,Y,zeros(size(S.S'))),hold off
      zlabel(zlbl_txt)
      title(ylbl3_txt)
      set(gca,'xticklabel','','yticklabel','')
      lighting phong
      %lighting gouraud
      %light
    elseif plotflag==5
      if (abs(thmax-thmin)<3*pi), % angle given in radians
	theta = S.theta*180/pi; % convert to degrees
	thmin = thmin*180/pi;
	thmax = thmax*180/pi; 
      else
	theta = S.theta
      end
      c = contour(freq,theta,S.S,Nlevel); % calculate levels
      if isempty(c)
	c = contour(freq,theta,S.S');%,Nlevel); % calculate levels
      end
      
      fixthetalabels(thmin,thmax,'y',2) % fix yticklabel and ylabel for theta
      title(ylbl3_txt)
      xlabel(xlbl_txt);
      if 0,
	[z_level] = findzlevel(c); % find contour levels
	% label the contour levels
	textstart_x = 0.06; textstart_y=0.94;
	cltext1(z_level,textstart_x,textstart_y) % a local variant of cltext
      else
	colormap('jet')
	hcb = colorbar;
	grid on
      end
    end
	
  otherwise, error('unknown spectral type')
end

%  The following two commands install point-and-click editing of
%   all the text objects (title, xlabel, ylabel) of the current figure:

set(findall(gcf,'type','text'),'buttondownfcn','edtext')
set(gcf,'windowbuttondownfcn','edtext(''hide'')')

return

function  xtxt = num2pistr(x,N)
% NUM2PISTR convert a scalar x to a text string in fractions of pi
%           if the numerator is less than 10 and not equal 0 
%           and if the denominator is less than 10.
%
%  CALL xtxt = num2pistr(x,n)
% 
% xtxt = a text string in fractions of pi
%  x   = a scalar
%  n   = maximum digits of precision. (default 3)


if nargin<2|isempty(N)
  N=3;
end
[num den]=rat(x/pi);
if (den<10)& (num<10)&(num~=0),
  if abs(den)==1, dtxt=[]; else, dtxt=['/' num2str(den)]; end % denominator
  if abs(num)==1, % numerator
    if num==-1, ntxt='-'; else, ntxt=[];end 
  else 
    ntxt=num2str(num); 
  end
  xtxt= [ntxt '\pi' dtxt];
else
  xtxt = num2str(x,N);
end
return  
  
function fixthetalabels(thmin,thmax,xy,dim)
% FIXTHETALABELS pretty prints the ticklabels and x or y labels for theta  
%  
% CALL fixthetalabels(thmin,thmax,xy )
%
%  thmin, thmax = minimum and maximum value for theta (wave direction)
%  xy           = 'x' if theta is plotted on the x-axis 
%                 'y' if theta is plotted on the y-axis 
%  dim          = specifies the dimension of the plot (ie number of axes shown 2 or 3)
%  If abs(thmax-thmin)<3*pi it is assumed that theta is given in radians 
%  otherwise degrees

ind = [('x' == xy)  ('y' == xy) ];
yx = 'yx';
yx = yx(ind);
if nargin<4|isempty(dim),
  dim=2;
end

if abs(thmax-thmin)<3*pi, %  radians given. Want xticks given as fractions  of pi
  set(gca,[xy 'tick'],pi*(thmin/pi:0.25:thmax/pi))
  set(gca,[xy 'ticklabel'],[])
  x=get(gca,[xy 'tick']);
  y=get(gca,[yx 'tick']);
  y=y(1)-(y(2)-y(1))/3;
  xlim=get(gca,[xy 'lim']);xlim=xlim(2);
  
  if xy=='x'
     for j=1:length(x)
      xtxt = num2pistr(x(j));
      figtext(x(j),y,xtxt,'data','data','center');
    end
    ax = [thmin thmax 0 inf];
    
    if dim<3,
      figtext(mean(x),2*y,'Wave directions (rad)','data','data','center')
    else
      ax = [ax  0 inf];
      xlabel('Wave directions (rad)')
    end
  else
    for j=1:length(x)
      xtxt = num2pistr(x(j));
      figtext(y,x(j),xtxt,'data','data','center');
    end
    ax = [0 inf thmin thmax];
   
    if dim<3,
      figtext(2*y,mean(x),'Wave directions (rad)','data','data','center')
    else
      ax = [ax 0 inf];
      ylabel('Wave directions (rad)')
    end
  end
  %xtxt = num2pistr(x(j));
  %for j=2:length(x)
  %  xtxt = strvcat(xtxt,num2pistr(x(j)));
  %end
  %set(gca,[xy 'ticklabel'],xtxt)
else % Degrees given
  set(gca,[xy 'tick'],thmin:45:thmax)
  if xy=='x'
    ax=[thmin thmax 0 inf];
    if dim>=3,      ax=[ax 0 inf];    end
    xlabel('Wave directions (deg)')
  else
    ax=[0 inf thmin thmax ];
    if dim>=3,      ax=[ax 0 inf];    end
    ylabel('Wave directions (deg)')
  end
end
axis(ax)
return

function [z_level, c] = findzlevel(c)
% FINDZLEVEL find the contour levels calculated by the contourc function
%
% CALL [z_level, c1] = findzlevel(c)
%
%  example:
%
%  c = contour(x,y,Z)
%  [z_level, c1] = findzlevel(c)
%  polar(c1(:,1),c1(:,2))
%
limit = size(c,2);
ix = 1; iy=1;
while(iy < limit)
  z_level(ix) = c(1,iy);
  npoints = c(2,iy);
  c(:,iy)=NaN;
  nexti = iy+npoints+1;
  iy = nexti;
  ix=ix+1;
end
% sort and remove duplicate levels
z_level = sort(z_level);
z_level = z_level(diff([0 z_level])~=0);

return

function cltext1(z_level,textstart_x,textstart_y)
% CLTEXT1 Places contour level text in the current window
%
%  [h,ax] = cltext1(levels,x,y)
%
%         x       = the starting x-coordinate of the text.
%         y       = the starting y-coordinate of the text.
%

%  textstart_x = -0.07; textstart_y=1.00;

if (nargin<4)|isempty(x_unit),
  x_unit='data';
end


delta_y     = 1/33;
%dir:
%h = figtext(textstart_x-0.07,textstart_y,' Level curves at:','norm');

% k2d:
h = figtext(textstart_x-0.05,textstart_y,' Level curves at:','norm');
set(h,'FontWeight','Bold')

%z_level = sort(z_level);
%z_level = z_level(diff([0 z_level])~=0)

textstart_y = textstart_y-delta_y;
for ix=1:length(z_level)
  textstart_y = textstart_y-delta_y;
  figtext(textstart_x,textstart_y,num2str(z_level(ix),4),'norm');
end
return