%PARONORMCDT0 normalize Paroscientific frequency output for Paros cals with C,D,T & 0 constants
%
%       press = paronormcdt0(freq, temp, sensor)
%       press = normalized pressure in db
%       freq = sensor frequency output in Hz
%       temp = temperature,
%               if a scalar, it is applied to all pressure points
%               if a vector, each pressure point is temperature compensated
%                       for each corresponding temperature
%		sensor = sensor serial number as a number
%		NOTES:  
%           9.75 db = approx. 14.7 psia at sea level
%           The sensor's calibration information must be entered or updated
%           in this m-file

%       written by Marinna Martini & Fran Hotchkiss
%               U.S. Geological Survey, Woods Hole, MA
%       8/25/03 version 0.0 based on paronorm

function press=paronormcdt0(freq, temp, sensor)
prog_ver = 'SVN $Revision: 1082 $';
disp(sprintf('%s version %s',mfilename,prog_ver))

G=0.68946498;     % conversion constant, should never change
                  % psi to decibars (approx meters)

% Calibration constants for temperature compensated calibrations
%	using the CDT0 equation
% ser3=[sensor C1 C2 C3 D1 D2 T1 T2 T3 T4 T5 C D T0 Date(y m d)]
ser3=[47037 175.1086 -5.70706E-3 -1.87896E-5 0.045638 0.0 24.37787 -1.35891E-4 1.19770E-6 1.23228E-9 0.0 174.9805 0.045638 24.37555 1991 10 28];
ser3=[52103 981.2742 -1.76254E-3 -9.66989E-5 0.026698 0.0 27.68403 -6.97677E-5 9.30850E-7 1.10183E-9 0.0 981.1946 0.026698 27.68299 1992 10 26; ser3];
ser3=[60005 590.4284 -4.05709E-4 -6.89342E-5 0.020903 0.0 27.83571 -1.12340E-4 1.02844E-6 1.30711E-9 0.0 NaN      NaN      NaN      1995 01 16; ser3];
ser3=[60006 587.4469 -1.77288E-3 -7.31431E-5 0.021239 0.0 27.92124 -7.23719E-5 9.83991E-7 1.45109E-9 0.0 NaN      NaN      NaN      1995 01 16; ser3];
ser3=[62493 624.4604 -2.27217e-3 -6.80591e-5 0.029671 0.0 27.70230 -6.21501e-5 1.01351e-6 1.32810e-9 0.0 NaN      NaN      NaN      1995 10 02; ser3];
ser3=[62495 589.1866  4.24177E-4 -6.25844E-5 0.026538 0.0 27.78682 -4.32422E-5 9.63356E-7 1.41121E-9 0.0 NaN      NaN      NaN      1995 10 02; ser3];
%ser3=[68022 1024.946 -3.76508E-3 -1.12975E-4 0.040873 0.0 27.74432 -5.64250E-5 9.08633E-7 1.68019E-9 0.0 NaN      NaN      NaN      1997 10 13; ser3];
ser3=[68022 1024.995 -3.76508E-3 -1.12975E-4 0.040873 0.0 27.74370 -5.64250E-5 9.08633E-7 1.68019E-9 0.0 NaN      NaN      NaN      2004 08 10; ser3];
ser3=[68093 2046.522 -4.52377E-3 -2.06479E-4 0.021522 0.0 27.79633 -6.83120E-5 9.86269E-7 1.26877E-9 0.0 NaN      NaN      NaN      1998 01 21; ser3];
%ser3=[69367 230.5027 -4.20275e-3 -7.58202e-6 0.036831 0.0 27.65536 -8.77608e-5 7.73802e-7 2.25430e-9 0.0 NaN      NaN      NaN      2003 09 05; ser3];
ser3=[69367 230.5070 -4.20275e-3 -7.58202e-6 0.036831 0.0 27.65511 -8.77608e-5 7.73802e-7 2.25430e-9 0.0 NaN      NaN      NaN      2008 11 14; ser3];
ser3=[69506 991.7378 -8.68840E-3 -1.23004E-4 0.036005 0.0 27.77617 -6.44352E-5 9.71216E-7 1.34311E-9 0.0 NaN      NaN      NaN      1998 02 06; ser3];
ser3=[70390 606.2662  1.58605E-3 -5.55501E-5 0.024865 0.0 27.75208 -7.33351E-5 1.06256E-6 9.14717E-10 0.0 NaN     NaN      NaN      2007 08 14; ser3];
ser3=[38147 2101.8740 -4.44603e-03 -3.40722e-04 -1.16430e-02 0.0 2.466958e+01 -1.07382e-04 1.27366e-06 0.0  0.0 NaN  NaN NaN  1998 02 06; ser3];
ser3=[45417 167.6087 -5.70595e-03 -1.57493E-5 0.052763 0.0 2.448838e+01 -1.12994e-04 1.12587e-06 9.29819E-10  0.0 NaN  NaN NaN  2004 08 09; ser3];

% determine if sensor contains valid data
if exist('sensor','var') ~= 1,
	sensor = ser3(1,1);
    disp(sprintf('Using a demonstration sensor, %d',sensor))
end
if ischar(sensor),
   sensor=str2double(sensor);
elseif iscell(sensor),
   sensor=str2double(char(sensor));
end
if find(ser3(:,1)==sensor), 
	sptr=find(ser3(:,1)==sensor);
else disp(['Unrecognized Sensor: ' int2str(sensor) ' using default constants.']),
	sptr=1;
end

if exist('temp','var') ~= 1,
        temp = 21.0;      % if not provided, make it up!
end

C1=ser3(sptr,2); C2=ser3(sptr,3); C3=ser3(sptr,4); 
D1=ser3(sptr,5); D2=ser3(sptr,6); 
T1=ser3(sptr,7); T2=ser3(sptr,8); T3=ser3(sptr,9); T4=ser3(sptr,10); T5=ser3(sptr,11); 
C=ser3(sptr,12); D=ser3(sptr,13); T0=ser3(sptr,14); 
% adjust calibration coefficients for temperature dependence
if length(temp) == 1,  % if temperature is scalar
    C=C1+(C2.*temp)+(C3.*(temp.^2));	%	psia
    D=D1+(D2.*temp);	% no units
    T0=T1+(T2.*temp)+(T3.*(temp.^2))+(T4.*(temp.^3))+(T5.*(temp.^4));	% usec
    % Calculate pressure in psi
    %T0=T0./(1000000);	% convert to sec
    K = 1.0-((T0.*freq./1000000).^2);
    %period=(1./freq).*1000000;
    %K = 1.0-((T0.^2)./(period.^2));
    press = (C.*K).*(1.0 - (D.*K));
elseif (length(temp) > 1) && (length(temp) == length(count)),
    C=C1+C2.*temp+C3.*(temp.^2);	%	psia
    D=D1+D2.*temp;	% no units
    T0=T1+T2.*temp+T3.*(temp.^2)+T4.*(temp.^3)+T5.*(temp.^4);	% usec
    % Calculate pressure in psi
    %T0=T0./(1000000);	% convert to sec
    K = 1.0-((T0.*freq./1000000).^2);
    press = C.*K.*(1-D.*K);
else
    % Calculate pressure in psi
    %T0=T0./(1000000);	% convert to sec
    K = 1.0-((T0.*freq./1000000).^2);
    press = C.*K.*(1-D.*K);
end

% convert psi to decibars
press=press.*G;

% replace trash with NaN so it won't be plotted
flags=find(isfinite(press)==0);
if length(flags) == 1,
    press(flags)=NaN;
elseif isempty(flags) == 0,
    press(flags)=ones(size(flags)).*NaN;
end