
function pHcal()

% Description of inputs:

% 1. Coeff's obtained during TC (temperature calibration)
%  (TC's are typically done following a 'SC' (sensor configuration), i.e.
%  if an electrode failed, etc. These are properties of the board itself.)
% 
%    VbatSF: battery voltage Scale Factor ('R')
%    slope: 
%    intercept: 

% 2. Coeff's obtained during PC (pH calibration)
%  (PC's are done every few months.)
% 
%    Ts: water temperature during calibration [K]
%    E0: a function of Vout, Ts, and pHstd:
%      E0 = Vout - (pHstd*R*Ts*log(10)/F), where:
%      Vout: voltage output (V1) of electrode at calibration -
%            a single m'ment or a few m'ments avg'd
%      pHstd: the pH of the standard (e.g., CRM Batch 133, 136, etc.,
%             using Total scale)

%constants
 R=8.31451; %gas constant (J mol-1 K-1)
 F=96487; %to convert between Coulombs and Faradays (1 F = 96487 C)

%user-select raw pH data file from logger
 [myfile,mydir,~] = uigetfile('*.csv','Select pH data file');
 disp(['pH Data File: ' mydir myfile]);

%open file
 fid=fopen([mydir myfile]);

%go thru header, get instrument S/N
 fline=[]; 
 ct=0;
 while isempty(strfind(fline,'"[End Display]"'))
   ct=ct+1;
   fline=fgetl(fid);
   header{ct}=fline;
  %S/N
   if ~isempty(strfind(fline,'"SerialNumber"'))
     f=strfind(fline,',');
     SN=fline(f(end)+2:end-1);
   end
 end
 clear f ct;
 %get last 3 (empty) header lines after "[End Display]"
  for n=1:3
    fgetl(fid);
  end
  clear n ans;
 disp(['Instrument Serial Number: ' SN]);

%user-define cal coeff's
 VbatSF = input(['Enter Temperature Calibration ' ...
   'Voltage Scale Factor (''R''): ']);
 slope = input('Enter Temperature Calibration Slope: ');
 intercept = input('Enter Temperature Calibration Intercept: ');
 E0 = input('Enter pH Calibration E0: ');
 Ts = input('Enter pH Calibration Water temperature (K): ');

%go thru data, get Time, V1, V2, V4
 time=[];
 V1=[]; %pH sensor voltage [-150:150] mV
 V2=[]; %sealed positive battery voltage [-3.6 3.6] V
 V4=[]; %Wheatstone bridge thermistor output
 while isempty(strfind(fline,'"[End Reading]"'))
   fline=fgetl(fid);
   f=strfind(fline,',');
   try %won't work if you're on the last row
     time=[time; datenum([1900 1 1]) + ...
       str2num(fline(f(1)+2:f(2)-2))-2]; % **DON'T CHANGE TIME ZONE**
     V1=[V1; str2num(fline(f(2)+2:f(3)-2))];
     V2=[V2; str2num(fline(f(5)+2:f(6)-2))];
     V4=[V4; str2num(fline(f(11)+2:f(12)-2))];
   end
 end
 fclose(fid);
 clear fid fline f;

%battery voltage w/ scale factor applied
 Vbat=V2*VbatSF;

%Eq'n for thermistor resistance based on voltage divider and wheatstone bridge
 Col5=V4+(Vbat*(15000/(1e6+15000)));
 Col6=Col5*1e6;
 Col7=Vbat-Col5;
 Rtherm=Col6./Col7;
 clear Vbat Col5 Col6 Col7;

%Eq'n for temperature based on factory constants and measured resistance
 %Steinhart Constants
  SH_A = 1.06329736674527E-03;
  SH_B = 2.51377462346306E-04;
  SH_C = 2.55455247726963E-08;
 k = 1 ./ ( SH_A + (SH_B * log(Rtherm)) + (SH_C * (log(Rtherm).^3)) );
 durafetTemp = k - 273.15;
 clear Rtherm SH_A SH_B SH_C k;

%apply temp cal coeff's
 calTemp=(slope*durafetTemp)+intercept;
 %NaN imaginary values
  for n=1:length(calTemp)
    if ~isreal(calTemp(n))
      calTemp(n)=NaN;
    end
  end
  clear n;

%calculate pH
 ST=(R*(calTemp+273.15)*log(10)) ./ F;
 E0T=E0-0.001*((calTemp+273.15)-Ts);
 pH=(V1-E0T)./ST;
 clear ST EOT;

%check
 figure;
 subplot(211); plot(pH,'b-'); ylabel('pH');
 title([datestr(time(1),'yyyy.mm.dd HH:MM') ' - ' ...
   datestr(time(end),'yyyy.mm.dd HH:MM')]);
 subplot(212); plot(calTemp,'k-'); ylabel('cal''d temp');

%--%

%write new 'converted' file
 s=strfind(myfile,'raw');
 if ~isempty(s) %should have 'raw' in filename, replace with 'conv'
   myfileCal=[myfile(1:s-1) 'conv' myfile(s+3:end)];
 else %otherwise append 'conv' to end of filename.
   myfileCal=[myfile(1:end-4) '_conv.csv'];
 end
 fid=fopen([mydir myfileCal],'w+');
 clear s;

%first, write old header to new 'conv' file
 for n=1:length(header)
   fprintf(fid,'%s\r\n',header{n});
 end
 clear n;

%now, add new header info (comma-delimited)
 %define new header
  newHeader={''; ...
    ['Date Processed, ' datestr(now,'mm/dd/yyyy HH:MM')]; ...
    'Calibration Values Used:'; ...
    ['VbatSF, ' num2str(VbatSF) ', Battery voltage Scale Factor']; ...
    ['slope, ' num2str(slope)]; ...
    ['intercept, ' num2str(intercept)]; ...
    ['E0, ' num2str(E0) ', pH calibration coefficient ''E0''']; ...
    ['Ts, ' num2str(Ts) ', temperature of calibration (K)']; ...
    ['R, ' num2str(R) ', Gas Constant (J mol-1 K-1)']; ...
    ['F, ' num2str(F) ', Faraday Constant (C mol-1)']; ...
    ''; ...
    ['GMT Year, Month, Day, Hour, Minute, Second, Voltage 1, ' ...
      'Voltage 2, Voltage 4, Temperature (C), pH']};
 %write new header
  for n=1:length(newHeader)
    fprintf(fid,'%s\r\n',newHeader{n});
  end

 %write data rows
  for n=1:length(V1)
    fprintf(fid,'%s\r\n',[datestr(time(n),'yyyy,mm,dd,HH,MM,SS') ',' ...
      num2str(V1(n)) ',' num2str(V2(n)) ',' num2str(V4(n)) ',' ...
      num2str(calTemp(n)) ',' num2str(pH(n))]);
  end
  clear n;

 fclose(fid);

%save figure
 set(gcf,'PaperPositionMode','auto');
 saveas(gcf,[mydir myfileCal(1:end-4)],'png');
