function  BP9_Field_GAUGE(filename1,filename2)
% -------------------------------------------------------------------------
%               NTHMP (National Tsunami Hazard Mitigation Program)
% -------------------------------------------------------------------------
%                          2010 NTHMP VALIDATION WORKSHOP
%              Venue: Galveston, TAMU TX
%                     Margh 30- Apr. 1
%                               ---- o ----
%               MATLAB SCRIPT:  Okushiri Island
%               TYPE:           Field Benchmarking Problem (Field)
%                               NUMBER: 09
%               TSUNAMI SOURCE: Coseismic Deformation (CS)                       
%               SCRIPT NAME   : BP9_Field_GAUGE.m
%                
%               AUTHOR(S)     :J. Horrillo, Gyeong-Bo Kim.  	        
%                              Copyright 2011 NTHMP/NOAA
%                              Revision: 0
%                              Date: Sep/2011
%               REFERENCE:
%               Synolakis, C.E. (1986): The Runup of Long Waves. Ph.D.
%                   Thesis, California Institute of Technology, Pasadena,
%                   California, 91125, 228 pp. 
% 
%                   NOAA/PMEL Center for Tsunami research
%               http://nctr.pmel.noaa.gov/benchmark/index.html
% 
%                   Synolakis, C.E. (1987): The runup of solitary waves.
%               J. Fluid Mech., 185, 523-545.  
% 
%               INCLUDED FILES WITH THIS SCRIPT: 
%
%               Essashi_Iwanai.mat (Field data file for time sierie at 
%               gages 
%
%                 
% 
%               CONTACT INFO:   Dr. Juan J. Horrillo, 
%                               Texas A&M University at Galveston 
%                               Maritime Systems Engineering (MASE) 
%                               200 Seawolf Park or P.O. Box 1675 
%                               Galveston, TX  77553-1675 
%                               Phone: (409) 740-4465
%                               horrillj@tamug.edu
%                               Website: http://www.tamug.edu/
% 
% -------------------------------------------------------------------------
%   >  BP9_Field_GAUGE(filename1,filename2)
%
%   BP9_Field_GAUGE reads input files 'filename1' and 'filename2' containing 
%   model results and compares them with field data recorded at gages Iwanai
%   and Esashi.
%
%   filename1 and filename2 are model generated files corresponding to
%   to Iwanai and Esashi respectively.
%    
%   Input files 'filename1' and 'filename2' need to be named as follows:
%   XXXXXX_BP9_IWANAI.dat and XXXXXX_BP9_ESASHI.dat, where XXXXXX is
%   the name of the numerical model that is tested,
%   eg. NEOWAVE_BP9_IWANAI.dat and NEOWAVE_BP9_ESASHI.dat.
% 
%   The input files, 'filename1' and 'filename2' must have the following
%   configuration:
%
%   1- Files 'filename1' is a matrix ASCII files of size
%      121x2(row x column) as shown below (no heading).
%      The number of rows could be different of 689, but it is
%      desirable to sample at least every 0.50.     
%      t ranges= [0 60]min  with icrements of 0.5.
%
%   filename1
%   +---------+-------------+
%   |  t(min)         eta(m)   
%   +---------+-------------+ 
%     0.00000        0.00000 
%     0.50000        0.00414 
%     1.00000        0.00829 
%     1.50000        0.01243 
%     2.00000        0.01657 
%           .              . 
%           .              .  
%           .              . 
%    57.00000        0.35074 
%    57.50000        0.35626 
%    58.00000        0.36179 
%    58.50000        0.36731 
%    59.00000        0.37257 
%    59.50000        0.37533 
%    60.00000        0.37809 
%   +--------+--------------+
%
%   2- Files 'filename2' is a matrix ASCII files of size
%      601x2(row x column) as it is shown below (no heading).
%      The number of rows could be different of 401, but it is
%      desirable to sample at least every 0.10.     
%      t ranges= [0 60]min with icrements of 0.1. 
%
%   filename2
%   +---------+-------------+
%   |  t(min)         eta(m)   
%   +---------+-------------+ 
%     0.00000       -0.01690 
%     0.10000       -0.01690 
%     0.20000       -0.01690 
%     0.30000       -0.01690 
%     0.40000       -0.01690 
%     0.50000       -0.01690 
%           .              . 
%           .              .  
%           .              . 
%    59.60000       -0.13305 
%    59.70000       -0.15487 
%    59.80000       -0.17551 
%    59.90000       -0.19370 
%    60.00000       -0.21190 
%   +-------+---------+---------+
%
%   NOTE:
%   This script is self-containing. It loads the lab. data automatically to
%   facilitate user's comparison, peer review and determine model 
%   performance.
%   The field data have been saved on files Iwanai_Esashi.mat and 
%   Be sure that the lab data files Iwanai_Esashi.mat 
%   resides in the same directory or subdirectory as the Matlab script
%   (BP9_Field_GAUGE.m).
%   The user's input files 'filename1' and 'filename2' must be present
%   in the same directory or subdirectory as well.
%
%   This script uses the normalized root mean square deviation (NRMSD)
%   or root mean square error to measure the numerical model precision.
%   NRMSD measures the differences between values predicted by the
%   numerical model and the field gauge data.
% 
%   Examples:
%   > BP9_Field_GAUGE
%     In this example the Matlab script  plots field data at Iwanai and
%     Esashi.
%   > BP9_Field_GAUGE('NEOWAVE_BP9_IWANAI.dat','NEOWAVE_BP9_ESASHI.dat')
%     In this example the Matlab script plots both  user's model
%     result cases against field data, (Iwanai and Esashi).
%   > BP9_Field_GAUGE('NEOWAVE_BP9_IWANAI.dat',[])
%     In this example the Matlab script just plots user's model results
%     against field data for Iwanai case.
%   > BP9_Field_GAUGE([],'NEOWAVE_BP9_ESASHI.dat')
%     In this example the Matlab script just plots user's model results
%     against field data for Esashi case.
% -------------------------------------------------------------------------
help BP9_Field_GAUGE
str3=[' '];
str4=[' '];
err=nan;
err1=nan;
if nargin < 1
   a1=ones(441,8)*nan;
   a2=ones(241,3)*nan;
   filename1=[];
   filename2=[];
else
   if isempty(filename1),
      a1=ones(441,8)*nan;
      filename1=[];
   else   
      load(filename1);
      len=strfind(filename1, '.');
      filename1=filename1(1:len-1);      
      eval(['a1=' filename1 ';']);
      eval(['clear ' filename1 ';']);
      len3=strfind(filename1, '_');
      str3=filename1(1:len3-1);    
   end
   if isempty(filename2),
      a2=ones(241,3)*nan;
      filename2=[];
   else 
      load(filename2);
      len=strfind(filename2, '.');
      filename2=filename2(1:len-1);      
      eval(['a2=' filename2 ';']);
      eval(['clear ' filename2 ';']);
      len4=strfind(filename2, '_');
      str4=filename2(1:len4-1);
   end
end
% read Field data for Iwanai and Esashi
load Iwanai_Esashi
% to_E  to_I  yo_E  yo_I
% BENCHMARK NAME
str = '- OKUSHIRI - July,12 1993';
% COLOR LAB. RESULT
cc=[1 0 0];
% -------------------------------------------------------------------------
% CASES  OKUSHIR TIDE GAUGE
% -------------------------------------------------------------------------
err=nan;
err1=nan;
fig1=figure;
set(fig1,'name',str,...
         'position',[150 100 1400 900]);
annotation(fig1,'textbox',...
    'String',[str],...
    'FontWeight','bold',...
    'FontSize',16,...
    'FitHeightToText','off',...
    'LineStyle','none',...
    'Position',[0.40 0.94 0.46 0.05]);
% -------------------------------------------------------------------------
%   IWANAI 
% -------------------------------------------------------------------------
if ~isempty(filename1),
    I=find(to_I>=0.0 & to_I<=40.0);
    J=find(a1(:,1)>=0.0 & a1(:,1)<=40.0);    
    to_In=to_I(I);
    yo_In=yo_I(I);
    Norm=max(yo_In)-min(yo_In);
    yo_Ir=interp1(a1(:,1),a1(:,2),to_In);
    ina= find(~isnan(yo_Ir));
%    figure
%    plot(to_In,yo_Ir,'r.')            % Model results at Exp. x location 
%    hold on 
%    plot(a1(:,1),a1(:,2),'g.')        % model results
%    plot(to_I,yo_I,'b.')              % Exp. Results
%    return
%--------------------------------------------------------------------------
   % GET NRMSD error
   N=length(ina);
   err = sqrt( sum( (yo_Ir(ina)-yo_In(ina)).^2) / N )/Norm*100.0;
   % GET MAX WAVE AMPLITUDE or RUNUP
   %err1=abs(max(yo_In(ina))-max(a1(J,2)))/max(yo_In(ina))*100.0;
   err1=abs(max(yo_In(ina))-max(a1(J,2)))/Norm*100.0;
   Norm
   [1 max(yo_In(ina)) max(a1(J,2))]
end
% ------------------------------------------------------------------------
%  Plot  Results for IWANAI
% ------------------------------------------------------------------------
subplot(2,1,1)
if ~isempty(filename1),
   plot(a1(:,1),a1(:,2),'k--','linewidth',3)
   hold on
end
%plot(to_I(1:1:end), yo_I(1:1:end),'.','MarkerSize',12,'color',cc)
plot(to_I(1:1:end), yo_I(1:1:end),'-','linewidth',2,'color',cc)
ax1=[0 60 -4.00 4.00];
axis(ax1)
ylabel('\zeta(m)','fontsize',16)
if ~isempty(filename1),
   legend(str3,'GAUGE DATA','location','best');
else
   legend('GAUGE DATA','location','best');
end
text(0.80*ax1(2),0.9*ax1(4),'IWANAI','fontsize',14)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str4 ' - NRMSD=' num2str(err,'%3.0f') '%'],'fontsize',14)
end
if ~isnan(err1)
text(0.65*ax1(2),0.40*ax1(4),[str4 ' - ERR Max Wave Amp =' num2str(err1,'%3.0f') '%'],'fontsize',14)
end
set(gca,'xtick',[0:10:60],...
        'xticklabel',[],...
        'Ytick',[-4.00:1.00:4.00],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on
% -------------------------------------------------------------------------
%   ESASHI 
% -------------------------------------------------------------------------
err=nan;
err1=nan;
if ~isempty(filename2),
    I=find(to_E>=0.0 & to_E<=40.0);
    J=find(a2(:,1)>=0.0 & a2(:,1)<=40.0);    
    to_En=to_E(I);
    yo_En=yo_E(I);
    Norm=max(yo_En)-min(yo_En);
    yo_Er=interp1(a2(:,1),a2(:,2),to_En);
    ina= find(~isnan(yo_Er));
%    figure
%    plot(to_En,yo_Er,'r.')            % Model results at Exp. x location 
%    hold on 
%    plot(a2(:,1),a2(:,2),'g.')        % model results
%    plot(to_E,yo_E,'b.')              % Exp. Results
%    return
%--------------------------------------------------------------------------
   % GET NRMSD error
   N=length(ina);
   err = sqrt( sum( (yo_Er(ina)-yo_En(ina)).^2) / N )/Norm*100.0;
   % GET MAX WAVE AMPLITUDE or RUNUP
   %err1=abs(max(yo_En(ina))-max(a2(J,2)))/max(yo_En(ina))*100.0;
   err1=abs(max(yo_En(ina))-max(a2(J,2)))/Norm*100.0;
   Norm
   [2 max(yo_En(ina)) max(a2(J,2))]
end
% ------------------------------------------------------------------------
%  Plot  Results for IWANAI
% ------------------------------------------------------------------------
subplot(2,1,2)
if ~isempty(filename2),
   plot(a2(:,1),a2(:,2),'k--','linewidth',3)
   hold on
end
%plot(to_E(1:1:end), yo_E(1:1:end),'.','MarkerSize',12,'color',cc)
plot(to_E(1:1:end), yo_E(1:1:end),'-','linewidth',2,'color',cc)
ax1=[0 60 -4.00 4.00];
axis(ax1)
ylabel('\zeta(m)','fontsize',16)
xlabel('Time(min)','fontsize',16)
text(0.80*ax1(2),0.9*ax1(4),'ESASHI','fontsize',14)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str4 ' - NRMSD=' num2str(err,'%3.0f') '%'],'fontsize',14)
end
if ~isnan(err1)
text(0.65*ax1(2),0.40*ax1(4),[str4 ' - ERR Max Wave Amp =' num2str(err1,'%3.0f') '%'],'fontsize',14)
end
set(gca,'xtick',[0:10:60],...
        'xticklabel',[0:10:60],...
        'Ytick',[-4.00:1.00:4.00],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on
