function c_table=OBScals_hatteras09()
% calculate NTU curve for OBS's in a deployment, from lab calibrations 
% these were from 12/12/08

%%% START USGS BOILERPLATE -------------% There is no published documentation for these programs- use the 
%  internal comments and help function.
% Programs written in Matlab- as early as 7.2.0.232 (R2006a) and 
% most recent modifications in  v7.8.0 (R2009a)
% Programs ran on PC with Windows XP Professional OS, and RHEL4 linux.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%%% END USGS BOILERPLATE --------------

 
clear OBS NTU
close all
OBS = [];

if 1,
    NTU = [0 50 100 333 533 588 660];
    % order: sn cal_vals and date (yyyy_mm_dd)
    OBS=[2277 -0.1474 0.853 1.852 6.348 9.88 10.959 11.1142 2008 12 12];
    OBS=[2184 -0.0865 0.957 2.003 6.663 10.516 11.0741 11.0768 2008 12 12; OBS];
    OBS=[2181 -0.198 0.78 1.744 6.098 9.598 10.54 11.0878 2008 12 12; OBS];
    OBS=[2183 -0.3856 0.804 1.973 7.273 11.0896 11.0896 11.0896 2008 12 12; OBS];
    OBS=[2278 -0.1455 0.989 2.095 7.092 11.0931 11.1083 11.1083 2008 12 12; OBS];
    OBS=[2182 -0.0954 1.034 2.103 7.152 11.0656 11.0742 11.0742 2008 12 12; OBS];
    OBSl=[7409 0.0019 0.039 0.0747 0.2372 0.3602 0.3971 0.4467 2008 12 12];
    OBS3h=[7409 0.00330 1525 0.2962 0.9405 1.444 1.5871 1.756 2008 12 12];
    % 2184 0 at 0.033V
    % OBS=[2184 0.435 0.833 2.421 4.732 2007 8 6; OBS];
end

% subtracting the zero offset should end up with a zero intercept, but in
% the two cases I'm testing, I still end up with something like -1.6 for
% both 2181 and 2184
%zoffs= [0.033; 0; 0; 0; 0; -0.030];
%for ik=1:length(zoffs)
%    OBS(ik,2:end)=OBS(ik,2:end)-zoffs(ik);
%end

if 0,
    % D&A two point
    NTU = [0.3 2000];
    OBS = [2181 0.003 5.002 2003 5 7; OBS];
end

[nOBS, ncol] = size(OBS);

dstring = datestr(now,1);

for iOBS = 1:nOBS,    
    coeff(iOBS,:) = polyfit(OBS(iOBS,2:end-3),NTU,1);
    fit(iOBS,:) = polyval(coeff(iOBS,:),OBS(iOBS,2:end-3));

    figure
    subplot(2,1,1)
    plot(OBS(iOBS,2:end-3),NTU,'-',OBS(iOBS,2:end-3),fit(iOBS,:),'x')
    ylabel('NTU')
    xlabel('OBS output VDC')
    title(sprintf('Two degree fit for OBS #%d',OBS(iOBS,1)));
    text('units','norm','pos',[0.05 0.9],'string',dstring)
    text('units','norm','pos',[0.05 0.8],'string',...
        sprintf('Calibration date %s',datestr(datenum(OBS(iOBS,end-2:end)))))
    subplot(2,1,2)
    plot(OBS(iOBS,2:end-3),NTU-fit(iOBS,:),'x')
    xlabel('NTU')
    ylabel('Difference - NTU')
    text('units','norm','pos',[0.05 0.95],'string',...
        sprintf('STD of difference = %4.3f',std(NTU-fit(iOBS,:))))
    title(sprintf('NTU = %6.4fV + %6.4f',coeff(iOBS,:)'))
end

%save the [OBS ID, slope and intercept] from polyfit
c_table=[OBS(:,1) coeff];
% format the output nicely, if you want.
% sprintf('%i %d %d\n',fliplr(c_table'))

% output NTU values for this cal [sn slope intercept]
%2277 5.550809e+01 -5.705454e-01 : NTU = 55.5081V + -0.5705
%2184 5.467069e+01 -6.184162e+00 : NTU = 54.6707V + -6.1842
%2181 5.652432e+01 3.260284e+00  : NTU = 56.5243V + 3.2603
%2183 5.223164e+01 3.075535e+00  : NTU = 52.2316V + 3.0755
%2278 5.330388e+01 -6.600147e+00 : NTU = 53.3039V + -6.6001
%2182 5.362673e+01 -9.115388e+00 : NTU = 53.6267V + -9.1154