 % after running the averaging and Mscript scripts, use this script to
% correct for the averaged readings by using the slopes and intercepts

%Averaged data
clear all

File=dir('AVGS_ST.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','delimiter',','); % looks for numbers (%f) and text strings (%s), mostly separated by commas 
    % extract columns to use
    %data averaged to 2 minutes
    raw_valve=data{1};
    avg_doy=data{2};
    avg_co2=data{3};
    avg_pH=data{4};
    avg_temp=data{5};
    avg_sal=data{6};
    avg_lat=data{7};
    avg_lon=data{8};
    fclose(fileID);
end

%now to extract the standard data

File2=dir('ALLSTDS.txt');
File2={File2.name};
for ia=1:length(File2)

    fileID = fopen(File2{ia});
    
    %read in data
    data = textscan(fileID,'%f%f%f%f%f%f%f%f%f','delimiter',',');
    doy=data{1};
    Avg_stdv3=data{2};
    Avg_stdv4=data{3};
    Avg_stdv5=data{4};
    stdev_stdv3=data{5};
    stdev_stdv4=data{6};
    stdev_stdv5=data{7};
    slope=data{8};
    intercept=data{9};
end
avg_std=[Avg_stdv3 Avg_stdv4 Avg_stdv5];
stdev_std=[stdev_stdv3 stdev_stdv4 stdev_stdv5];
ALLCO2=[ raw_valve avg_doy avg_co2 avg_pH avg_temp avg_lat avg_lon];
ALLcalib=[ avg_std stdev_std slope intercept];
%dv_all = datevec(avg_doy);
%Date=datetime(dv_all);
%Dn=datenum(Date);
Dn=avg_doy;
isample= find(raw_valve == 1);
avg_carbon=avg_co2(isample);
iatm= find(raw_valve == 2);
avg_atmos=avg_co2(iatm);
%%

IDXX=knnsearch(doy,avg_doy);% uses k-nearest neighbor to find the closest values to apply the calibration
slope_intercept=[slope(IDXX) intercept(IDXX)];
num_rows = size(avg_carbon,1);
for ie = 1:1:num_rows
    for ig = 1:1:1
         % uses slope_intercept to predict y and prevents you from having to type yfit =  p(1) * x + p(2)
        cor_co2(ie,ig) = polyval(slope_intercept(ie,:),avg_carbon(ie,ig)); 
        
    end    
end
num_rowss = size(avg_atmos,1);
for iee = 1:1:num_rowss
    for igg = 1:1:1
         % uses slope_intercept to predict y and prevents you from having to type yfit =  p(1) * x + p(2)
        cor_atmos(iee,igg) = polyval(slope_intercept(iee,:),avg_atmos(iee,igg)); 
        
    end    
end

% for ir=1:length(cor_atmos) %getting R-squared
%      yresid = avg_co2(ir,:) - cor_co2(ir,:);
%      SSresid = sum(yresid.^2);
%      SStotal = (length(avg_co2(ir,:))-1) * var(avg_co2(ir,:));
%   rsq_pco2(ir) = 1 - SSresid/SStotal;
% end
    yresid = avg_carbon - cor_co2;
    SSresid = sum(yresid.^2);
    SStotal = (length(avg_carbon)-1) * var(avg_carbon);
    rsq = 1 - SSresid/SStotal;




%%

