function hplstat(fname,HAL)
%*     Copyright c 1998 The board of trustees of the Leland Stanford     *
%*                      Junior University. All rights reserved.          *
%*     This script file may be distributed and used freely, provided     *
%*     this copyright notice is always kept with it.                     *
%*                                                                       *
%*     Questions and comments should be directed to Todd Walter at:      *
%*     walter@relgyro.stanford.edu                                       *
%HPLSTAT
%   HPLSTAT(FNAME, HAL) reads the given file of TMS Horizontal 
%   Protection Limit statistics and plots the color coded histogram.  
%   In addition the regions of integrity failures and unavailability
%   are shown.  The optional HAL argument is given in meters (default 30m).
%   Example code for generating the file is included at the bottom of hplstat.m

%   Todd Walter May 4th 1998
%   7 July 1998 AJHansen, added a HAL argument
%   14 July 1998 Todd Walter added variable ranges

%NOTE: To plot on a non-color printer type:
% colormap('gray')

if nargin < 1
    error('Must input a filename as a string')
end
if nargin < 2
    HAL = 30;
end
if ~isstr(fname) 
    error('Input must be a string variable')
end

c = 1;
colors = 'brygcmw';
clf

%load in data
fid = fopen(fname,'r');
err_bin=fread(fid,[100 1],'double');
sig_bin=fread(fid,[100 1],'double');
data=fread(fid,[100 100],'uint32');
data=data';
fclose(fid);

%determine the number of points and axis ranges
n_pts=sum(sum(data));

d_err_bin = mean(diff(err_bin));
x_lo_bnd = min(err_bin) - d_err_bin/2;
x_up_bnd = max(err_bin) + d_err_bin/2;

d_sig_bin = mean(diff(sig_bin));
y_lo_bnd = min(sig_bin) - d_sig_bin/2;
y_up_bnd = max(sig_bin) + d_sig_bin/2;

z_lo_bnd = 1;
z_up_bnd = max(max(data));

%clear plot
clf;

