function qo=towyo_parse(t,lat,lon,z,q,K,mincast)

% function qo=towyomap(t,lat,lon,z,q,xg,yg,K,MC)
%
% t   - time of observation
% lat - latitude of each observation
% lon - longitude of each observation
% z   - depth of observation
% q   - quantity to be mapped
%
% Optional
% K   - Kernel whose value will change sign when direction of
%       travel changes sign ([1 1 1 1 -1 -1 -1 -1])
% MC  - minimum depth range of valid towyo (5)

if nargin<6
  K=[ones(4,1);-1*ones(4,1)];
end

if nargin<7
  mincast=5;
end

ks=length(K);
kt=floor(ks/2);

% Find towyo turns

zint=cumsum(sign(diff([0;z])+eps));
d=sign(conv(zint(:),K));
d=d(kt:length(d)-kt+1);

inds=[1;find(diff(d))];
inde=[find(diff(d))+1;length(z)];

% Parse and integrate data

for in=1:length(inds)

  li=inds(in):inde(in);
  
  qo.t(in,1)=mean(t(li));
  qo.lat(in,1)=mean(lat(li));
  qo.lon(in,1)=mean(lon(li));
  qo.d=nan;
  
  qo.shal(in,1)=max(z(li));
  qo.deep(in,1)=min(z(li));
  qo.totd(in,1)=abs(min(z(li))-max(z(li)));
  qo.dirn(in,1)=sign(z(inds(in))-z(inde(in)));
  
  dz=qo.dirn(in,1)*diff(z(inds(in):inde(in)));
  dq=(q(li(1:length(li)-1))+q(li(2:length(li))))/2;
  qo.sum(in,1)=abs(sum(dz.*dq));
  qo.mean(in,1)=qo.sum(in)/-sum(dz);
  [qo.max(in,1),qo.dzmax(in,1)]=max(q(li));
  qo.dzmax(in,1)=z(li(qo.dzmax(in,1))); 

end

indg=find(qo.totd>mincast);

qo.t=qo.t(indg);
qo.lat=qo.lat(indg);
qo.lon=qo.lon(indg);
qo.d=cumsum([0;geodist(qo.lat,qo.lon)]);

qo.shal  = qo.shal(indg);
qo.deep  = qo.deep(indg);
qo.totd  = qo.totd(indg);
qo.dirn  = qo.dirn(indg);

qo.sum   = qo.sum(indg);
qo.mean  = qo.mean(indg);
qo.max   = qo.max(indg);
qo.dzmax = qo.dzmax(indg);
