
% This is for post-processing, not real-time processing.

% veh is SV2 or SV3
% dep is yyyymmdd of deployment date
% pCO2cal is 5-element vector [A B C D E] of pCO2 cal coeff's
% pHcal is 2-element vector [E0 Ts]
% BBcal is 3x2 matrix from WetLabs *.dev file:
%  rows are per wavelength
%  columns are [scale-factor dark-value]

veh='SV3';
dep='20150501';
pCO2cal=[-30.2789 0.43639 1.092 -0.002425 1];
pHcal=[-0.397897 293.12];
BBcal=[1.08E-05	47;	3.95E-06	46;	0.0121	49];

if ismac
  pre='/Volumes/';
elseif ispc
  pre='//atlas.shore.mbari.org/';
else
  pre='/mbari/';
end

%working directory
 mydir=[pre 'WaveGlider/deployment_data/' veh '/' dep '/'];
 myfile=[veh '_' dep];

%add path to m-files/functions called from here
 addpath(genpath([pre 'WaveGlider/code/']));
 clear pre;

%import csv data into Matlab
 fid=fopen([mydir myfile '.csv']);

%skip past header
 hd=1; ft=[];
 while hd
   fline=fgetl(fid);
   s=strfind(fline,',');
   hd=strcmp(fline(s(1)-1),'*');
   ft=[ft; ftell(fid)];
 end
 clear fline s hd;

%data
 dn=datenum([1970 1 1]);
 ct=0; ctWx=0; ctCTD=0; ctCO2=0; ctBB=0;
 fseek(fid,ft(end-1),'bof'); %go back one line
 while ~feof(fid)
   fline=fgetl(fid);
   f=strfind(fline,',');
   s=str2num(fline(f(1)+1:end));
   if ~isempty(strfind(fline,'Waveglider')) && length(s)==11
       ct=ct+1;
       glider.lat(ct)=s(1);
       glider.lon(ct)=s(2);
       glider.time(ct)=dn+s(3)/86400/1000; %ms since epoch
       glider.tempSub(ct)=s(9);
   elseif ~isempty(strfind(fline,'CTD')) && length(s)==11
       ctCTD=ctCTD+1;
       gliderCTD.lat(ctCTD)=s(1);
       gliderCTD.lon(ctCTD)=s(2);
       gliderCTD.time(ctCTD)=dn+s(3)/86400/1000;
       gliderCTD.cond(ctCTD)=s(4);
       gliderCTD.o2(ctCTD)=s(5);
       gliderCTD.o2Hz(ctCTD)=s(6);
       gliderCTD.o2sol(ctCTD)=s(7);
       gliderCTD.press(ctCTD)=s(8);
       gliderCTD.sal(ctCTD)=s(9);
       gliderCTD.temp(ctCTD)=s(10);
       gliderCTD.sid(ctCTD)=s(11); %serial num?
   elseif ~isempty(strfind(fline,'pCO2')) && length(s)==38
       ctCO2=ctCO2+1;
       gliderCO2.lat(ctCO2)=s(1);
       gliderCO2.lon(ctCO2)=s(2);
       gliderCO2.time(ctCO2)=dn+s(3)/86400/1000; %'sample time'
