function psdp(varargin)
% ------------------------------------------------------------------
% Power Spectral Density of log pings.
%
% SYNOPSIS:
%
%   psdp(x, unit, sampling, vars, vars, ...)
%
%      Calculate and plot PSD for variables 'vars'. This argument
%      must be a column vector containing the numbers associated to the
%      variables. 
%
%      The 'x' parameter is the time variable which is used for
%      resampling of data and rescaling of frequencies.
%
%      The 'unit' parameter specifies the unit of the time variable.
%      Possible unit string are:
%
%         'yd' : year day
%
%      The 'sampling' parameter specifies the resampling interval in
%      seconds.
%
% Only data satisfying the limit tests defined in LIM are used to calculate
% the spectral density.
% ------------------------------------------------------------------
%Time:<Last updated on 98/10/06 at 21:56:50 by veins@TORSK>
%File:<c:/codas3/gfi/matlab/matlab50/scan_bin/psdp.m>

scanglob;

% Initializations:
if ~isempty(PSDP_FIG),
  delete(PSDP_FIG);
  PSDP_FIG = [];
  PSDP_SUB = [];
end
PSDP_VARS = [];

% Get input arguments:
PSDP_VARS(1) = varargin{1};
unit = varargin{2};
dx = varargin{3};
ii = 1:length(varargin{4});
j  = 2;
for i=ii,
  PSDP_VARS(j) = varargin{4}(i);
  j = j+1;
end

% Process unit argument:
switch unit
  case 'yd',
	 timetf = 'psdp_yd2s';
  otherwise,
	 error([  'Unknown unit ' unit ]);
end

% Process sampling argument:
if dx < 0,
  dx = abs(dx);
  if ~isempty(STAT),
	 dx = STAT(dx,5);
  else
	 error('No statistic array available for sampling rate');
  end
end

% Initializations:
nv = length(PSDP_VARS);
x = PSDP_VARS(1);
PSDP_DATA = zeros(nv, NPSDP);
startp = 1;
err = 0;

% Loop through the pings:
while (err == 0),
  
  % Apply psdp for a section of NPSDP pings:
  [startp, err] = sub_psdp(startp, unit, dx);
  
  % Plot current psd results:
  if err == 0,
	 psdp_plot;
	 pause(0.5);  % to allow plotting current results
  else
	 if err == -10,
		err = 0;   % Bad data, no updates
	 end
  end
  
end

% Normalize the spectrum
s = psdp_normalize;

% Replot normalized values:
psdp_plot;

% Put some additionnal information on the plot:
nextps = next_psf(PSFILE);
psdp_plot_info(dx, s, nextps);

% Save current plot:
psdp_save_plot(unit, s, nextps);
	
return;	
  
% ------------------------------------------------------------------
% Apply psdp for a subsection of NPSDP pings, starting from
% ping number sp
function [sp, err] = sub_psdp(sp, unit, dx)

scanglob;

% Get log pings:
[err, DATA] = getp(sp, NPSDP);
if err > 0,
  error([ 'Error while reading pings' ]);
end
if isempty(DATA),
  err = -1;
  return;
end
sp = sp + size(DATA,2);

% Some variables:
nv = length(PSDP_VARS);
fs = 1/dx;

% Flag log pings:
[DATA, flags] = flagdata(DATA, LIM, '-hdg', HDG);

% Resample data:
idata = psdp_resample(unit, dx);
if isempty(idata),
  err = -10;
  return;
end

% Calculate psd:
nfft = min(256, size(idata, 2));
for i=2:nv,
  [psd_data, psd_freq] = psd(idata(i,:), nfft, fs, 'mean');
  nf = length(psd_freq);
  psd_data = psd_data';
  psd_freq = psd_freq';
  % Resample to same frequencies:
  if length(find(PSDP_DATA(1,:) ~= 0)) > 0,
	 nf1 = min(find(PSDP_DATA(1,2:NPSDP) == 0));
	 psd_data = interp1(psd_freq, psd_data, PSDP_DATA(1,1:nf1));
  else
	 nf1 = nf;
	 PSDP_DATA(1,1:nf1) = psd_freq;
  end
  % Warning: if summing a `NaN' value, the result is `NaN'
  psd_data(find(isnan(psd_data))) = 0.0;
  PSDP_DATA(i,1:nf1) = PSDP_DATA(i,1:nf1) + psd_data;
end

return;

% ------------------------------------------------------------------
function data = psdp_resample(unit, dx)

scanglob;

% Some variables:
nv = length(PSDP_VARS);
np = size(DATA,2);

% Units transformations:
cx = psdp_convert(DATA(PSDP_VARS(1),:), unit);
xmin = min(cx);
xmax = max(cx);

% New time samples:
ix = xmin:dx:xmax;
ix = ix(find(~isnan(ix)));
data = [];
if isempty(ix),
  return;
end
data = ones(nv, length(ix)) * NaN;

% Resamples variables:
for i=2:nv,
  good = find(~isnan(DATA(PSDP_VARS(i),:)));
  if ~isempty(good),
	 data(i,:) = interp1(cx(good), DATA(PSDP_VARS(i),good), ix); 
  end
end
data(1,:) = ix;

return;

% ------------------------------------------------------------------
function d = psdp_convert(d, unit)

scanglob;

