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

% Loads Wave Glider .csv data from //atlas/WaveGlider/deployment_data/
%  and makes some very basic plots.

%define vehicle, deployment, data directory
 veh = 'wgTiny';
 dep = '20170412';
 mydir = ['//atlas.shore.mbari.org/WaveGlider/deployment_data/' ...
    veh '/' dep '/realTime/'];
%  mydir = ['//atlas.shore.mbari.org/WaveGlider/deployment_data/' ...
%     veh '/' dep '/QC/'];
%  mydir = ['//atlas.shore.mbari.org/WaveGlider/deployment_data/' ...
%     veh '/' dep '/postCruise/'];

%read-in .csv files
 ctd = csvread([mydir 'CTD.csv'],1,0);
%  ctd1 = csvread([mydir 'CTD1.csv'],1,0);
%  ctd2 = csvread([mydir 'CTD2.csv'],1,0);
% ctdd = csvread([mydir 'CTDd.csv'],1,0);
 co2 = csvread([mydir 'CO2.csv'],1,0);
 ep = csvread([mydir 'EcoPuck.csv'],1,0);
 met = csvread([mydir 'MET.csv'],1,0);
 metw = csvread([mydir 'METwind.csv'],1,0);
 pos = csvread([mydir 'Position.csv'],1,0);
 dis = csvread([mydir 'Distance.csv'],1,0);

%replace zeros w/ NaNs
 ctd(ctd==0) = NaN;
%  ctd1(ctd1==0) = NaN;
%  ctd2(ctd2==0) = NaN;
 co2(co2==0) = NaN;
 ep(ep==0) = NaN;
 met(met==0) = NaN;
 metw(metw==0) = NaN;
 pos(pos==0) = NaN;
 dis(dis==0) = NaN;

%convert time in Unix to Matlab datenum
 dn = datenum([1970 1 1]);
 ctd(:,1) = dn + ctd(:,1)/86400/1000;
%  ctd1(:,1) = dn + ctd1(:,1)/86400/1000;
%  ctd2(:,1) = dn + ctd2(:,1)/86400/1000;
%  ctdd(:,1) = dn + ctdd(:,1)/86400/1000;
 co2(:,1) = dn + co2(:,1)/86400/1000;
 ep(:,1) = dn + ep(:,1)/86400/1000;
 met(:,1) = dn + met(:,1)/86400/1000;
 metw(:,1) = dn + metw(:,1)/86400/1000;
 pos(:,1) = dn + pos(:,1)/86400/1000;
 dis(:,1) = dn + dis(:,1)/86400/1000;
 clear dn;

%interpolate to get lon/lat for each sensor
 ctd(:,5) = interp1(pos(:,1),pos(:,2),ctd(:,1));
 ctd(:,6) = interp1(pos(:,1),pos(:,3),ctd(:,1));
%  ctd1(:,5) = interp1(pos(:,1),pos(:,2),ctd1(:,1));
%  ctd1(:,6) = interp1(pos(:,1),pos(:,3),ctd1(:,1));
%  ctd2(:,5) = interp1(pos(:,1),pos(:,2),ctd2(:,1));
%  ctd2(:,6) = interp1(pos(:,1),pos(:,3),ctd2(:,1));
%  ctdd(:,5) = interp1(pos(:,1),pos(:,2),ctdd(:,1));
%  ctdd(:,6) = interp1(pos(:,1),pos(:,3),ctdd(:,1));
 co2(:,end+1) = interp1(pos(:,1),pos(:,2),co2(:,1));
 co2(:,end+1) = interp1(pos(:,1),pos(:,3),co2(:,1));
 ep(:,end+1) = interp1(pos(:,1),pos(:,2),ep(:,1));
 ep(:,end+1) = interp1(pos(:,1),pos(:,3),ep(:,1));
 met(:,end+1) = interp1(pos(:,1),pos(:,2),met(:,1));
 met(:,end+1) = interp1(pos(:,1),pos(:,3),met(:,1));
 metw(:,end+1) = interp1(pos(:,1),pos(:,2),metw(:,1));
 metw(:,end+1) = interp1(pos(:,1),pos(:,3),metw(:,1));

%simple time-series plots
 figure; plot(ctd1(:,1),ctd1(:,3),'b.-');
% figure; plot(ctd2(:,1),ctd2(:,2),'b.');
 figure; plot(ep(:,1),ep(:,4),'b.-');
 figure; plot(co2(:,1),co2(:,2),'b.-');
 figure; plot(co2(:,1),co2(:,4),'b.-');

%map plot
 figure; hold on;
 plot(pos(:,2),pos(:,3),'b.-');

%colored map plot
 figure; 
%  scatter(ctdd(:,5),ctdd(:,6),3,ctdd(:,2),'filled');
 scatter(ctd1(:,6),ctd1(:,7),3,ctd1(:,2),'filled');
 colorbar;

%pCO2 vs pH
 %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
  x = deNaN(co2(:,1),co2(:,3));
  y = co2(:,2);
  figure;
  plot(x,y,'.');
  hold on;
  plot(xPred(:,1),pH_th,'-','Color',[.4 .4 .4],'LineWidth',2);
  plot(xPred(:,2),pH_th,'k-','LineWidth',2);
  xlabel('pCO2 [ppm]');
  ylabel('pH');
  scatter(x,y,3,co2(:,1),'filled');
  colormap(jet);
