procit=1;

 %% Initial data processing
 clear all
 figstr='giffig'; 

 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);

 %% 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)];


figure(1); clf; giffig
 idx1=find(dZdT<-.3);
 idx2=find(dZdT>.3);
 idx3=find(abs(dZdT)<.2);
subplot(211); plot(dZdT,ph,'k.'); hold on; 
 plot(dZdT(idx1),ph(idx1),'ro');
 plot(dZdT(idx2),ph(idx2),'bo');
 plot(dZdT(idx3),ph(idx3),'mo');
 xlabel('depth rate (m/s)');  ylabel('ph');
 grid
subplot(212); plot(dZdT,Depth,'k.'); hold on; 
 plot(dZdT(idx1),Depth(idx1),'ro');
 plot(dZdT(idx2),Depth(idx2),'bo');
 plot(dZdT(idx3),Depth(idx3),'mo');
 xlabel('depth rate (m/s)');  ylabel('Depth');
 grid 

figure(2); clf; giffig

pHcor = ph;	% corrected pH

%%% Bin pH ranges and define pH as a function of dZ/dT
% Zrange = [0 50; 50 100;100 150;150 200;200 250;250:300];
phrange = [7.9 8.1; 7.8 7.9; 7.7 7.8; 7.6 7.7; 7.5 7.6]; nrange=length(phrange);
raterange = [-.8 -.4;-.4 -.2; -.2 0; 0 .2; .2 .4; .4 .8];
for R = 1:length(phrange);
 midx = find(ph >=phrange(R,1) & ph < phrange(R,2));
 cph = ph(midx); cdZdT=dZdT(midx); 
 for N=1:length(raterange);
  crange = raterange(N,:);
  idx=find(cdZdT > crange(1) & cdZdT <= crange(2));
  mdepthrate(N)=mean(cdZdT(idx)); mph(N)=mean(cph(idx));
 end

 eval(['subplot(' int2str(nrange) '1' int2str(R) ')']);
 plot(mdepthrate,mph,'*'); hold on
 ylabel('mean pH'); 
 title(['pH range ' num2str(phrange(R,1)) ' to ' num2str(phrange(R,2))]);
 axis([-.6 .6 min(mph)-.01 max(mph)+.01]);
 hold on;
 plot([0 0],[min(mph)-.01 max(mph)+.01],'k--');
end

%% The key problem is in the pH range 7.9 and above
R=1;
midx = find(ph >=phrange(R,1) & ph < phrange(R,2));
cph = ph(midx); cdZdT=dZdT(midx);
for N=1:length(raterange);
 crange = raterange(N,:);
 idx=find(cdZdT > crange(1) & cdZdT <= crange(2));
 mdepthrate(N)=mean(cdZdT(idx)); mph(N)=mean(cph(idx));
end

figure(3);
 plot(mdepthrate,mph,'*'); hold on
 idx = find(mdepthrate <0); crate = mdepthrate(idx); cpH = mph(idx);
 %% Anchor at lowest depth rate (min pH)
 fitrate = [repmat(max(crate),1,10) crate];
 fitmph = [repmat(min(cpH),1,10) cpH];
 %% Model II regression
 [m,b,r,sm,sb] = lsqfitma(fitrate,fitmph);
 prate = linspace(min(crate),max(crate),10); pmph = m*prate+b;
 plot(prate,pmph,'r'); 
 ylabel('mean pH'); 
 title(['pH range ' num2str(phrange(R,1)) ' to ' num2str(phrange(R,2))]);
 axis([-.6 .6 min(mph)-.01 max(mph)+.01]);
 hold on;
 plot([0 0],[min(mph)-.01 max(mph)+.01],'k--');

 %% Now you have a slope and intercept for defining average depth error
 %% from "true" pH in the range R, with true taken as that measured at
 %% the slowest rate.  This is in a mean sense, and you apply it to the
 %% data before gridding.
 %% Define points that require adjustment
 idx1 = find(cdZdT < max(crate)); 
 wrate=cdZdT(idx1); wph = cph(idx1);
 wX=max(crate)-wrate;	% positive values mean fast depth rate requiring
  			% adjustment
 %% Adjust values that have excessive depth rate
 idx=find(wX<=0); wX(idx)=NaN;
 offset = m*wX; 	% offsets are all negative
 idx = find(finite(offset)); wph(idx) = wph(idx)+offset;

 ccph = cph; ccph(idx1) = wph; clf; plot(cph,ccph,'*') 
 ylabel('Corrected pH'); xlabel('Current pH');
 axis([7.9 max(cph) 7.9 max(cph)]);
 keyboard

 %% Place back into cph
 cph(idx1) = wph;
 %% Place back in ph
 pHcor(midx) = cph;

disp('Using griddata...'); meth='linear';
 gP = griddata(tdist,Depth,ph,idist,iZ,meth);
 pidx = find(abs(dZdT)<0.3);
 gPc = griddata(tdist,Depth,pHcor,idist,iZ,meth);
 idx=find(finite(dZdT));
 gdepthrate = griddata(tdist(idx),Depth(idx),dZdT(idx),idist,iZ,meth);

 figure(3); clf; giffig
 subplot(311);
 pcolor(idist,iZ,gdepthrate); shading interp;
 colorbar; set(gca,'ydir','reverse'); axis(AX);
 xlabel('Distance along transect (km)'); ylabel('Depth (m)');
 title(['depth rate']);
 hold on; p=plot(tdist,Depth,'k.'); set(p,'markersize',4);

 subplot(312);
 pcolor(idist,iZ,gP); shading interp;
 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);
 cax=caxis;

 subplot(313);
 pcolor(idist,iZ,gPc); shading interp;
 caxis(cax);
 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(pidx),Depth(pidx),'k.'); set(p,'markersize',4);

save pHcor pHcor
