procit=1;
figstr='figure(1); clf; paper'; % homefig(1)'; 

if procit
%% Initial data processing
 clear all
 figstr='figure(1); clf; paper'; % homefig(1)'; 
 % figstr='figure(1); clf; paper; figure(2); 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,:) 'auv']);
  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;

 % Depth is positive
 dZdt = [NaN; diff(Depth)./diff(time)]; 
 yoyo_discriminator = 0.1;	% m/s
 going_up=find(dZdt<-1*yoyo_discriminator); 
 going_down=find(dZdt>yoyo_discriminator); 

 phup=ph(going_up); phdown=ph(going_down);
 o2up=o2(going_up); o2down=o2(going_down);
 tup=t(going_up); tdown=t(going_down);
 sup=s(going_up); sdown=s(going_down);
 Zup = Depth(going_up); Zdown=Depth(going_down);
 Zbin = [15:10:max(Depth)-5]; buf=5;
 for N = 1:length(Zbin);
  uidx=find(Zup>=Zbin(N)-buf & Zup<Zbin(N)+buf);
  didx=find(Zdown>=Zbin(N)-buf & Zdown<Zbin(N)+buf);
  pHu(N)=mean(phup(uidx)); pHd(N)=mean(phdown(didx));
  o2u(N)=mean(o2up(uidx)); o2d(N)=mean(o2down(didx));
  tu(N)=mean(tup(uidx)); td(N)=mean(tdown(didx));
  su(N)=mean(sup(uidx)); sd(N)=mean(sdown(didx));
 end
end

eval(figstr);
YL = [11 250]; figure(1); P=[.2 .2 .6 .6];

%% Oxygen
axes('position',P);
dat=o2u-o2d; p=plot(dat,Zbin,'k--'); set(p,'linewidth',2);
set(gca,'Ydir','reverse','fontsize',12); 
set(gca,'Xlim',[-1*max(abs(dat)) max(abs(dat))],'Ylim',YL);
hold on; plot([0 0],YL,'k:');
set(gca,'color','none','box','off','YTick',[],'XTick',[]);
AX=axis; axes('position',P+[0 -.07 0 0]);
axis(AX); set(gca,'YTick',[],'Ycolor','w','color','none','tickdir','out');
set(gca,'fontsize',12); xlabel('O_2 (ml/l)');
ylabel('Depth');

%% pH 
axes('position',P);
dat=pHu-pHd; p=plot(dat,Zbin,'k'); set(p,'linewidth',2);
set(gca,'Ydir','reverse','fontsize',12); 
set(gca,'Xlim',[-1*max(abs(dat)) max(abs(dat))],'Ylim',YL);
set(gca,'Tickdir','out','box','off','color','none')
xlabel('pH');

%% Salinity
axes('position',P);
dat=su-sd; p=plot(dat,Zbin,'k--'); 
set(gca,'Ydir','reverse','fontsize',12); 
set(gca,'Xlim',[-1*max(abs(dat)) max(abs(dat))],'Ylim',YL);
set(gca,'color','none','Xaxislocation','top','box','off','XTick',[],'YTick',[]);
AX=axis; axes('position',P+[0 .07 0 0]);
axis(AX); 
set(gca,'Xaxislocation','top','fontsize',12); xlabel('Salinity (psu)');
set(gca,'YTick',[],'Ycolor','w','color','none','tickdir','out');

axes('position',P);
dat=tu-td; p=plot(dat,Zbin,'k'); 
set(gca,'Ydir','reverse','fontsize',12); 
set(gca,'Xlim',[-1*max(abs(dat)) max(abs(dat))],'Ylim',YL);
set(gca,'color','none','Xaxislocation','top');
set(gca,'Tickdir','out','box','off')
xlabel(['Temperature (' setstr(176) 'C)']);

axes('position',P); axis([0 1 YL]);
set(gca,'Ydir','reverse');
t=text(.2,72,'O_2'); set(t,'fontsize',12);
t=text(.8,50,'pH'); set(t,'fontsize',12);
t=text(.9,103,'S'); set(t,'fontsize',12);
t=text(.1,228,'T'); set(t,'fontsize',12);
set(gca,'box','on','color','none','XTick',[],'YTick',[]);

