procit=1; 

if procit
%% Initial data processing
 clear all; 
 output4predict=1; filtit=1; savit=1; plotit=1;
 figstr='figure(1); clf; paper'; 
 jgbgrid=0; delaunay=1;

missions = ['2000.056.04/'; '2000.056.05/'; '2000.056.06/'; '2000.056.07/'; '2000.056.09/'];
 
 [N x] = size(missions);

 for M = 1:N;
  eval(['load ' missions(M,:) 'ctd']);
  if M==1;
   o2=O2; s=S; t=T; ph=pH; Depth=depth; time=epochsec; 
   lat=startlat; lon=startlon; 
  else
   o2=[o2;O2]; s=[s;S]; t=[t;T]; ph=[ph;pH]; Depth=[Depth;depth]; 
   time=[time;epochsec];
   lat=[lat;startlat]; lon=[lon;startlon]; 
  end
 idx(M) = length(T);	% idx contains endpoint of segment
 end

 %% Salinity has one glitch
 bidx = find(s>35); disp([int2str(length(idx)) ' salinity values > 35']);
 s(bidx)=NaN;

 %% Calculate interpolated lat and lon
 %% Add last lat/lon based on next to last 
 npt = length(lat); 
  alat = lat(npt) + (lat(npt) - lat(npt-1)); 
  alon = lon(npt) + (lon(npt) - lon(npt-1)); 
  lat=[lat;alat]; lon=[lon;alon];
 
 %% Interpolate lat/lon
 for M = 1:N;
  ilat = [linspace(lat(M),lat(M+1),idx(M))]';
  ilon = [linspace(lon(M),lon(M+1),idx(M))]';
  if M==1;
   ILAT = ilat; ILON=ilon;
  else
   ILAT = [ILAT;ilat]; ILON=[ILON;ilon];
  end
 end
 
 %% Estimate transect distance
 addpath '/u/ryjo/matlab/CSIRO'
 [dist phas]=sw_dist(ILON([1 length(ILAT)]),ILON([1 length(ILON)]),'km');
 % dist = dist*1000;	% m
 tdist = linspace(0,dist,length(ILON));
 idist = linspace(0,dist,300);
 iZ = [10:1:max(Depth)]';
 AX = [0 dist min(iZ) max(iZ)];

 keyboard
 % disp(['Gridding full resolution pH...']); meth='linear';
 % gP = griddata(tdist,Depth,ph,idist,iZ,meth);

%% Try combining up and down cast into one with a mean distance 
%% Use Depth < 10 as cutoff
Zcutoff=10;
% figure(1); clf; bigfig;
% plot(time,Depth,'.'); set(gca,'ydir','reverse');
% [x y]=ginput(18); timebins = x([[1:2:17]' [2:2:18]']); 
% save timebins timebins x
load timebins

eval(figstr);

if yoyoV
 fidx = find(Depth>=Zcutoff); fidx=fidx(1);
 didx = find(Depth<Zcutoff); ddidx=[NaN;diff(didx)];
 idx = find(ddidx>50); ends = didx(idx);
 starts = [fidx;ends([1:length(ends)-1])+1];
 plot(tdist,Depth,'k.'); set(gca,'ydir','reverse'); hold on
 plot(tdist(ends),Depth(ends),'ro');
 plot(tdist(starts),Depth(starts),'go');
 pause
end

if yoyoUV
  plot(time,Depth,'k.'); set(gca,'ydir','reverse'); hold on
  plot(timebins(:,1),repmat(0,length(timebins)),'go'); 
  plot(timebins(:,2),repmat(0,length(timebins)),'rx'); 
 pause
end

%% Ensure that no shallow data enters means
idx = find(Depth<10 | Depth>250); Depth(idx)=NaN;

%% Rebin to 2 m, combined down/up
RZ = [11:2:max(Depth)]';

