%PAROPLOTHYDRA check paros operating range against Hydra limits
% function [x, y] = paroplothydra(T, serial, depth, hydraoffset);
%	generate a figure which shows the operating ranges of
%	various specific Paroscientific pressure sensors
%   T = temperature
%	serial = sensor serial number to get data, as a char string
%           if serial is provided, only that sensor is plotted
%   depth = deployment depth, db (m)
%   hydraoffset = the offset setting in the Hydra, will be plotted
%   [x, y] = the actual data used to generate the graph... for the sensor requested
%	
% Note that the sensor's calibration information must be entered or updated
% in parpnormcdt0.m

% Written by Marinna Martini
% U.S. Geological Survey Woods Hole Science Center

function [x, y] = paroplothydra(T, serial, depth, hydraoffset)
prog_ver = 'SVN $Revision: 1082 $';
disp(sprintf('%s version %s',mfilename,prog_ver))

if exist('T','var')~=1, T = 20;	end % deg C
if exist('depth','var')~=1, depth=650; end
if exist('hydraoffset','var')~=1, 
    hydraoffset=21; 
    disp(sprintf('Hydraoffset is assumed to be %d',hydraoffset))
else
    disp(sprintf('Hydraoffset given by user %d',hydraoffset))
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}=750000:100:790000; c{i} = 'g';  % 130 m sensor
i=i+1; serials{i} = {'1274'};  s{i}=750000:100:790000; c{i} = 'r';  % 130 m sensor
i=i+1; serials{i} = {'47037'}; s{i}=770000:100:950000; c{i} = 'c'; % 20 m sensor
i=i+1; serials{i} = {'60005'}; s{i}=840000:100:950000; c{i} = 'm'; % 60 m sensor?
i=i+1; serials{i} = {'60006'}; s{i}=840000:100:960000; c{i} = 'b'; % 60 m sensor
i=i+1; serials{i} = {'69367'}; s{i}=840000:100:970000; c{i} = 'g'; % 20 m sensor
i=i+1; serials{i} = {'68022'}; s{i}=840000:100:990000; c{i} = 'r'; % 130 m sensor
i=i+1; serials{i} = {'68093'}; s{i}=840000:100:950000; c{i} = 'c'; % 270 m sensor
i=i+1; serials{i} = {'70390'}; s{i}=840000:100:950000; c{i} = 'm'; % 60 m sensor
i=i+1; serials{i} = {'571'};   s{i}=750000:100:850000; c{i} = 'b'; % 130 m sensor
i=i+1; serials{i} = {'1045'};  s{i}=750000:100:790000; c{i} = 'g';  % 130 m sensor
i=i+1; serials{i} = {'1557'};  s{i}=750000:100:790000; c{i} = 'r';  % 130 m sensor
i=i+1; serials{i} = {'45417'}; s{i}=750000:100:790000; c{i} = 'c'; % 20 m sensor
i=i+1; serials{i} = {'6763'};  s{i}=750000:100:1000000; c{i} = 'm'; % 600 m sensor
i=i+1; serials{i} = {'1558'};  s{i}=750000:100:800000; c{i} = 'b'; % 130 m sensor
i=i+1; serials{i} = {'62493'}; s{i}=840000:100:960000; c{i} = 'g'; % 60 m sensor
i=i+1; serials{i} = {'62495'}; s{i}=840000:100:970000; c{i} = 'r'; % 60 m sensor
i=i+1; serials{i} = {'69506'}; s{i}=840000:100:980000; c{i} = 'c'; % 130 m sensor
i=i+1; serials{i} = {'454171'}; s{i}=750000:100:800000; c{i} = 'm'; % 20 m sensor
i=i+1; serials{i} = {'38196'}; s{i}=750000:100:810000; c{i} = 'b';  % 130 m sensor
i=i+1; serials{i} = {'38147'}; s{i}=750000:100:820000; c{i} = 'g';  % 130 m sensor

if exist('serial','var'),
    for idx = 1:length(serials),
        if strcmp(serial,char(serials{idx})), 
            data = s{idx}(:);
            serials = cell(1);
            s = cell(1);
            c = cell(1);
            break; 
        end
    end
    s{1} = data;
    serials{1} = serial;
    c{1} = 'b';
end
%offsets = zeros(length(serials),1);

% translate to frequency
N=6208;        % number of pressure pulses counted, typical
Fcl=2457600;    % clock frequency, typical Tattletale
for i = 1:length(serials)
    freq=(2*N*Fcl)./s{i}(:); % freq = sensor output frequency, hz
    s{i} = freq;
end

p = cell(length(serials));
sens = cell(length(serials));
for i = 1:length(serials),
    p{i} = paronormcdt0(s{i}, T, serials{i});
    sens{i}=(p{i}(length(p{i}))-p{i}(1))/(s{i}(length(p{i}))-s{i}(1));
    disp(sprintf('Paros #%s has sensitivity of %g db/Hz',...
        char(serials{i}(:)),sens{i}(:)));
%     if exist('hydraoffset','var'),
%         hydralim{i}(1) = paronormcdt0(s{i}, T, serials{i});
%     end
end

rot=270;		% text rotation
% x axis is in psia
%SD = (depth./10).*14.5038;
%surface = 14.7;
% x axis is in db
%SD = depth;
surface = 10.135;

% check the offset
surfidx = find(p{i}>surface,1,'first');
freqatsurf = s{i}(surfidx);
if freqatsurf-hydraoffset*1000 < 16000,
    disp(sprintf('Hydraoffset of %d is OK',hydraoffset))
else
    newoffset = ceil((freqatsurf-16000)/1000);
    disp(sprintf('recommend changing hydraoffset from %d to %d or more',hydraoffset,newoffset))
end

clf
hold on
for i=1:length(serials),
    plot(p{i},s{i},c{i})
end
if exist('hydraoffset','var'),
    xl = get(gca,'xlim');
    plot(xl, [freqatsurf-16000 freqatsurf-16000],'k-')
    tfreq=text(diff(xl)/2,freqatsurf-16000,...
        sprintf('Hydra frequency measurement limit for pressfreqoffset = %d',hydraoffset));
    set(tfreq,'vert','baseline','horizontalalignment','center')
end
yl=get(gca,'ylim');
plot([surface surface],yl,'k-')
xl=get(gca,'xlim');
set(gca,'xlim',[0 xl(2)]);
ylabel('Output, Hz')

% fix x axis display
xlabel('Input pressure, db (absolute)')

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

yl=get(gca,'ylim');
if exist('depth','var'),
    plot([depth+surface depth+surface],yl,'k-')
    tbottom=text(depth+surface,yl(2),'Deployed Depth');
    set(tbottom,'rot',rot,'color','k','vert','baseline')
end
tsurface=text(surface,yl(2),'  Sea Surface');
set(tsurface,'rot',rot,'color','k','vert','baseline')
if exist('serial','var'),
    title(sprintf('Paros Sensor %s sensitivity %f db/Hz',serials{1}(:),sens{1}))
    %xl=get(gca,'xlim');
    %set(gca,'xlim',[0 depth+0.1*depth]);
else
    title('USGS Woods Hole Paros Sensor Behavior')
end
hold off

if exist('serial','var'),
    x = p{1}(:);
    y = s{1}(:);
else
	x=[];
	y=[];
end