function [verticalVelocity, time_dummy] = verticalVelocity(pTime, pPressure, pInterval, pMethod, pNop)
% vervelo calculates vertical velocity from pressure
% nop = number of points over which average runs
% assumption: deleted values are represented by NaNs

% 28 Oct 03, HHFurey.  I'm confused.  Arguement 4 is method, and 5 is nop, above - 
% but below doesn't see it that way.  I'll change this and see what happens.
% 29 Oct 03, HHF.  OK, here it is: this routine is fed pressure and date strings that
% are clipped to be the length of the tracked segments - we should be using the entire
% pressure record available, so that we have all vertical velocities possible.
% Also, passed out variable time_dummy because this is only record we
% have of times of vertical velocity calculations.
% 19 September 2012, hhf.  I woudl liek to change hte method to interpolate
% the data onto the time string intrinsic to teh incoming data, not create
% a new time string.  This solves problem of infinite division.

% Dhruv Balwada 9 July 2012 : added lines as time_dummy was not assigned in
% methods other than "SPLINE"

if nargin < 4
  pMethod = 'spline';
end

if nargin < 5
  pNop = 1;
end


verticalVelocity = ones(size(pPressure)) * NaN;

      
if strcmp(pMethod, 'average')
  for ix = 1:length(pPressure)
    if ~((ix - pNop < 1) || (ix + pNop > length(pPressure)))
      if ~any(isnan(pPressure([ix - pNop, ix + pNop])))
        verticalVelocity(ix) = -(pPressure(ix + pNop) - pPressure(ix - pNop)) / (pInterval * 2 * pNop);
      end
    end
  end
  time_dummy = pTime;
elseif strcmp(pMethod, 'backward')
  for ix = 1:length(pPressure)
    if ~(ix - pNop < 1)
      if ~any(isnan(pPressure([ix - pNop, ix])))
        verticalVelocity(ix) = -(pPressure(ix) - pPressure(ix - pNop)) / (pInterval * pNop);
      end
    end
  end
  time_dummy = pTime;
elseif strcmp(pMethod, 'forward')
  for ix = 1:length(pPressure)
    if ~(ix + pNop > length(pPressure))
      if ~any(isnan(pPressure([ix, ix + pNop])))
        verticalVelocity(ix) = -(pPressure(ix + pNop) - pPressure(ix)) / (pInterval * pNop);
      end
    end
  end
  time_dummy = pTime;
elseif strcmp(pMethod, 'spline')
   [time_dummy, ~, w] = artoa.data.splineWithDummy( ...
       pTime, ...
       pPressure, ...
       'NaN', ...
       pInterval / 86400 ... % [dbar per day]
   );
                                                        % interval shoudl
                                                        % be in decimal
                                                        % day, not seconds.
   w = w * 1000 / (3600 * 24) ; % conversion to mm/s
   verticalVelocity = -w(:);  % vertical velocities are positive upward
end