if yoyoV
 Rt = repmat(NaN,length(RZ),length(starts)); Rs=Rt; Ro2=Rt; Rph=Rt;
 for N=1:length(starts);
 Rdist(N) = mean(tdist(starts(N):ends(N)));
 cs=s(starts(N):ends(N)); ct=t(starts(N):ends(N));
 cph=ph(starts(N):ends(N)); co2=o2(starts(N):ends(N));
 cZ = Depth(starts(N):ends(N)); cmxZ=max(cZ);
 for k=1:length(RZ);
  idx=find(cZ >= RZ(k)-1 & cZ < RZ(k)+1);
  if ~isempty(idx);
   cdat=ct(idx); cdat(isnan(cdat))=[];
   if ~isempty(cdat); Rt(k,N)=mean(cdat); end
   cdat=cs(idx); cdat(isnan(cdat))=[];
   if ~isempty(cdat); Rs(k,N)=mean(cdat); end
   cdat=co2(idx); cdat(isnan(cdat))=[];
   if ~isempty(cdat); Ro2(k,N)=mean(cdat); end
   cdat=cph(idx); cdat(isnan(cdat))=[];
   if ~isempty(cdat); Rph(k,N)=mean(cdat); end
  end
 end
end
end

if yoyoUV
 Rt = repmat(NaN,length(RZ),length(timebins)); Rs=Rt; Ro2=Rt; Rph=Rt;
 for N=1:length(timebins);
 a=find(time>=timebins(N,1) & time<=timebins(N,2));
 starts(N)=min(a); ends(N)=max(a);
 Rdist(N) = mean(tdist(starts(N):ends(N)));
 cs=s(starts(N):ends(N)); ct=t(starts(N):ends(N));
 cph=ph(starts(N):ends(N)); co2=o2(starts(N):ends(N));
 cZ = Depth(starts(N):ends(N)); cmxZ=max(cZ);
 for k=1:length(RZ);
  idx=find(cZ >= RZ(k)-1 & cZ < RZ(k)+1);
  if ~isempty(idx);
   cdat=ct(idx); cdat(isnan(cdat))=[]; 
   if ~isempty(cdat); Rt(k,N)=mean(cdat); end
   cdat=cs(idx); cdat(isnan(cdat))=[]; 
   if ~isempty(cdat); Rs(k,N)=mean(cdat); end
   cdat=co2(idx); cdat(isnan(cdat))=[]; 
   if ~isempty(cdat); Ro2(k,N)=mean(cdat); end
   cdat=cph(idx); cdat(isnan(cdat))=[]; 
   if ~isempty(cdat); Rph(k,N)=mean(cdat); end
  end
 end
end
end

if output4predict
 RZg = repmat(RZ,1,length(timebins)); 
 t=Rt(:); s=Rs(:); o2=Ro2(:); ph=Rph(:); Depth =RZg(:);
 idx = find(isnan(t)|isnan(s)|isnan(o2)|isnan(ph));
 t(idx)=[]; s(idx)=[]; o2(idx)=[]; ph(idx)=[]; Depth(idx)=[];
 save UD4predict t s o2 ph Depth
 return
end

% figure(1);
% set(gcf,'units','pixels','position',[20 100 1200 900],'color','w'); 
% subpl t(211);
% pcolor(idist,iZ,gP); shading interp; 
% set(gca,'ydir','reverse');
%  colorbar; hold on
% contour(idist,iZ,gP,[7.5:.1:8],'k');
% [cc,h]=contour(idist,iZ,gP,[7.9,7.98],'k'); set(h(:),'linewidth',2);
% plot(tdist,Depth,'k.'); plot(tdist(ends),Depth(ends),'ro');
% plot(tdist(starts),Depth(starts),'go');
% for N=1:length(starts); plot(repmat(Rdist(N),size(RZ)),RZ,'k.'); end
% title('Gridded full resolution');
% AX=axis;