switch unit
  case 'yd',
	 good = find(~isnan(d));
	 if ~isempty(good),
		d(good) = d(good) * 86400.0;
	 end
end

return;

% ------------------------------------------------------------------
function psdp_init_plot

scanglob;

% Initialize figure:
PSDP_FIG = figure('NumberTitle', 'off', 'Name', 'PSDP', ...
	 'PaperType', 'a4letter');
PSDP_SUB = [];

% Calculate size of plots:
[pos, np, m] = psdp_plot_sizes;

% Create plots:
for n = 1:np,
  pos = psdp_plot_sizes(n);
  l = pos(1); b = pos(2); w = pos(3); h = pos(4);
  PSDP_SUB(n) = subplot('position', pos);
end

return;

% ------------------------------------------------------------------
function varargout = psdp_plot_sizes(varargin)

scanglob;

no_plot = 0;
if nargin > 0,
  no_plot = varargin{1};
end

% Calculate size of plots:
np = length(PSDP_VARS)-1;
m = 0.08;
h = (1-2*m)/np - 1.5*m;
b = 1 - 2*m - h;

% Number of variables per plot:
nvmax = 1;

% Calculate plot indentation:
l = (nvmax+1)*m*0.75;
w = 1 - l - 1.5*m;

% Get sizes for current subplot:
if no_plot > 0,
  b = b - (no_plot-1)*(1.5*m+h);
end

% Output arguments:
if nargout > 0,
  varargout{1} = [l, b, w, h];
end
if nargout > 1,
  varargout{2} = np;
end
if nargout > 2,
  varargout{3} = m;
end

return;

% ------------------------------------------------------------------
function psdp_plot_labels

scanglob;

figure(PSDP_FIG);

old_axes = gca;
aa = axes('position', [0 0 1 1]);
axis('off');

% Get plot sizes:
[pos, np, m] = psdp_plot_sizes;

for n=1:np,
  % Get plot sizes:
  pos = psdp_plot_sizes(n);
  l = pos(1); b = pos(2); w = pos(3); h = pos(4);  
  % Put labels:
  text(l+m, b-m, 'Frequency (Hz)', 'FontSize', 8, 'Units', 'normalized', ...
		'Interpreter', 'none');
  y = PSDP_VARS(n+1);
  text(m/2, b, HDR{y}, 'FontSize', 8, ...
		'Units', 'normalized', 'Rotation', 90, 'Interpreter', 'none');
end

axes(old_axes);

return;
  
% ------------------------------------------------------------------
function psdp_plot

scanglob;

if isempty(PSDP_FIG),
  psdp_init_plot;
end

np = length(PSDP_SUB);
nx = size(PSDP_DATA, 2);
ix = min(find(PSDP_DATA(1,2:nx) == 0))-1;

for n=1:np,
  subplot(PSDP_SUB(n));
  iy = n+1;
  plot(PSDP_DATA(1,1:ix), PSDP_DATA(iy,1:ix))
  set(PSDP_SUB(n), 'FontSize', 8);
  grid on;
end

psdp_plot_labels;

return;

% ------------------------------------------------------------------
function s = psdp_normalize

scanglob;

% Some usefull values:
nv = length(PSDP_VARS);
nf = min(find(PSDP_DATA(1,2:NPSDP) == 0));

l = PSDP_DATA(1,nf) - PSDP_DATA(1,1);
s = trapz(PSDP_DATA(1,1:nf), PSDP_DATA(2:nv, 1:nf), 2)/l;
for n=2:nv,
  PSDP_DATA(n,:) = PSDP_DATA(n,:)/s(n-1);
end

return;

% ------------------------------------------------------------------
function psdp_plot_info(dx, s, nextps)

scanglob;

% Activate plot:
figure(PSDP_FIG);

note = [ ' - BIN File: ' BINFILE ' - PS File: ' nextps ];
info = sprintf( '%s\nNPSDP = %-d - Fs = %-10.5f', ...
	 'Normalized Power Spectral Density', NPSDP, 1/dx);
info = sprintf( '%s\nN =', info);
np = length(PSDP_SUB);
for i=1:np,
  info = sprintf('%s %10.5f', info, s(i));
end
plotstuf('-title', CRUISE_ID, '-subtitle', TIME_RANGE, '-note', note, ...
	 '-info', info, '-institute', INSTITUTE);

return;

% ------------------------------------------------------------------
function psdp_save_plot(unit, s, nextps)

scanglob;

% Some usefull variables:
nv = length(PSDP_VARS);

% Get file name and save plot:
figure(PSDP_FIG);
orient tall
eval([ 'print -dps ' nextps ]);

% Add here some info to be output in ASCFILE
hdr = sprintf(' PSDP PLOT: %s\n', nextps);
hdr = sprintf('%s\n\n TIME VARIABLE: %d', hdr, PSDP_VARS(1));
hdr = sprintf('%s\n UNIT FOR TIME: %s', hdr, unit);
hdr = sprintf('%s\n PINGS NUMBER : %d', hdr, NPSDP);
hdr = sprintf('%s\n VARIABLES    : \n\n\t ..y  Normalization', hdr);
for n=2:nv,
  hdr = sprintf('%s\n\t %3d  %10.5f', hdr, PSDP_VARS(n), s(n-1));
end
logp(1, hdr);
fprintf('\n\n File is saved into \n\t%s\n\n', nextps);

return