%        gliderCO2.time(ctCO2)=dn+s(4)/86400/1000; %'time'
       gliderCO2.pHtime(:,ctCO2)=gliderCO2.time(ctCO2)+...
         [0:10:50]/1440;
       gliderCO2.pHV1(:,ctCO2)=s(30:35)/1000;
      %Equil
       gliderCO2.EquilPumpOnCO2(ctCO2)=s(5);
       gliderCO2.EquilPumpOnTemp(ctCO2)=s(6);
       gliderCO2.EquilPumpOnPress(ctCO2)=s(7);
       gliderCO2.EquilPumpOffCO2(ctCO2)=s(8);
       gliderCO2.EquilPumpOffTemp(ctCO2)=s(9);
       gliderCO2.EquilPumpOffPress(ctCO2)=s(10);
       gliderCO2.EquilPumpOffHum(ctCO2)=s(11);
      %Zero
       gliderCO2.ZeroPumpOnCO2(ctCO2)=s(12);
       gliderCO2.ZeroPumpOnTemp(ctCO2)=s(13);
       gliderCO2.ZeroPumpOnPress(ctCO2)=s(14);
       gliderCO2.ZeroPumpOffCO2(ctCO2)=s(15);
       gliderCO2.ZeroPumpOffTemp(ctCO2)=s(16);
       gliderCO2.ZeroPumpOffPress(ctCO2)=s(17);
      %Air
       gliderCO2.AirPumpOnCO2(ctCO2)=s(18);
       gliderCO2.AirPumpOnTemp(ctCO2)=s(19);
       gliderCO2.AirPumpOnPress(ctCO2)=s(20);
       gliderCO2.AirPumpOffCO2(ctCO2)=s(21);
       gliderCO2.AirPumpOffTemp(ctCO2)=s(22);
       gliderCO2.AirPumpOffPress(ctCO2)=s(23);
       gliderCO2.AirPumpOffHum(ctCO2)=s(24);
      %Standard
       gliderCO2.StandardFlowOnPress(ctCO2)=s(25);
       gliderCO2.StandardFlowOffCO2(ctCO2)=s(26);
       gliderCO2.StandardFlowOffTemp(ctCO2)=s(27);
       gliderCO2.StandardFlowOffPress(ctCO2)=s(28);
       gliderCO2.StandardFlowOffHum(ctCO2)=s(29);
      %Can Humidity
       gliderCO2.CanHumidity(ctCO2)=s(36);
      %Two extra columns of only zeros (?)
   elseif ~isempty(strfind(fline,'Weather')) && length(s)==10
       ctWx=ctWx+1;
       gliderWx.lat(ctWx)=s(1);
       gliderWx.lon(ctWx)=s(2);
       gliderWx.time(ctWx)=dn+s(3)/86400/1000;
       gliderWx.wdirAvg(ctWx)=s(4);
       gliderWx.wspdAvg(ctWx)=s(5)*0.514444; %knots to m/s
       gliderWx.wspdMax(ctWx)=s(6)*0.514444;
       gliderWx.airPress(ctWx)=s(7);
       gliderWx.airTemp(ctWx)=s(10);
   elseif ~isempty(strfind(fline,'ECO-Puck'))
       s=str2num(fline(f(1)+1:f(end-1))); %last col is garb txt
       ctBB=ctBB+1;
       gliderBB.lat(ctBB)=s(1);
       gliderBB.lon(ctBB)=s(2);
       gliderBB.time(ctBB)=dn+s(3)/86400/1000;
       gliderBB.beta1(ctBB)=BBcal(1,1) * (s(5) - BBcal(1,2));
       gliderBB.beta2(ctBB)=BBcal(2,1) * (s(7) - BBcal(2,2));
       gliderBB.chl(ctBB)=BBcal(3,1) * (s(9) - BBcal(3,2));
       gliderBB.wavel(:,ctBB)=s([4 6 8]); %is constant
       gliderBB.therm(ctBB)=s(10); %thermistor (not used)
   end
 end
 fclose(fid);
 clear dn fid ft n ct* fline f s ans;

%get unique time stamps only
 %glider
  fn=fieldnames(glider);
  [~,i]=unique(glider.time);
  for n=1:length(fn)
    e=eval(['glider.' fn{n}]);
    glider=setfield(glider,fn{n},e(i));
  end
 %gliderCO2
  if exist('gliderCO2')
    [~,i]=unique(gliderCO2.time);
    fn=fieldnames(gliderCO2);
    for n=1:length(fn)
      e=eval(['gliderCO2.' fn{n}]);
      gliderCO2=setfield(gliderCO2,fn{n},e(:,i));
    end
  end
 %gliderCTD
  if exist('gliderCTD')
    [~,i]=unique(gliderCTD.time);
    fn=fieldnames(gliderCTD);
    for n=1:length(fn)
      e=eval(['gliderCTD.' fn{n}]);
      gliderCTD=setfield(gliderCTD,fn{n},e(i));
    end
  end
 %gliderWx
  if exist('gliderWx')
    [~,i]=unique(gliderWx.time);
    fn=fieldnames(gliderWx);
    for n=1:length(fn)
      e=eval(['gliderWx.' fn{n}]);
      gliderWx=setfield(gliderWx,fn{n},e(i));
    end
  end
 %gliderBB
  if exist('gliderBB')
    [~,i]=unique(gliderBB.time);
    fn=fieldnames(gliderBB);
    for n=1:length(fn)
      e=eval(['gliderBB.' fn{n}]);
      gliderBB=setfield(gliderBB,fn{n},e(:,i));
    end
  end
 clear i fn n e;

