% -- 
% 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 'OA1_201401_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 = 'OA1';
 dep = '201401';
 mydir = ['//atlas.shore.mbari.org/OA_Moorings/deployment_data/' ...
   buoy '/' dep '/'];
 myfile = 'DATALOG';

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

%cal info
 WScal = [0.056, 15.4]; %WetStar factory cal CWO, SF
 pHcal = [datenum([2014 1 23 23 46 0; 2015 1  8 18 32 0])', -0.369461, 280.66; ...
          datenum([2015 1  8 18 32 0; 2015 9 16 17  0 0])', -0.371961, 283.44; ...
          datenum([2015 9 16 17  0 0; 2016 5  4 16 35 0])', -0.373990, 279.13; ...
          datenum([2016 5  3 16 36 0; 2016 7 20 16 20 0])', -0.351839, 293.74];
 CO2cal = [-14.3968, 0.4257, 1.0147, -0.00022448];
 METcal = [1 0]; %Air Temp slope, intercept

%tBeg, tEnd
 tBeg = datenum([2014 1 23 23 46 0]);
 tEnd = datenum([2016 7 20 16 20 0]);
 
%CTD
 f = find(MICROCAT(:,1)>tBeg & MICROCAT(:,1)<tEnd);
 CTD = MICROCAT(f,[1,2,5]);
 clear MICROCAT f;
 %QC salinity
  %fixed floor
   f = find(CTD(:,3) < 30);
   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
   [~,i] = min(abs(CTD(:,1) - datenum([2015 12 2 20 9 20])));
   CTD(i,3) = NaN;
   clear i;

%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 = interp1(CTD(:,1),CTD(:,3),AANDERAA_O2(:,1));
 [~,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([2014 2 6 20 46 0]) & ...
    O2(:,1) < datenum([2014 2 24 22 17 0]) ); %pump issues [diver bubbles/jelly]
  f2 = find( O2(:,1) > datenum([2014 8 30 0 16 0]) & ...
    O2(:,1) < datenum([2014 9 4 20 2 0]) );
  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( Fluor(:,1) < datenum([2014 2 25 1 32 0]) );
  f2 = find( Fluor(:,1) > datenum([2014 8 29 17 47 0]) & ...
    Fluor(:,1) < datenum([2014 9 4 19 1 0]) );
  Fluor([f1; f2],2) = NaN;
  clear f;

%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([2014 2 6 23 2 0]) & ...
    pH(:,1) < datenum([2014 2 25 0 1 0]) ); %pump issues [diver bubbles/jelly]
  f2 = find( pH(:,1) > datenum([2014 8 29 18 2 0]) & ...
    pH(:,1) < datenum([2015 1 8 19 1 0]) ); %battery death
  f3 = find( pH(:,1) > datenum([2015 6 30 3 47 0]) & ...
    pH(:,1) < datenum([2015 9 16 17 49 0]) ); %battery death
  f4 = find( pH(:,1) > datenum([2016 3 31 20 25 0]) ); %wires twisted inside can
  pH([f1; f2; f3; f4],2) = NaN;
  clear f1 f2 f3 f4;

%index CO2
 f = find(PCO2(:,1)<tBeg | PCO2(:,1)>tEnd);
 PCO2(f,:) = [];
 f = find(PCO2(1:end-5,9)==1 & PCO2(6:end,9)==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 = interp1(CTD(:,1),CTD(:,3),pCO2(:,1));
   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) < 100 | CO2equil(:,2) > 103 );
    f2 = find( pCO2(:,1) > datenum([2014 8 4 23 55 0]) & ...
      pCO2(:,1) < datenum([2014 9 4 15 1 0]) );
    f3 = find( pCO2(:,1) > datenum([2015 11 18 19 3 0]) ); %chain snagged on rock, filters plugged
    [~,f4] = min(abs(pCO2(:,1) - datenum([2015 5 19 4 54 33])));
    f5 = find(pCO2(:,2) < 0);
    pCO2([f1; f2; f3; f4; f5],2) = NaN;
    clear f1 f2 f3 f4 f5;
   %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) < 95);
   f2 = find(pCO2(:,3) < 0);
   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;
 %get rid of one NaN
  f = find(isnan(VR2C(:,2)));
  VR2C(f,:) = [];
  clear f;

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