%PAROPLOT.M check paros operating range against TT6 limits
% function [x, y] = paroplot1(type, N, Fcl, T, Timeout, serial, depth);
%	generate a figure which shows the operating ranges of
%	various specific Paroscientific pressure sensors
%	type = counts | freq | freqdb | countdb
%		counts = psia vs counts
%		freq = psia vs freq
%		freqdb = db vs freq
%		countdb = db vs counts
%	N = period count
%   Fcl = clock frequency
%   T = temperature
%	Timeout = TTBASIC timeout setting
%	serial = sensor serial number to get data, as a char string
%   depth = deployment depth, db (m)
%   [x, y] = the actual data used to generate the graph... for the sensor requested
%	
function [x, y] = paroplot(type, N, Fcl, T, Timeout, serial, depth)

if ~exist('N','var'),	
	N=6208;	
end
if exist('Timeout')~=1,
	Timeout=43;		% setting used in the TTBASIC period statement
end
if exist('Fcl')~=1, Fcl=2457600;	end % Hz
if exist('T')~=1, T = 20;	end % deg C
G=0.68946498;     % conversion constant, should never change
                  % psi to decibars (approx meters)
if exist('type')~=1, type='freqdb'; end
if exist('depth')~=1, depth=650; end

% comment out a sensor you don't want displayed
i=0;
i=i+1; serials{i} = {'52103'}; s{i} = 840000:100:950000; c{i} = 'b';  % 130 m sensor
i=i+1; serials{i} = {'1703'}; s{i}=755000:100:850000; c{i} = 'g';  % 130 m sensor
i=i+1; serials{i} = {'1274'}; s{i}=730000:100:800000; c{i} = 'r';  % 
i=i+1; serials{i} = {'47037'}; s{i}=770000:100:830000; c{i} = 'c'; % 20 m sensor
i=i+1; serials{i} = {'60005'}; s{i}=840000:100:930000; c{i} = 'm'; % 60 m sensor?
i=i+1; serials{i} = {'60006'}; s{i}=840000:100:940000; c{i} = 'b'; % 60 m sensor
i=i+1; serials{i} = {'69367'}; s{i}=850000:100:920000; c{i} = 'g'; % 20 m sensor
i=i+1; serials{i} = {'68022'}; s{i}=840000:100:940000; c{i} = 'r'; % 130 m sensor
i=i+1; serials{i} = {'68093'}; s{i}=850000:100:950000; c{i} = 'c'; % 270 m sensor
i=i+1; serials{i} = {'70390'}; s{i}=840000:100:920000; c{i} = 'm'; % 60 m sensor
i=i+1; serials{i} = {'571'}; s{i}=770000:100:870000; c{i} = 'b'; % 130 m sensor
i=i+1; serials{i} = {'1045'}; s{i}=740000:100:820000; c{i} = 'g'; % 130 m sensor
i=i+1; serials{i} = {'1557'}; s{i}=740000:100:825000; c{i} = 'r'; % 130 m sensor
i=i+1; serials{i} = {'45417'}; s{i}=770000:100:840000; c{i} = 'c'; % 20 m sensor
i=i+1; serials{i} = {'6763'}; s{i}=800000:100:910000; c{i} = 'm'; % 600 m sensor
i=i+1; serials{i} = {'1558'}; s{i}=760000:100:840000; c{i} = 'b'; % 130 m sensor
i=i+1; serials{i} = {'62493'}; s{i}=840000:100:910000; c{i} = 'g'; % 60 m sensor
i=i+1; serials{i} = {'62495'}; s{i}=840000:100:950000; c{i} = 'r'; % 60 m sensor
i=i+1; serials{i} = {'69506'}; s{i}=840000:100:940000; c{i} = 'c'; % 130 m sensor
i=i+1; serials{i} = {'454171'}; s{i}=770000:100:840000; c{i} = 'm'; % 20 m sensor
i=i+1; serials{i} = {'38196'}; s{i}=730000:100:800000; c{i} = 'b'; % 270 m sensor
offsets = zeros(length(serials),1);