return
figure(1);
subplot(411);
plot(pHd,Zbin,'b',pHu,Zbin,'r');
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
title('Binned pH'); legend('downcast','upcast');
subplot(412);
plot(o2d,Zbin,'b',o2u,Zbin,'r');
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
title('Binned o2');
subplot(413);
plot(td,Zbin,'b',tu,Zbin,'r');
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
title('Binned T'); 
subplot(414);
plot(sd,Zbin,'b',su,Zbin,'r');
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
title('Binned S'); 

figure(3); thx_fig;
dat=(pHd-pHu); p=plot(dat,Zbin,'k'); set(p,'linewidth',2);
set(gca,'Ydir','reverse','fontsize',12);
hold on;
return
%% 
set(gca,'Ydir','reverse','Xaxislocation','bottom','box','off','Yaxislocation','right','fontsize',12,'color','none');
% set(gca,'Xlim',[-1*max(abs(dat)) max(abs(dat))]);
% set(gca,'Ylim',[0 max(Zbin)]); plot([0 0],[0 max(Zbin)],'k:');
% xlabel('Average O_2 difference, downcast-upcast (ml/l)');

dat=o2d-o2u; plot(dat,Zbin,'k');
set(gca,'Ydir','reverse','Xaxislocation','bottom','box','off','Yaxislocation','right','fontsize',12,'color','none');
hold on;
set(gca,'Xlim',[-1*max(abs(dat)) max(abs(dat))]);
set(gca,'Ylim',[0 max(Zbin)]); plot([0 0],[0 max(Zbin)],'k:');
xlabel('Average O_2 difference, downcast-upcast (ml/l)');


return

subplot(411);
p=plot(t(going_down),Depth(going_down),'b.',t(going_up),Depth(going_up),'r.');
% set(p,'markersize',4);
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
set(gca,'Xlim',[min(t) max(t)]);
title('Temperature'); legend('downcast','upcast');
ylabel('Depth (m)');

subplot(412);
p=plot(o2(going_down),Depth(going_down),'b.',o2(going_up),Depth(going_up),'r.');
% set(p,'markersize',4);
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
set(gca,'Xlim',[min(o2) max(o2)]);
title('Oxygen'); % legend('downcast','upcast');
ylabel('Depth (m)');

subplot(413);
p=plot(ph(going_down),Depth(going_down),'b.',ph(going_up),Depth(going_up),'r.');
% set(p,'markersize',4);
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
set(gca,'Xlim',[min(ph) max(ph)]);
title('pH'); % legend('downcast','upcast');
ylabel('Depth (m)');

s(bidx)=[]; sDepth=Depth; sDepth(bidx)=[]; 
dZ = [NaN; diff(sDepth)]; going_down=find(dZ>0); going_up=find(dZ<0);
subplot(414);
p=plot(s(going_down),Depth(going_down),'b.',s(going_up),Depth(going_up),'r.');
% set(p,'markersize',4);
set(gca,'Ydir','reverse','Xaxislocation','top','box','on','fontsize',12);
set(gca,'Xlim',[33 max(s)]);
title('Salinity'); % legend('downcast','upcast');
ylabel('Depth (m)');

return
 %% 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,1000);
 iZ = [10:1:max(Depth)]';
 AX = [0 dist min(iZ) max(iZ)];
end