% Dnc=Dn(isample);
% Dna=Dn(iatm);
% AAA=datestr(Dn,'mm/dd/yyyy');
% % plotyy(Dn,avg_co2,Dn,slope_intercept(:,1),'plot','scatter')
% % hold on
% % plotyy(Dn,cor_co2,Dn,slope_intercept(:,1),'plot','scatter')
% % datetick('x','mm/dd','keepticks');
% hold on
% scatter(Dna,cor_atmos)
% fittedXa=linspace(min(Dna), max(Dna), 200);
% fita=polyfit(Dna,cor_atmos,1);
% vala=polyval(fita,fittedXa)
% plot(fittedXa,vala,'linewidth',5)
% hold off
% xlabel ('time (days)')
% ylabel patm
% title ('patm anomalies over the whole cruise')
% datetick('x','mm/dd','keepticks');
% figure(2)
% hold on
% 
% scatter(Dnc,cor_co2)
% fittedX = linspace(min(Dnc), max(Dnc), 200);
% fit1=polyfit(Dnc(1:2500),cor_co2(1:2500),1);
% val1=polyval(fit1,fittedX(1:50))
% fit2=polyfit(Dnc(2501:5000),cor_co2(2501:5000),1);
% val2=polyval(fit2,fittedX(51:100));
% fit3=polyfit(Dnc(5001:7500),cor_co2(5001:7500),1);
% val3=polyval(fit3,fittedX(101:150));
% fit4=polyfit(Dnc(7501:end),cor_co2(7501:end),1);
% val4=polyval(fit4,fittedX(151:200));
% 
% plot(fittedX(1:50),val1,'linewidth',5)
% plot(fittedX(51:100),val2,'linewidth',5)
% plot(fittedX(101:150),val3,'linewidth',5)
% plot(fittedX(151:200),val4,'linewidth',5)
% hold off
% 
% datetick('x','mm/dd','keepticks');
% % model shows a big discrepancy in co2 values in the first 24 hours
% 
% 
% %legend('corrected CO2','corrected CO2')
% xlabel ('time (days)')
% ylabel pCO2
% title ('pCO2 anomalies over the whole cruise')
% 
% hold off
% % line([98 98], [min(cor_co2) max(cor_co2)]);
% % line([392 392], [min(cor_co2) max(cor_co2)]);
% % line([196 196], [min(cor_co2) max(cor_co2)]);
% % line([294 294], [min(cor_co2) max(cor_co2)]);
% 
% figure (4)
% hold on
% scatter(Dn,slope_intercept(:,1))
% ylabel slope
% datetick('x','mm/dd','keepticks');
% figure(5)
% scatter(Dn,slope_intercept(:,2))
% datetick('x','mm/dd','keepticks');
% xlabel time
% ylabel intercept
% DDD=datestr(Date);
% hold off
% %dlmwrite('corr.csv', all, 'delimiter', ',', 'precision', 15); 
% %%
% [hAx,hLine1,hLine2]=plotyy(Dnc,cor_co2,Dna,cor_atmos)
% ylabel(hAx(1),'pCO2') % left y-axis
% ylabel(hAx(2),'patm') % right y-axis
% datetick('x','mm/dd','keepticks');
% 
% %%
%   mdl=fitlm(avg_co2,cor_co2)
%  resid=mdl.Residuals.Raw;
%  find((mdl.Diagnostics.CooksDistance)>50*mean(mdl.Diagnostics.CooksDistance));
%  plotResiduals(mdl,'fitted')
%  meanval=mean(cor_co2);
%  stdDev=std(cor_co2);
%  outliers = abs(cor_co2-meanval) > (5 * stdDev); 
%  XXX=find(outliers>0)
%   cor_co2(XXX)=[];
%   cor_atmos(XXX)=[];
%   Dna(XXX)=[];
%   Dnc(XXX)=[];
% 
% % %  avg_co2(XXX)=[];
% %  figure(2)
% % scatter(time,avg_co2,'r','.')
% % hold on
% % scatter(time,cor_co2,'b','.')
% % fitlm(avg_co2,cor_co2)
% % figure(3)
% % t=time';
% % plot(t,resid)
%%

%to convert latitude and longitude from dec minutes to decimal degrees

% removing unwanted CO2

% cor_co2(H)=[];
% LAT(H)=[];
% LON(H)=[];
% avg_co2(H)=[];
% slope_intercept(H,:)=[];
% Dn(H,:)=[];

%dlmwrite('correctedtest.csv', all, 'delimiter', ',', 'precision', 15);
%to convert latitude and longitude from dec minutes to decimal degrees


% removing unwanted CO2
%H=find(cor_co2>800);
% cor_co2(H)=[];
% LAT(H)=[];
% LON(H)=[];
% avg_co2(H)=[];
% slope_intercept(H,:)=[];
% Dn(H,:)=[];
Valve1=ones(1,length(cor_co2))';
Valve2=ones(1,size(cor_atmos,1)).*2;
pCO2 = [Valve1 cor_co2];
pATM = [Valve2' cor_atmos];
Combined=[pCO2;pATM];
all=[Dn Combined avg_pH avg_temp avg_sal avg_lat avg_lon];
dlmwrite('calibrated_ST.csv', all, 'delimiter', ',', 'precision', 15); 