%NaN placeholder CO2 standard values = 0
 if exist('gliderCO2')
   f=find(gliderCO2.StandardFlowOffCO2==0);
   gliderCO2.StandardFlowOnPress(f)=NaN;
   gliderCO2.StandardFlowOffCO2(f)=NaN;
   gliderCO2.StandardFlowOffTemp(f)=NaN;
   gliderCO2.StandardFlowOffPress(f)=NaN;
   gliderCO2.StandardFlowOffHum(f)=NaN;
   clear f;
 end

%derive CTD variables
 if exist('gliderCTD')
  %calc sigma-t
   [~,gliderCTD.sigt]=swstate(gliderCTD.sal,gliderCTD.temp,...
     gliderCTD.press);
  %convert oxygen from ml/liter to micromoles/kg (ryjo)
   gliderCTD.o2umolkg = (gliderCTD.o2 * 1.4276) .* ...
     (1e6./((gliderCTD.sigt+1000)*32));
  %calc O2 perc saturation?
 end

%derive bb from beta
if exist('gliderBB')
  %get beta and bb from pure sw (Morel, 1974)
   th=117*pi/180; %all these WL sensors measure at 117 degrees
   lambda=[mode(gliderBB.wavel(1,:)) mode(gliderBB.wavel(2,:))];
   d=0.09;
   if exist('gliderCTD')
     S=repmat(interp1(gliderCTD.time,gliderCTD.sal,...
       gliderBB.time),2,1);
     S(isnan(S))=33.5; %in case CTD got switched off, etc.
   else
     S=33.5; %assume constant
   end
   le=length(gliderBB.time);
   term1=repmat(1.38 * ((lambda/500).^-4.32),le,1)';
   term3=repmat(1 + ((cos(th)^2)*(1-d)/(1+d)),le,2)';
   beta_w = term1 .* ((1+(0.3*S/37))*10^-4) .* term3;
   bb_w = 8*pi/3*beta_w*(2+d)/(1+d);
   clear th d S le term1 term3;
  %calc betap, bbp, and estimate bb
  % ONLY KEEPING BB FOR NOW, COULD ADD BETAP, BBP LATER
   X=1.1;
   %470
    betap1=gliderBB.beta1-beta_w(1,:);
    bbp1=2*pi*X*betap1;
    gliderBB.bb1=bbp1+bb_w(1,:);
   %650 | 700
    betap2=gliderBB.beta2-beta_w(2,:);
    bbp2=2*pi*X*betap2;
    gliderBB.bb2=bbp2+bb_w(2,:);
   clear bb_w beta_w X betap* bbp*;
  %rename '1','2' according to wavel
   gliderBB=setfield(gliderBB,['beta' num2str(lambda(1))],...
     gliderBB.beta1);
   gliderBB=setfield(gliderBB,['beta' num2str(lambda(2))],...
     gliderBB.beta2);
   gliderBB=setfield(gliderBB,['bb' num2str(lambda(1))],...
     gliderBB.bb1);
   gliderBB=setfield(gliderBB,['bb' num2str(lambda(2))],...
     gliderBB.bb2);
   gliderBB=rmfield(gliderBB,{'beta1','beta2','bb1','bb2'});
   clear lambda;
end

