function c_table=OBScals_etm()
% calculate NTU curve for OBS's in a deployment


%%% 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 = [50 100 290 587];
    % 2181 0 at -0.030V
    OBS=[2181 0.345 0.715 2.212 4.368 2007 8 6; OBS];
    OBS=[2277 0.357 0.726 2.211 4.338 2007 8 6; OBS];
    OBS=[2278 0.425 0.843 2.519 4.946 2007 8 6; OBS];
    OBS=[2182 0.392 0.798 2.412 4.753 2007 8 6; OBS];
    OBS=[2183 0.331 0.763 2.505 5.007 2007 8 6; OBS];
    % 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'))