clear all
File=dir('original2016-07-26_081758.txt'); 
File={File.name};
for ia=1:length(File)

    fileID = fopen(File{ia});
    
    %read in data
    data = textscan(fileID,'%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f','delimiter',','); % looks for numbers (%f) and text strings (%s), mostly separated by commas 
    
    fclose(fileID);
    
    %Extract columns to use:

    raw_doy              = data{1};
    raw_co2          = data{2};
    raw_co2_abs      = data{3};
    raw_h2o          = data{4};
    raw_h2o_abs      = data{5};
    raw_valve        = data{6};
    raw_cell_t       = data{7};
    raw_cell_p       = data{8};
    raw_820          = data{9};
    raw_teq          = data{10};% equilibrator
    raw_peq          = data{11};
    raw_temp1        = data{12};
    raw_ph1          = data{13};
    raw_temp2        = data{14};
    raw_cond2        = data{15};
    raw_sfromc       = data{16};%thermosalinograph
    raw_tsg_t        = data{17};
    raw_tsg_c        = data{18};
    raw_tsg_s        = data{19};
    raw_bs           = data{20};
    raw_bnumber      = data{21};
    raw_cpudate      = data{22};
    raw_cputime      = data{23};
    raw_lat          =  data{24};
    raw_lon          = data{25};
    raw_speed        = data{26};
    
end
dn_all = datenum(2016,0,raw_doy,0,0,0); % extracts the timestamp and puts it in a numerical format.
dv_all = datevec(dn_all);
 X=1:5;
Data_all=[raw_doy raw_co2 raw_co2_abs raw_h2o raw_h2o_abs raw_valve raw_cell_t raw_cell_p raw_820 raw_teq raw_peq raw_temp1 raw_ph1 raw_temp2 raw_cond2 raw_sfromc raw_tsg_t raw_tsg_c raw_tsg_s raw_bs raw_bnumber raw_cpudate raw_cputime raw_lat raw_lon raw_speed];


%% to flag/ignore the first 5 readings at every valve change:
ZZ=find(Data_all(:,6) == 0); % removes unwanted readings where the valve is '0'
Data_all(ZZ,:)=[];
DN_std=datenum(2016,0,Data_all(:,1),0,0,0);

isample = find(Data_all(:,6) == 1); % finds the indices where valve position=1
inewsample = find(diff(isample)~=1); % indices of when the valve changes position
iatm = find(Data_all(:,6) == 2); % same methodology as isample
inewatm = find(diff(iatm)~=1); 
b=[1];
inewatm=[b;inewatm;length(iatm)];
inewsample=[b;inewsample];
istd1 = find(Data_all(:,6) == 3); 
inewstd1 = find(diff(istd1)~=1); 
inewstd1=[inewstd1;length(istd1)];
istd2 = find(Data_all(:,6) == 4);
inewstd2 = find(diff(istd2)~=1);
inewstd2=[inewstd2;length(istd2)];
istd3 = find(Data_all(:,6) == 5);
inewstd3 = find(diff(istd3)~=1);
inewstd3=[inewstd3;length(istd3)];
residual = ones(length(inewstd3),12);
  inewstd1(2,:)=[];
  inewstd2(2,:)=[];
  inewstd3(2,:)=[];

CO2=Data_all(:,2);

% The following loop disregards the first 7 values of when the valve
% changes position, finds the average of each standard per cycle of
% calibration, as well as the standard deviation and the timestamp.
for ii = 1:length(inewstd1)
    avg_std(ii,1) = nanmean(CO2(istd1(inewstd1(ii)-15):istd1(inewstd1(ii)))); % disregards the first 7 vaules when the valve has changed to 3 and gets the average % if you want to subtract the frist 8 (instead of 7) then decrease the number subtracted (i.e. 15 to 14)
    stdev_std(ii,1) = nanstd(CO2(istd1(inewstd1(ii)-15):istd1(inewstd1(ii))));
    avg_dn(ii,1) = nanmean(DN_std(istd1(inewstd1(ii)-15):istd1(inewstd1(ii))));%  disregards the first 5 vaules when the valve has changed to 3 and gets the stdev
    avg_std(ii,2) = nanmean(CO2(istd2(inewstd2(ii)-15):istd2(inewstd2(ii))));
    stdev_std(ii,2) = nanstd(CO2(istd2(inewstd2(ii)-15):istd2(inewstd2(ii))));
    avg_dn(ii,2) = nanmean(DN_std(istd2(inewstd2(ii)-15):istd2(inewstd2(ii))));
    avg_std(ii,3) = nanmean(CO2(istd3(inewstd3(ii)-15):istd3(inewstd3(ii))));
    stdev_std(ii,3) = nanstd(CO2(istd3(inewstd3(ii)-15):istd3(inewstd3(ii))));
    avg_dn(ii,3) = nanmean(DN_std(istd3(inewstd3(ii)-15):istd3(inewstd3(ii))));
end

AVG_time=mean(avg_dn,2);

std_values(1:93,1:3) = 1;% std_values is a matrix that contains the standard gas pCO2 readings
std_values(:,1) = 202.7;
std_values(:,2) = 800.0;
std_values(:,3) = 1503.0;

sz=size(avg_std);
% linear fit of pCO2 readings from valves 3-5, versus the standard values
% of the gases
for ic =1:sz(:,1)
   slope_intercept(ic,:) = polyfit(std_values(ic,1:2), avg_std(ic,1:2),1); % calculates the slope and y-intercept
end
format long
standards = [AVG_time avg_std stdev_std slope_intercept];
% Write the standards matrix to a .csv file
dlmwrite('standards_0726_1.csv', standards, 'delimiter', ',', 'precision', 15); 