%plot each data point as a pixel
[i,j]=find(data);
face_mat=[[1 2 6 5]' [2 3 7 6]' [3 4 8 7]' [4 1 5 8]' [1 2 3 4]' [5 6 7 8]']';
colors=colormap;
for idx=1:length(i)
  z=log10(data(i(idx),j(idx)));
  vtx_mat=[err_bin(j(idx))+[-0.5 0.5 0.5 -0.5 -0.5 0.5 0.5 -0.5]'*d_err_bin sig_bin(i(idx))+[-0.5 -0.5 0.5 0.5 -0.5 -0.5 0.5 0.5]'*d_sig_bin [0 0 0 0 z z z z]'];
  c_idx = ceil(63*(log10(data(i(idx),j(idx)))/log10(z_up_bnd))) + 1;
  patch('Vertices', vtx_mat, 'Faces', face_mat, 'FaceColor',colors(c_idx,:), 'EdgeColor','none');
end

%determine availability and # of integrity failures

i_fail1=find(err_bin(j) >= HAL & sig_bin(i) < HAL);
n_fail1=sum(sum(diag(data(i(i_fail1),j(i_fail1)))));
i_fail2=find(err_bin(j)./sig_bin(i) >=1.0 & err_bin(j) < HAL);
n_fail2=sum(sum(diag(data(i(i_fail2),j(i_fail2)))));
i_fail3=find(err_bin(j)./sig_bin(i) >=1.0 & sig_bin(i) >= HAL);
n_fail3=sum(sum(diag(data(i(i_fail3),j(i_fail3)))));
i_cont=find(sig_bin(i) >= HAL);
n_cont=sum(sum(diag(data(i(i_cont),j(i_cont)))));
i_avail = find(err_bin(j) < HAL & sig_bin(i) < HAL);
n_avail=sum(sum(diag(data(i(i_avail),j(i_avail)))));


%set the axes limits and color values
set(gca,'XLim',[x_lo_bnd x_up_bnd]);
set(gca,'YLim',[y_lo_bnd y_up_bnd]);
set(gca,'CLim', [z_lo_bnd z_up_bnd]);


%show the region of normal operation
text(.37*(HAL - x_lo_bnd) + x_lo_bnd, .93*(HAL - y_lo_bnd) + y_lo_bnd, 0.75*z_lo_bnd, '{\fontsize{14pt}Normal Operation}');
if(n_avail/n_pts >= .999995)
  HT=text(.33*(HAL - x_lo_bnd) + x_lo_bnd, .86*(HAL - y_lo_bnd) + y_lo_bnd, 0.75*z_lo_bnd, '> 99.999%');
else
  HT=text(.37*(HAL - x_lo_bnd) + x_lo_bnd, .86*(HAL - y_lo_bnd) + y_lo_bnd, 0.75*z_lo_bnd, [num2str(100.0*n_avail/n_pts,'%6.3f'), '%']);
end
set(HT, 'FontSize', 14);


%outline the region of integrity failures
patch([HAL HAL x_up_bnd x_up_bnd],[y_lo_bnd HAL HAL y_lo_bnd],-[0.5 0.5 0.5 0.5],[1 0.1 0.1]);
HT=text(0.5*(x_up_bnd - HAL) + HAL,  .55*(HAL - y_lo_bnd) + y_lo_bnd, 0.75*z_lo_bnd, ['Integrity']);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
HT=text(0.5*(x_up_bnd - HAL) + HAL,  .45*(HAL - y_lo_bnd) + y_lo_bnd, 0.75*z_lo_bnd, ['Failures: ', int2str(n_fail1)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

%outline the region of HPL failures
patch([x_lo_bnd HAL HAL],[y_lo_bnd HAL y_lo_bnd],-[0.5 0.5 0.5],[1 0.55 0.55]);
HT=text(.67*(HAL - x_lo_bnd) + x_lo_bnd, .35*(HAL - y_lo_bnd) + y_lo_bnd, 0.75*z_lo_bnd, ['{\fontsize{14pt}HPL_{WAAS}}']);
set(HT,'HorizontalAlignment','Center');
HT=text(.67*(HAL - x_lo_bnd) + x_lo_bnd, .25*(HAL - y_lo_bnd) + y_lo_bnd, 0.75*z_lo_bnd, ['Failures: ', int2str(n_fail2)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

%outline the region of unavailability
patch([x_lo_bnd x_up_bnd x_up_bnd x_lo_bnd],[HAL HAL y_up_bnd y_up_bnd],-[0.5 0.5 0.5 0.5],[1 1 0.5]);
HT=text(.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, .70*(y_up_bnd - HAL) + HAL, 0.75*z_lo_bnd, '{\fontsize{14pt}System Unavailable}');
set(HT,'HorizontalAlignment','Center');
HT=text(.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, .55*(y_up_bnd - HAL) + HAL, 0.75*z_lo_bnd, ['Alarm Epochs: ', int2str(n_cont)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

%outline the region where integrity failures and unavailability overlap
patch([HAL x_up_bnd x_up_bnd],[HAL y_up_bnd HAL],z_lo_bnd*-[0.45 0.45 0.45],[1 .55 0.3]);
HT=text(.70*(x_up_bnd - HAL) + HAL, .325*(y_up_bnd - HAL) + HAL, 0.75*z_lo_bnd, ['{\fontsize{14pt}HPL_{WAAS}}']);
set(HT,'HorizontalAlignment','Center');
HT=text(.70*(x_up_bnd - HAL) + HAL, .175*(y_up_bnd - HAL) + HAL, 0.75*z_lo_bnd, ['Failures: ', int2str(n_fail3)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

hold on;
grid off;
%make the grid visible over the patch regions (integrity, unavailability)
ytick=get(gca,'YTick');
xtick=get(gca,'XTick');
nytick=length(ytick);
nxtick=length(xtick);
for i=1:nytick
    plot3([x_lo_bnd x_up_bnd],ytick(i)*[1 1], [0.6 0.6], 'k:');
end
for i=1:nxtick
    plot3(xtick(i)*[1 1], [y_lo_bnd y_up_bnd], [0.6 0.6], 'k:');
end

%label the axes and add a title
xlabel('{\fontsize{18pt}Error (m)}');
set(get(gca,'XLabel'),'Position',[.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, -.05*(y_up_bnd - y_lo_bnd) + y_lo_bnd z_lo_bnd]);
ylabel('{\fontsize{18pt}HPL_{WAAS} (m)}');
set(get(gca,'YLabel'),'Position',[-.06*(x_up_bnd - x_lo_bnd) + x_lo_bnd, .5*(y_up_bnd - y_lo_bnd) + y_lo_bnd z_lo_bnd]);
title (['Horizontal Performance (',int2str(n_pts),' epochs)']);
set(get(gca,'Title'),'FontSize',14);
set(get(gca,'Title'),'Position',[.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, 1.02*(y_up_bnd - y_lo_bnd) + y_lo_bnd z_lo_bnd]);


%put the sigma_H_major scale on the right hand y-axis
K_h_pa = 6.18;
for i =ceil(y_lo_bnd/K_h_pa):floor(y_up_bnd/K_h_pa)
  if(abs(i*K_h_pa - (0.5*(y_up_bnd - y_lo_bnd) + y_lo_bnd)) > 0.05*(y_up_bnd - y_lo_bnd))
    plot3([.99*(x_up_bnd - x_lo_bnd) + x_lo_bnd x_up_bnd], i*K_h_pa*[1 1],[0.65 0.65],'k');
    text(1.02*(x_up_bnd - x_lo_bnd) + x_lo_bnd, i*K_h_pa, 0.65, int2str(i));
  end
end
HT=text(1.04*(x_up_bnd - x_lo_bnd) + x_lo_bnd,0.5*(y_up_bnd - y_lo_bnd) + y_lo_bnd, 0.65,'{\fontsize{14pt}\sigma}_{H_{major}} (m)');
set(HT,'HorizontalAlignment','Center');
set(HT,'Rotation',90);

%put the color scale up on the right hand side
H=colorbar('vert');
set(get(H,'Ylabel'),'String','{\fontsize{14pt}Number of Points per Pixel}');
set(H,'YScale','log');
set(H,'YLim',[z_lo_bnd z_up_bnd]);
set(H,'CLim', [z_lo_bnd z_up_bnd]);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%EXEMPLAR C CODE for generating hpl_????.hst files

%Initialization:
%#define HIST_BINS 100
%unsigned long hpl_stat[HIST_BINS][HIST_BINS];
%memset(hpl_stat, 0, HIST_BINS*HIST_BINS*sizeof(unsigned long));

%Accumulation:
%
%int j,k;
%double N_err, E_err;
%double N_cov, E_cov;
%double  NEx_cov, sig_H_major;
%
%  j = (int)floor(2.0*sqrt(N_err*N_err + E_err*E_err);
%  if(j >= HIST_BINS)
%      j = HIST_BINS - 1;
%  else if(j < 0)
%      j = 0;
% /* sig_H_major is the major axis of the horizontal error ellipse. 
%  * It is defined in Appendix J of the WAAS MOPS */
% sig_H_major =  sqrt(0.5*(E_cov + N_cov) + 
%                       sqrt(0.25*(E_cov - N_cov)*(E_cov - N_cov) + 
%                       NEx_cov*NEx_cov));
%  k = (int)floor(2.0*6.18*sig_H_major);
%  if(k >= HIST_BINS)
%      k = HIST_BINS - 1;
%  else if(k < 0)
%      k = 0;
%  hpl_stat[k][j]++;

%Writing to file:
%char filename[256];
%int i,ID=0x????;
%double herr_bin_value[HIST_BINS];
%double hpl_bin_value[HIST_BINS];
%FILE *fp;
%
% for(i=0; i<HIST_BINS;i++){
%    verr_bin_value[i] = ((double)(i) + 0.5)/2.0;
%    vpl_bin_value[i]  = ((double)(i) + 0.5)/2.0;
% }
% sprintf(file_name,"hpl_%04X.hst",ID);
% fp = fopen(file_name,"wb");
% fwrite(herr_bin_value, sizeof(double), HIST_BINS,fp);
% fwrite(hpl_bin_value, sizeof(double), HIST_BINS,fp);
% fwrite(hpl_stat, sizeof(unsigned long), HIST_BINS*HIST_BINS, fp);
% fclose(fp);