function spcpltjm2(P,f,cpdlim,ciph,cic,pt)
%
%  function spcpltjm2(P,f,cpdlim,ciph,cic,[pt])
%
%  jim manning's modification of MATLAB's specplot:
%
%  1993 - minor axis labeling additions and confidence interval calculations
%  1995 - option to subplot data
%
%  P is the output of Spectrum.m
%  f=(1:n)/n*sampd/2 where sampd is sample per day.
%  cpdlim is the cycle per day limit on x-axis
%  ciph   is the confidence interval of phase
%  cic    is the confidence interval of coherence
%  pt     pause between plots in seconds (optional)
%         pt=0 provides a "keyboard" stop after each plot
%         pt=99 provides one plot with 4 subplots
%  
%
%SPECPLOT Plot the output of the SPECTRUM function.
%       SPcpltjm(P,Fs), uses P, the output of SPECTRUM, and Fs, the
%       sample frequency, to successively plot:
%
%               Pxx - X Power Spectral Density & confidence.
%               Pyy - Y Power Spectral Density & confidence.
%               abs(Txy) - Transfer Function Magnitude.
%               angle(Txy) - Transfer Function Phase.
%               Cxy - Coherence Function.
%
%       The 95% confidence intervals are displayed on the power
%       spectral density curves.
%
%       SPECPLOT(P) uses normalized frequency, Fs = 2, so that 1.0 on
%       the frequency axis is half the sample rate (the Nyquist 
%       frequency).

%       J.N. Little 7-9-86
%       Revised 4-25-88 C.R. Denham.
%       Copyright (c) 1986 by the MathWorks, Inc.

[n,m] = size(P);
if nargin < 2
   Fs = 1;
end
%f = (1:n-1)/n*Fs/2;
iday=find(f<cpdlim);
if iday==[]; return; end;

if nargin<6;
  pt=input('Enter time to pause in seconds between plots: ');
  disp(' Note: pause between plots can now be passed to spcpltjm2.m');
end;

if m == 2
   c = [P(:,1)+P(:,2) P(:,1)-P(:,2)];
  else
   c = [P(:,1)+P(:,6) P(:,1)-P(:,6)];
end
c = c .* (c > 0);
plot(f(iday),P(iday,1),'g',f(iday),c(iday,1),'g--',f(iday),c(iday,2),'g--'), ...
 ylabel(' Power Spectral Density '), ...
 xlabel('cycles/day'), ...


% generate a axis along the top in period (days)
% First, decide how frequent to post days
 vec=axis;
 sy=(vec(4)-vec(3))/10.;
 text(f(1),vec(4)-sy/2,'Period (days)')
 for i=1:ceil(length(iday)/5):floor(length(iday));
   text(f(i),vec(4)-sy,sprintf('%4.1f',1/f(i)),'horizontalalignment','center')
 end;
%
%
if m == 2
   return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%% Power Spectral Density %%%%%%%%%%%%%%%%%%%%%

% plot spectral density of both a & y

if pt==99;
  subplot(2,2,1);
end;
 
c = [P(iday,1)+P(iday,6) P(iday,1)-P(iday,6)];
c = c .* (c > 0);
plot(f(iday),P(iday,1),'g',f(iday),c(iday,1),'g--',f(iday),c(iday,2),'g--',...
  'linewidth',1), ...
 ylabel(' Power Spectral Density'), ...
 hold on;

c = [P(iday,2)+P(iday,7) P(iday,2)-P(iday,7)];
c = c .* (c > 0);
plot(f(iday),P(iday,2),'r',f(iday),c(iday,1),'r--',f(iday),c(iday,2),'r--',...
  'linewidth',3), ...

disp(' differentiate lines (red is 2nd series) & add title with gtext')
%keyboard


% generate a axis along the top in period (days)
 vec=axis;
 sy=(vec(4)-vec(3))/10.;
 text(f(1),vec(4)-sy/2,'Period (days)')
 for i=1:ceil(length(iday)/5):floor(length(iday));
   text(f(i),vec(4)-sy,sprintf('%4.1f',1/f(i)),'horizontalalignment','center')
 end;

if pt==0;
 xlabel('frequency (cycles/day)')
  keyboard;
elseif pt~=99;
 xlabel('frequency (cycles/day)')
 pause(pt)
end;
hold off;



%%%%%%%%%%%%%%%%%%%%%%%%%  TRANSFER FUNCTION %%%%%%%%%%%%%%%%%
if pt==99;
  subplot(2,2,2);
end;
plot(f(iday),abs(P(iday,4))), ...
 ylabel('Gain'), ...

% generate a axis along the top in period (days)
 vec=axis;
 sy=(vec(4)-vec(3))/10.;
 text(f(1),vec(4)-sy/2,'Period (days)')
 for i=1:ceil(length(iday)/5):floor(length(iday));
   text(f(i),vec(4)-sy,sprintf('%4.1f',1/f(i)),'horizontalalignment','center')
 end;

if pt==0;
 xlabel('frequency (cycles/day)')
  keyboard;
elseif pt~=99;
 xlabel('frequency (cycles/day)')
 pause(pt)
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PHASE  %%%%%%%%%%%%%%%%%%%%%%
if pt==99;
  subplot(2,2,3);
end;

plot(f(iday),180/pi*angle(P(iday,4)),f(iday),180/pi*angle(P(iday,4))+ciph(iday)','--',...
f(iday),180/pi*angle(P(iday,4))-ciph(iday)','--');
 ylabel('Phase (degrees)'), ...
 xlabel('frequency (cycles/day)')

% generate a axis along the top in period (days)
 vec=get(gca,'ylim');
 sy=(vec(2)-vec(1))/10.;
 text(f(1),vec(2)-sy/2,'Period (days)')
 for i=1:ceil(length(iday)/5):floor(length(iday));
   text(f(i),vec(2)-sy,sprintf('%4.1f',1/f(i)),'horizontalalignment','center')
 end;

if pt==0;
  keyboard;
elseif pt~=99;
 pause(pt)
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  COHERENCE  %%%%%%%%%%%%%%%%%%%%%
if pt==99;
  subplot(2,2,4);
end;

% according to the definition coh= abs(Pxy)^2/[(Pxx)(Pyy)]
%
pp5=P(iday,5);
plot(f(iday),pp5), ...
 hold on;
 plot([min(f(iday)) max(f(iday))],[cic cic],'r');
 ylabel('Coherence'), ...
 xlabel('frequency (cycles/day)')

% generate a axis along the top in period (days)
 vec=get(gca,'ylim');
 sy=(vec(2)-vec(1))/10.;
 text(f(1),vec(2)-sy/2,'Period (days)')
 for i=1:ceil(length(iday)/5):floor(length(iday));
   text(f(i),vec(2)-sy,sprintf('%4.1f',1/f(i)),'horizontalalignment','center')
 end;

if pt==0;
  keyboard;
elseif pt~=99;
 pause(pt)
end;

 hold off;
