
function pH = load_M1_pH(buoy,dep,ctd,pHcal)

% buoy = 'm1-slave1';
% dep = '201507';
% ctd = load_M1_CTD('m1',dep);
% pHcal = [ ...
%   datenum([2015 10 23]), datenum(now+1), -0.384977, 279.15 ]; %[t1 t2 E0 Ts]

%read cumulative file
 data1 = read_M1_pH_full(buoy,dep);

%read daily file(s)
 data2 = read_M1_pH_daily(buoy,dep);

%append
 data = [data1; data2];
 clear data1 data2;

%get rid of nonsense time values
 f = find( data(:,1)<datenum([1900 1 0]) | ...
   data(:,1)>now +30 );
 data(f,:)=[];
 clear f;

%take unique time values only
 [~,i]=unique(data(:,1));
 data=data(i,:);
 clear i;

%sort by time (may not be nec)
 [~,i]=sort(data(:,1));
 data=data(i,:);
 clear i;

%assign to structure
 pH.time = data(:,1);
 pH.V1 = data(:,2)/1000; %mV to V
 clear data;

%apply pH calibration
 %interp ctd temp to pH time stamp
  tempi=interp1(ctd.time,ctd.temp,pH.time);
  if isnan(tempi(end)) && ~isnan(tempi(end-1))
    tempi(end) = tempi(end-1); %don't want to intersperse NaN's
  end
 %constants
  R=8.31451; %gas constant (J mol-1 K-1)
  F=96487; %converts between Coulombs and Faradays (1 F = 96487 C)
 %from pHcal.m
  pH.pH=NaN(size(pH.time));
  for n=1:size(pHcal,1)
    f=find(pH.time>=pHcal(n,1) & pH.time<=pHcal(n,2));
    ST=(R*(tempi(f)+273.15)*log(10)) ./ F;
    E0T=pHcal(n,3)-0.001*((tempi(f)+273.15)-pHcal(n,4));
    pH.pH(f)=(pH.V1(f)-E0T)./ST; 
  end
  clear tempi R F ST E0T n f;