if jgbgrid
 %% Use JGB's potgrid to grid data
 %% Size of output grid
 nx = 100; ny = 50;
 %% Weighting for gridding;
 wi = 1; wj = round(dist*1000/max(Depth));
 % wi = 1; wj = 100;
 disp(['Weighting for Distance : ' int2str(wj) ' km']);
 
 figure(1); clf; giffig
 %% Temperature
 [xg,yg,zg] = potgridJR(tdist',Depth,t,[nx ny],[wi wj]);
 % idx = find(zg < min(t) | zg > max(t)); zg(idx)=NaN*idx;
 subplot(411); pcolor(xg,yg,zg); shading interp; colorbar; set(gca,'Ydir','rev');
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);
 title(['Temperature, observed [min max] = [' num2str([min(t) max(t)]) ']']);
 ylabel('Depth (m)');
 
 %% Salinity 
 s(bidx)=[]; sDepth=Depth; sDepth(bidx)=[]; sdist=tdist; sdist(bidx)=[];
 [xg,yg,zg] = potgridJR(sdist',sDepth,s,[nx ny],[wi wj]);
 % idx = find(zg < min(s) | zg > max(s)); zg(idx)=NaN*idx;
 subplot(412); pcolor(xg,yg,zg); shading interp; colorbar; set(gca,'Ydir','rev');
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);
 title(['Salinity, observed [min max] = [' num2str([min(s) max(s)]) ']']);
 ylabel('Depth (m)');
 
 %% O2 
 [xg,yg,zg] = potgridJR(tdist',Depth,o2,[nx ny],[wi wj]);
 % idx = find(zg < min(o2) | zg > max(o2)); zg(idx)=NaN*idx;
 subplot(413); pcolor(xg,yg,zg); shading interp; colorbar; set(gca,'Ydir','rev');
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);
 title(['O_2, observed [min max] = [' num2str([min(o2) max(o2)]) ']']);
 ylabel('Depth (m)');
 
 %% pH 
 [xg,yg,zg] = potgridJR(tdist',Depth,ph,[nx ny],[wi wj]);
 % idx = find(zg < min(ph) | zg > max(ph)); zg(idx)=NaN*idx;
 subplot(414); pcolor(xg,yg,zg); shading interp; colorbar; set(gca,'Ydir','rev');
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);
 title(['pH, observed [min max] = [' num2str([min(ph) max(ph)]) ']']);
 xlabel('Distance along transect (km)');
 ylabel('Depth (m)');

 frame
 txt=text(.5,.97,'Gridded AUV observations, using potgrid.m');
 set(txt,'fontweight','bold','horizontalalignment','center');
 screen2jpeg('potgrid.jpg');
end

if delaunay
 disp('Using griddata...'); meth='linear';
 gT = griddata(tdist,Depth,t,idist,iZ,meth);
  s(bidx)=[]; sDepth=Depth; sDepth(bidx)=[]; sdist=tdist; sdist(bidx)=[];
 gS = griddata(sdist,sDepth,s,idist,iZ,meth);
 gO = griddata(tdist,Depth,o2,idist,iZ,meth);
 gP = griddata(tdist,Depth,ph,idist,iZ,meth);

 figure(1); clf; giffig
 subplot(411); 
 pcolor(idist,iZ,gT); shading interp; 
 % caxis([8.9 12.7]); 
 colorbar; set(gca,'ydir','reverse'); axis(AX);
 % xlabel('Distance along transect (km)'); 
 ylabel('Depth (m)'); 
 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,iZ,gS); shading interp; 
 caxis([33.1 34.1]); 
 colorbar; set(gca,'ydir','reverse'); axis(AX);
 % xlabel('Distance along transect (km)'); 
 ylabel('Depth (m)'); 
 title(['Salinity, observed [min max] = [' num2str([min(s) max(s)]) ']']);
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);
 
 subplot(413); 
 pcolor(idist,iZ,gO); shading interp; 
 % caxis([8.9 12.7]); 
 colorbar; set(gca,'ydir','reverse'); axis(AX);
 % xlabel('Distance along transect (km)'); 
 title(['O_2, observed [min max] = [' num2str([min(o2) max(o2)]) ']']);
 ylabel('Depth (m)'); 
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);
 
 subplot(414); 
 pcolor(idist,iZ,gP); shading interp; 
 % caxis([8.9 12.7]); 
 colorbar; set(gca,'ydir','reverse'); axis(AX);
 xlabel('Distance along transect (km)'); ylabel('Depth (m)'); 
 title(['pH, observed [min max] = [' num2str([min(ph) max(ph)]) ']']);
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);
 frame
 txt=text(.5,.97,'Gridded AUV observations, using Delaunay triagulation');
 set(txt,'fontweight','bold','horizontalalignment','center');
 screen2jpeg('delaunay.jpg');
end