% subplot(212);
% pcolor(Rdist,RZ,Rph); shading interp; set(gca,'ydir','reverse'); colorbar
% hold on;
% contour(Rdist,RZ,Rph,[7.5:.1:8],'k');
% [cc,h]=contour(Rdist,RZ,Rph,[7.9,7.98],'k'); set(h(:),'linewidth',2);
% for N=1:length(starts); plot(repmat(Rdist(N),size(RZ)),RZ,'k.'); end
% plot(tdist,Depth,'k.'); plot(tdist(ends),Depth(ends),'ro');
% plot(tdist(starts),Depth(starts),'go');
% title('Gridded averaged up and downcast');
% axis(AX);

%% Use griddata to grid to higher resolution
Rdist = Rdist-min(Rdist);
idist = linspace(min(Rdist),max(Rdist),150);
[gdist gZ]=meshgrid(Rdist,RZ);
disp('Using griddata...'); meth='linear';
%% T
idx = find(finite(Rt)); 
gT = griddata(gdist(idx),gZ(idx),Rt(idx),idist,RZ,meth);
%% S
idx = find(finite(Rs)); 
gS = griddata(gdist(idx),gZ(idx),Rs(idx),idist,RZ,meth);
%% O2 
idx = find(finite(Ro2)); 
gO2 = griddata(gdist(idx),gZ(idx),Ro2(idx),idist,RZ,meth);
%% pH 
idx = find(finite(Rph)); 
gpH = griddata(gdist(idx),gZ(idx),Rph(idx),idist,RZ,meth);

end	%procit

if filtit
 disp(['Filtering']);
 gS=medfiltJR(gS,5); gO2=medfiltJR(gO2,5); 
 gpH=medfiltJR(gpH,5); gT=medfiltJR(gT,5);
end

if savit
 save avgUD_GF Rdist RZ idist gT gS gpH gO2 
end

keyboard

if plotit
eval(figstr);

subplot(411);  pcolor(idist,RZ,gT); shading interp; hold on;
hold on; contour(idist,RZ,gT,[7.2:1:13.2],'k'); contour(idist,RZ,gT,[7.5:13.5],'k--');
axis([min(Rdist) max(Rdist) min(RZ) max(RZ)]);
t=text(1,250,['Temperature (' setstr(176) 'C)']); 
set(t,'fontweight','bold','fontsize',12); % caxis([8.9 12.7]); 
set(gca,'fontsize',12);
colorbar; set(gca,'ydir','reverse'); ylabel('Depth (m)'); 
% xlabel('Distance along transect (km)'); 
% title(['Temperature, observed [min max] = [' num2str([min(t) max(t)]) ']']);
% hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);

subplot(412);  pcolor(idist,RZ,gS); shading interp; hold on;
hold on; contour(idist,RZ,gS,[32:.05:35],'k'); 
axis([min(Rdist) max(Rdist) min(RZ) max(RZ)]);
t=text(1,250,['Salinity (psu)']); 
set(t,'fontweight','bold','fontsize',12); % caxis([8.9 12.7]); 
set(gca,'fontsize',12);
colorbar; set(gca,'ydir','reverse'); ylabel('Depth (m)'); 

subplot(413);  pcolor(idist,RZ,gO2); shading interp; hold on;
hold on; contour(idist,RZ,gO2,[1:.3:6],'k'); 
axis([min(Rdist) max(Rdist) min(RZ) max(RZ)]);
t=text(1,250,['Oxygen (ml/l)']); 
set(t,'fontweight','bold','fontsize',12); % caxis([8.9 12.7]); 
set(gca,'fontsize',12);
colorbar; set(gca,'ydir','reverse'); ylabel('Depth (m)'); 

subplot(414);  pcolor(idist,RZ,gpH); shading interp; hold on;
hold on; contour(idist,RZ,gpH,[7:.1:8.2],'k'); 
axis([min(Rdist) max(Rdist) min(RZ) max(RZ)]);
t=text(1,250,['pH']); 
set(t,'fontweight','bold','fontsize',12); % caxis([8.9 12.7]); 
set(gca,'fontsize',12);
colorbar; set(gca,'ydir','reverse'); ylabel('Depth (m)'); 

frame
t=text(.5,.95,titl);
set(t,'fontweight','bold','fontsize',15,'horizontalalignment','center');
end
