function [verticalVelocity, time_dummy] = vervelo(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 % was 3
  pMethod='spline';
end

if nargin < 5 % was 4
  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