%derive pCO2 and pH, calc CO2 flux
 if exist('gliderCO2') %big if statement here

  %derive pH
   %get CTD T at pH sample times
    if exist('gliderCTD')
     %temp
      tempi=interp1(gliderCTD.time,gliderCTD.temp,...
        gliderCO2.pHtime);
      [~,b]=find(isnan(tempi));
      u=unique(b);
      tempi(:,u)=repmat(gliderCO2.EquilPumpOffTemp(u)-0.5,6,1);
     %sal (for later)
      sali=interp1(gliderCTD.time,gliderCTD.sal,gliderCO2.pHtime);
      sali(isnan(sali))=33.5;
    else
      tempi=interp1(gliderCO2.time,...
        gliderCO2.EquilPumpOffTemp-0.5,gliderCO2.pHtime);
      sali=33.5;
    end
    clear b u;
   %derive
    R=8.31451; %gas constant (J mol-1 K-1)
    F=96487; %converts between Coulombs and Faradays (1 F = 96487 C)
    ST=(R*(tempi+273.15)*log(10)) ./ F;
    E0T=pHcal(1)-0.001*((tempi+273.15)-pHcal(2));
    gliderCO2.pH=(gliderCO2.pHV1-E0T)./ST; 
    clear R F ST E0T;

  %estimate pCO2 from pH
   %estimate alk from sal
    alk=2150.0 + 44.0*(sali-31.25); %(GF)
   %CO2SYS
    sil=0;
    po4=1;
    ne=numel(gliderCO2.pH);
    A=CO2SYS(reshape(alk,1,ne),reshape(gliderCO2.pH,1,ne),1,3,...
      reshape(sali,1,ne),reshape(tempi,1,ne),...
      NaN,0,NaN,sil,po4,1,4,1);
    gliderCO2.pCO2waterEst=reshape(A(:,4),...
      size(gliderCO2.pHtime,1),size(gliderCO2.pHtime,2)); %(uatm)
    clear tempi sali ne A;

  %derive pCO2
  %get CTD T,S at pCO2 sample times
   if exist('gliderCTD')
     tempi=interp1(gliderCTD.time,gliderCTD.temp,gliderCO2.time);
     f=find(isnan(tempi));
     tempi(f)=gliderCO2.EquilPumpOffTemp(f)-0.5;
     sali=interp1(gliderCTD.time,gliderCTD.sal,gliderCO2.time);
     sali(isnan(sali))=33.5;
     clear f;
   else
     tempi=gliderCO2.EquilPumpOffTemp-0.5;
     sali=33.5;
    end
  %derive  
   gliderCO2.pCO2air=calcCO2(pCO2cal,...
     gliderCO2.AirPumpOffCO2,...
     gliderCO2.ZeroPumpOffCO2,...
     gliderCO2.AirPumpOffTemp,tempi,sali);
   gliderCO2.pCO2water=calcCO2(pCO2cal,...
     gliderCO2.EquilPumpOffCO2,...
     gliderCO2.ZeroPumpOffCO2,...
     gliderCO2.EquilPumpOffTemp,tempi,sali);
   gliderCO2.pCO2stnd=calcCO2(pCO2cal,...
     gliderCO2.StandardFlowOffCO2,...
     gliderCO2.ZeroPumpOffCO2,...
     gliderCO2.StandardFlowOffTemp,tempi,sali);

  %estimate pH from pCO2
   %estimate alk from sal
    alk=2150.0 + 44.0*(sali-31.25); %(GF)
   %CO2SYS
    A=CO2SYS(alk,gliderCO2.pCO2water,1,4,...
      sali,tempi,NaN,0,NaN,sil,po4,1,4,1);
    gliderCO2.pHest=A(:,3)';
    clear sil po4 alk A;

  %calc pCO2 flux (from Reiko, via Dorota)
   if exist('gliderWx') %flux is a fn of wind spd
    %get wind speed at pCO2 sample times
     wSpdi=interp1(gliderWx.time,gliderWx.wspdAvg,gliderCO2.time);
    %normalize to 10 m height
     f=find(~isnan(wSpdi) & wSpdi<100);
     wSpdi_10=NaN(size(wSpdi));
     [~,wSpdi_10(f)]=cdnlp(wSpdi(f),1);
     clear f;
    %calc Schmidt number (Sc) from Wanninkhof 1992
     Sc = 2073.1 - (125.62 * tempi) + ...
       (3.6276 * (tempi.^2)) - (0.043219 * (tempi.^3));
    %calc k from Wanninkhof 1992
     k = 0.31 * (wSpdi_10.^2) .* ((Sc/660).^-0.5); %wind spd is m/s
     clear wSpdi wSpdi_10 Sc;
    %calc solubility from Weiss 1974 APPENDIX IN WANNINKHOF 1992 IS KG
     T = tempi + 273.15; %T in K, S in ppt
     K0 = -58.0931 + (90.5069 * (100./T)) + ...
       (22.2940 * log(T/100)) + ...
       ( sali .* (0.027766 - (0.025888 * (T/100)) + (0.0050578 * ...
       (T/100).^ 2)) );
     CO2Sol = exp(K0); %in mol L-1 atm-1
     clear tempi sali T K0;
   %calc flux
    gliderCO2.delta_pco2=gliderCO2.pCO2water-gliderCO2.pCO2air;
   %convert from cm mol uatm hr-1 m-3 atm-1 to mol m-2 yr-1 
    gliderCO2.flux = k .* CO2Sol .* gliderCO2.delta_pco2 * ...
      24 * 365.25 * (10^-6) * (10^-2); %JS
    clear CO2Sol k;
   end %flux

 end %end of big if statement

%save mat file to working directory
 w=whos('glider*');
 save([mydir myfile '.mat'],w.name);