disp(sprintf('Paros #, Sensitivity & accuracy'));
for i = 1:length(serials),
    p{i} = pnorm(s{i},N,Fcl,T,serials{i});
    acc{i} = (0.01/100)*(p{i}(end)-p{i}(1));
    disp(sprintf('%s  accuracy %5.3f db = %5.1f mm over %5.0f db range',char(serials{i}(:)),acc{i}(:),acc{i}(:)*100,p{i}(end)-p{i}(1)));
    sens{i}=(p{i}(length(p{i}))-p{i}(1))/(s{i}(length(p{i}))-s{i}(1));
    disp(sprintf('\t  sensitivity %g db/count for period count',sens{i}(:)));
    %disp(sprintf('\t  output signal range: %d to %d Hz',s{i}(1),s{i}(end)));
    % take the same pressure range and distribute over 0-5v
    vres{i} = 5/(p{i}(end)-p{i}(1)); % V/db
    disp(sprintf('\t  sensitivity: %g V/db for analog conversion',vres{i}));
    % for 16 bit A/D 2^15 counts per 5 volts
    % 2^15 counts/5 v 
    disp(sprintf('\t  sensitivity: %g db/count for 16 bit (2^15) analog conversion',(5/(2^15))*(1/vres{i})));
    disp(sprintf('\t  freq range: %g to %g',(2*N*Fcl)./s{i}(1),(2*N*Fcl)./s{i}(end))); % freq = sensor output frequency, hz
end

% Find the MIDAS upper limit due to max input frequency in counts
% TT6 with fast xtal, input freq to period circuits 
% must be < 20 KHz
% MIDAS currently divides frequencies by 2 to stay below this limit
% therefore: Max Paros Freq = 20 * 2 = 40 KHz
MCU = (2*N*Fcl)/40000;	% = 762839 counts

% Find the MIDAS lower limit due to the period command's time out feature
% Typical timeout setting is Timeout
% MIDAS currently divides frequencies by 2 
% 2*N/(timeout in sec) = lower sensor limit in Hz
MCL = 2*N/((Timeout-6)*10/1000);	
MCL = (2*N*Fcl)/MCL;	% = 909305 counts

% fix y axis labels
if strcmp(type,'freq') || strcmp(type,'freqdb'),
	MCU=(2*N*Fcl)/MCU;	% convert to Hz
	MCL=(2*N*Fcl)/MCL;	% convert to Hz
	rot=270;		% text rotation
	% convert counts to Hz
    for i=1:length(serials),
        s{i} = (N*Fcl*2)./s{i};
    end
elseif strcmp(type,'counts') || strcmp(type,'countdb'),
	rot=-90;		% text rotation
end

% fix x axis display
if strcmp(type,'freq') || strcmp(type,'counts'),
	% x axis is in psia
	SD = (depth./10).*14.5038;
	surface = 14.7;
	% pnorm returns decibars, so convert to psia
    for i=1:length(serials),
        p{i} = p{i}./G;
    end
elseif strcmp(type,'freqdb') || strcmp(type,'countdb'),
	% x axis is in decibars
	SD = depth;
	surface = 10.135;
end


clf
hold on
for i=1:length(serials),
    plot(p{i},s{i},c{i})
end
yl=get(gca,'ylim');
plot([SD SD],yl,'k-')
plot([surface surface],yl,'k-')
xl=get(gca,'xlim');
plot(xl,[MCU MCU],'k-')
plot(xl,[MCL MCL],'k-')
if type == 'freqdb', set(gca,'xlim',[0 xl(2)]); end
% fix y axis labels
if strcmp(type,'freq') | strcmp(type,'freqdb'),
	ylabel('Output, Hz')
else,
	ylabel('Output, count')
end

% fix x axis display
if strcmp(type,'freq') | strcmp(type,'counts'),
	xlabel('Input pressure, psia')
else,
	xlabel('Input pressure, db (absolute)')
end

for i = 1:length(serials),
    tptr = text(p{i}(length(p{i})),s{i}(length(s{i})),char(serials{i}(:)));
    set(tptr,'horiz','left','vertic','baseline','color',c{i})
end

tmidas=text((xl(2)-xl(1))/2,MCU,'Max period input limit');
set(tmidas,'color','k','vert','baseline')
tmidas=text((xl(2)-xl(1))/2,MCL,'Period timeout limit');
set(tmidas,'color','k','vert','baseline')

tbottom=text(SD,MCU,'Deployed Depth');
set(tbottom,'rot',rot,'color','k','vert','baseline')
tsurface=text(surface,MCU,'Sea Surface');
set(tsurface,'rot',rot,'color','k','vert','baseline')
title('USGS Woods Hole Paros Sensor Behavior')

hold off

if exist('serial')==1,
	% user wants the data on a sensor
    for i=1:length(serials),
        sidx(i) = strcmp(char(serials{i}(:)),serial);
    end
    sidx = find(sidx);
    x = p{sidx};
    y = s{sidx};
else
	x=[];
	y=[];
end