% -- 
% Jeff Sevadjian
% Research Technician
% Monterey Bay Aquarium Research Institute
% 7700 Sandholdt Rd
% Moss Landing, CA 95039, USA
% ph: +1 831-775-1888

% Applies calibrations to DATALOG.mat file (parsed from OASIS controller
%  using 'OA2_201505_parse_datalog.m').
% Equations taken off makeOAwebQC_oa1.py
% Calibrations taken off OA_cal_config.xlsx
% Also does a fairly thorough manual removal of bad data.

%define buoy, deployment, data directory
 buoy = 'OA2';
 dep = '201505';
 mydir = ['//atlas.shore.mbari.org/OA_Moorings/deployment_data/' ...
   buoy '/' dep '/'];
 myfile = 'DATALOG';

%load raw cntl data
 load([mydir 'raw/cntl/' myfile '.mat']);

%cal info
 WScal = [0.058 15.0]; %WetStar factory cal CWO, SF
 pHcal = [datenum([2015 5 13 20 28 0; 2016 1 26 17 35 0])', -0.351839, 293.74; ...
          datenum([2016 1 26 17 35 0; 2016 9 27 16  4 0])', -0.388613, 294.17];
 CO2cal = [-25.2633 0.83015 1.0495 -0.0019443];
 METcal = [1 0]; %Air Temp slope, intercept

%tBeg, tEnd
 tBeg = datenum([2015 5 13 20 28 0]);
 tEnd = datenum([2016 9 27 16  4 0]);

%CTD
 f = find(MICROCAT(:,1)>tBeg & MICROCAT(:,1)<tEnd);
 CTD = MICROCAT(f,[1,2,5]);
 clear MICROCAT f;
 %QC salinity
  %pump failed
   f = find( CTD(:,1) > datenum([2016 3 6 1 59 0]) );
   CTD(f,3) = NaN;
   clear f;
  %variance filter
   win = 100; %window
   nst = 5; %number of stdev's above median
   sm = slidingmed(CTD(:,3), win);
   ss = slidingstdev(CTD(:,3), win);
   fs = find(abs(CTD(:,3) - sm) > nst*ss);
   CTD(fs,3) = NaN;
   clear win nst sm ss fs;
  %manual outliers
   f = find( CTD(:,1) > datenum([2015 9 30 16 14 0]) & ...
     CTD(:,1) < datenum([2015 10 16 3 59 0]) ); %prolonged dip in S; T vs S out-of-whack
   CTD(f,3) = NaN;
   clear f;

%convert Aanderaa O2 conc to seawater, and then to umol/kg
 f = find(AANDERAA_O2(:,1)<tBeg | AANDERAA_O2(:,1)>tEnd);
 AANDERAA_O2(f,:) = [];
 o2uMS0 = AANDERAA_O2(:,4);
 T = interp1(CTD(:,1),CTD(:,2),AANDERAA_O2(:,1));
 S = 33.5; %better to keep fixed then try to change halfway
 [~,sigt] = swstate(S,T,1);
 B0 = -6.24097e-3;
 B1 = -6.93498e-3;
 B2 = -6.90358e-3;
 B3 = -4.29155e-3;
 C0 = -3.11680e-7;
 Ts = log((298.15-T) ./ (273.15+T));
 o2uM = o2uMS0 .* exp( S.* (B0 + B1*Ts + B2*(Ts.^2) + B3*(Ts.^3)) + C0*(S.^2) );
 o2umolkg = o2uM*1000./(sigt+1000);
 O2 = [AANDERAA_O2(:,1) o2umolkg];
 clear AANDERAA_O2 f o2uMS0 T S sigt B0 B1 B2 B3 C0 Ts o2uM o2umolkg;
 %QC O2 (manual outliers)
  f1 = find( O2(:,1) > datenum([2016 1 26 17 30 0]) & ...
    O2(:,1) < datenum([2016 2 23 19 59 0]) ); %plumbing messed up
  f2 = find( O2(:,1) > datenum([2016 9 8 2 16 0]) & ...
    O2(:,1) < datenum([2016 9 15 14 9 0]) ); %not sure, near end of depl
  O2([f1; f2],2) = NaN;
  clear f1 f2;

%convert WetStar voltages to ug/L
 f = find(EXT_ANALOG3(:,1)<tBeg | EXT_ANALOG3(:,1)>tEnd);
 EXT_ANALOG3(f,:) = [];
 volts = EXT_ANALOG3(:,2);
 WSout = WScal(2) * (volts - WScal(1));
 Fluor = [EXT_ANALOG3(:,1) WSout];
 clear EXT_ANALOG3 WScal volts WSout;
 %QC fluor (manual outliers)
  f1 = find( O2(:,1) > datenum([2016 1 26 17 30 0]) & ...
    O2(:,1) < datenum([2016 2 23 19 59 0]) ); %plumbing messed up
  Fluor(f1,2) = NaN;
  clear f f1;

%convert pH sensor voltage 'V1' to pH units
% input T should be from CTD
 f = find(EXT_ANALOG2(:,1)<tBeg | EXT_ANALOG2(:,1)>tEnd);
 EXT_ANALOG2(f,:) = [];
 pHV1 = EXT_ANALOG2(:,2)/1000;
 T = interp1(CTD(:,1),CTD(:,2),EXT_ANALOG2(:,1));
 R = 8.31451; %gas constant (J mol-1 K-1)
 F = 96487; %converts between Coulombs and Faradays (1 F = 96487 C)
 ST = ( R * (T+273.15) * log(10) ) / F;
 pHconv = NaN(size(pHV1));
 for n=1:size(pHcal,1)
   f = find(EXT_ANALOG2(:,1) > pHcal(n,1) & EXT_ANALOG2(:,1) < pHcal(n,2));
   E0T = pHcal(n,3)-0.001 * ( (T(f)+273.15) - pHcal(n,4) );
   pHconv(f) = (pHV1(f) - E0T)./ST(f);
 end
 pH = [EXT_ANALOG2(:,1) pHconv];
 clear EXT_ANALOG2 pHcal pHV1 T R F ST pHconv n f E0T;
 %QC pH (manual outliers)
  f1 = find( pH(:,1) > datenum([2015 12 31 17 59 0]) & ...
    pH(:,1) < datenum([2016 1 26 18 28 0]) ); %batteries dying
  f2 = find( pH(:,1) > datenum([2016 8 21 1 43 0]) ); %batteries dying
  pH([f1; f2],2) = NaN;
  clear f1 f2;

%index CO2
 f = find(PCO2(:,1)<tBeg | PCO2(:,1)>tEnd);
 PCO2(f,:) = [];
 f = find(PCO2(1:end-5,6)==1 & PCO2(6:end,6==6));
 pCO2 = PCO2(f+1,1); %time-off
 CO2equil = PCO2(f+1,2:4);
 CO2zero = PCO2(f+3,2:4);
 CO2air = PCO2(f+5,2:4);
 clear PCO2 f;
%apply CO2 temperature-cal
% input T should be from CO2 system, not CTD
% no standard was used for this deployment
 %water
  xCO2 = CO2equil(:,3) - CO2zero(:,3);
  T = CO2equil(:,1);
  xCO2water = CO2cal(1) + CO2cal(2)*T + CO2cal(3)*xCO2 + CO2cal(4)*T.*xCO2;
  clear xCO2 T;
  %convert sea water xCO2 to pCO2
  % input T, S here should be from CTD
   T = interp1(CTD(:,1),CTD(:,2),pCO2(:,1));
   S = 33.5; %better to keep fixed then try to change halfway
   TempK = T + 273.15;
   VPWP = exp(24.4543 - 67.4509 * (100./T) - 4.8489 .* log(TempK/100));
   VPCorrWP = exp(-0.000544 * S);
   VPSWWP = VPWP .* VPCorrWP;
   VPFac = 1 - VPSWWP;
   pCO2(:,2) = xCO2water .* VPFac;
   clear xCO2water T S TempK VPWP VPCorrWP VPSWWP VPFac;
 %air
  xCO2 = CO2air(:,3) - CO2zero(:,3);
  T = CO2air(:,1);
  pCO2(:,3) = CO2cal(1) + CO2cal(2)*T + CO2cal(3)*xCO2 + CO2cal(4)*T.*xCO2;
  clear xCO2 T;

 %QC pCO2
  %water
   %manual outliers
    f1 = find( CO2equil(:,2) < 99 | CO2equil(:,2) > 103 ); %plugged after winter storms
    [~,f2] = min(abs(pCO2(:,1) - datenum([2015 11 25 17 6 47])));
    [~,f3] = min(abs(pCO2(:,1) - datenum([2015 11 25 18 6 48])));
    f4 = find( pCO2(:,1) > datenum([2015 12 6 15 6 0]) ); %jump in press, quest'ble after
    pCO2([f1; f2; f3; f4],2) = NaN;
    clear f1 f2 f3 f4;
   %compare pH and pCO2w
    %index pCO2 to pH
     dt = 15/1440;
     xt = round(pCO2(:,1)/dt)*dt;
     yt = round(pH(:,1)/dt)*dt;
     [ti,ia,ib] = intersect(xt,yt);
     x = pCO2(ia,2);
     y = pH(ib,2);
     clear dt xt yt ia ib;
    %calc theoretical relation
     sali=[33 34];
     tempi=[8 15];
     alk=2150.0 + 44.0.*(sali-31.25); %(GF)
    %CO2SYS
     sil=0;
     po4=1;
     pH_th=[7.5:0.01:8.5];
     for n=1:length(sali)
       A=CO2SYS(alk(n),pH_th,1,3,sali(n),tempi(n),NaN,0,NaN,...
         sil,po4,1,4,1);
       xPred(:,n)=A(:,4); %(uatm)
     end
     clear sali tempi alk sil po4 n A;
%     %plot
%      figure; hold on;
%      plot(xPred(:,1),pH_th,'-','Color',[.4 .4 .4],'LineWidth',2);
%      plot(xPred(:,2),pH_th,'k-','LineWidth',2);
%      plot(x,y,'b.');
    %manual outliers (none; already removed above)
    clear ti x y xPred pH_th;

  %air (manual outliers)
   f1 = find(CO2air(:,2) < 99);
   f2 = find( pCO2(:,1) > datenum([2015 12 22 21 6 0]) ); %seems to be OK for a while after
   pCO2([f1; f2],3) = NaN;
   clear f1 f2;
  clear CO2* EXT_ANALOG1;

%MET
 %correct PB200 air temp
  f = find(AIRMAR_PB200(:,1)>tBeg & AIRMAR_PB200(:,1)<tEnd);
  MET = AIRMAR_PB200(f,:);
  MET(:,3) = MET(:,3) * METcal(1) + METcal(2);
  clear METcal f AIRMAR_PB200;
 %QC wind speed
  f = find(MET(:,7) > 50);
  MET(f,7) = NaN;
  clear f;

%charging system 
 clear SPECIAL1;

%VR2C
 f = find(SPECIAL2(:,1)>tBeg & SPECIAL2(:,1)<tEnd);
 VR2C = SPECIAL2(f,:);
 clear SPECIAL2 f;

%save
 save([mydir 'QC/' myfile '_QC.mat'],'CTD', 'Fluor', 'MET', 'O2', 'pCO2', ...
   'pH', 'VR2C');