function  BP4_Lab_SW01(filename1,filename2)
% -------------------------------------------------------------------------
%               NTHMP (National Tsunami Hazard Mitigation Program)
% -------------------------------------------------------------------------
%                          2010 NTHMP WORKSHOP
%              Venue: TBA (Catalina Island, CA;
%                     Santa Monica, CA; Irvine, CA; Corvallis, OR;
%                     Boulder, CO.
%                            NOV. 16th-18th 2010
%                               ---- o ----
%               MATLAB SCRIPT:  Solitary Wave on a Simple Beach
%               TYPE:           Laboratory Benchmarking Problem (Lab)
%                               NUMBER: 01
%               TSUNAMI SOURCE: Incident Solitary Wave (SW)                       
%               SCRIPT NAME   : BP4_Lab_SW01.m
%                
%               AUTHOR(S)     : D. Nicolsky, C. Synolakis, S Grilli
%                               and J. Horrillo.
%   	        
%                              Copyright 2010 NTHMP/NOAA
%                              Revision: 0
%                              Date: Sep/2010
%               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: 
%
%               SW_01_00185.mat (lab. data file for Time evolution 
%               of H= 0.0185 initial wave over a sloping beach with
%               cot(beta)=19.85. Constant depth-segment starts at
%               Xo=19.85. Markers show experimental results of Synolakis 
%               (1986, 1987). Experimental data is provided from t=30 to
%               t=70 with 10 increments.) 
%
%               SW_01_03.mat(lab. data file for Time evolution of H=0.30
%               initial wave over a sloping beach with  cot(beta)=19.85 
%               Constant-depth segment starts at Xo=19.85. Markers show
%               experimental results of Synolakis (1987).
%               Experimental data is provided from t=15 to t=30 with 5 
%               increments.
%               
%               Note: Above files (both) are needed to run successfully this
%                     script.
% 
%               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/
% 
% -------------------------------------------------------------------------
%   > BP4_Lab_SW01(filename1,filename2)
%
%   BP4_Lab_SW01 reads input files 'filename1' and 'filename2' containing 
%   user model results and compares them with laboratory data described
%   in the Benchmark "Solitary Wave on a simple Beach", NOAA/PMEL,
%   (http://nctr.pmel.noaa.gov/benchmark/index.html).
%
%   filename1 and filename2 are model generated files corresponding to
%   Case H=0.0185m and H=0.30m respectively.
%    
%   Input files 'filename1' and 'filename2' need to be named as follows:
%   XXXXXX_SW_01_00185.dat and XXXXXX_SW_01_03.dat, where XXXXXX is
%   the name of the numerical model that is tested,
%   eg. NEOWAVE_SW_01_00185.dat and NEOWAVE_SW_01_03.dat.
% 
%   The input files, 'filename1' and 'filename2' must have the following
%   configuration:
%
%   1- Files 'filename1' is a matrix ASCII files of size
%      441x6(row x column) as shown below (no heading).
%      x ranges= [-2 20]  with icrements of 0.05.
%      Snapshot of the free surface are
%      recorded at time t=[ 30, 40, 50, 60, 70]  
%
%   filename1
%   +-------+---------+---------+----------+---------+---------+
%   |  x()  Amp(t=30) Amp(t=40)  Amp(t=50) Amp(t=60) Amp(t=70)
%   +-------+---------+---------+----------+---------+---------+ 
%    -2.00      0.000     0.297     0.000    -0.069     -0.260
%    -1.95      0.000     0.396     0.000    -0.970     -0.564
%    -1.90      0.000     0.495     0.000     0.772     -0.743
%    -1.85      0.000     0.594     0.000     0.178     -0.056
%        .          .         .         .         .          .
%        .          .         .         .         .          .
%        .          .         .         .         .          .
%    19.85	    0.546     0.495     0.783    -0.762     -0.760
%    19.90	    0.456     0.786     0.438    -0.678     -0.768
%    19.95	   -0.679    -0.569     0.907     0.322     -0.529
%    20.00	    0.456    -0.045     0.000     0.772     -0.067
%   +-------+---------+---------+----------+---------+---------+
%
%   2- Files 'filename2' is a matrix ASCII files of size
%      601x5(row x column) as shown below (no heading).
%      x ranges= [-10 20] with icrements of 0.05.
%      Snapshot of the free surface are
%      recorded at time t=[ 15, 20, 25, 30] 
%
%   filename2
%   +-------+---------+---------+----------+---------+
%   |  x()  Amp(t=15) Amp(t=20)  Amp(t=25) Amp(t=30)
%   +-------+---------+---------+----------+---------+ 
%   -10.00      0.000     0.297     0.000     3.069   
%    -9.95      0.000     0.396     0.000     2.970    
%    -9.90      0.000     0.495     0.000     2.772  
%    -9.85      0.000     0.594     0.000     2.178
%        .          .         .         .         .   
%        .          .         .         .         .  
%        .          .         .         .         .  
%    19.85	    0.546     0.495     0.783    -1.762 
%    19.90	    0.456     0.786     0.438    -2.678 
%    19.95	   -0.679    -0.569     0.907     0.322 
%    20.00	    0.456    -0.045     0.000     0.772 
%   +-------+---------+---------+----------+---------+
%
%   NOTE:
%   This script is self-containing. It loads the lab. data automatically to
%   facilitate user's comparison, peer review and determine model 
%   performance.
%   The Experimental data have been saved on files SW_01_00185.mat and 
%   SW_01_03.mat for case H=0.0185m and H=0.30m respectively.
%   Be sure that the lab data files SW_01_00185.mat and SW_01_03.mat reside
%   in the same directory or subdirectory as the Matlab script
%   (BP4_Lab_SW01.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 values actually observed in the experiment
%   or in the field.
% 
%   Examples:
%   > BP4_Lab_SW01
%     In this example the Matlab scrip just plots laboratory results.
%   > BP4_Lab_SW01('NEOWAVE_SW_01_00185.dat','NEOWAVE_SW_01_03.dat')
%     In this example the Matlab scrip plots both cases laboratory results
%     and user's model results.
%   > BP4_Lab_SW01('NEOWAVE_SW_01_00185.dat',[])
%     In this example the Matlab scrip only plots laboratory results and
%     and user's model result for case H=0.0185.
%   > BP4_Lab_SW01([],'NEOWAVE_SW_01_03.dat')
%     In this example the Matlab scrip plots only laboratory results and
%     user's model result for case H=0.30.
% -------------------------------------------------------------------------
help BP4_Lab_SW01
str3=''
if nargin < 1
   a1=ones(441,6)*nan;
   a2=ones(601,5)*nan;
   filename1=['       '];
   filename2=['       '];
else
   if isempty(filename1),
      a1=ones(441,6)*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(601,5)*nan;
      filename2=['       '];
   else 
      load(filename2);
      len=strfind(filename2, '.');
      filename2=filename2(1:len-1);      
      eval(['a2=' filename2 ';']);
      eval(['clear ' filename2 ';']);
      len3=strfind(filename2, '_');
      str3=filename2(1:len3-1);
   end
end
% read lab data for CASE H=0.0185m
load SW_01_00185.mat
% read lab data for CASE H=0.3m
load SW_01_03.mat
% BENCHMARK NAME
str='Solitary Wave on a Canonical Beach';
% BENCHMARK CASE
str1='CASE     H = 0.0185';
str2='CASE     H = 0.30  ';
% COLOR LAB. RESULT
cc=[1 0 0];
% Bathy
% ------------
X=[-10:1: 20];
slp=1/19.85;
Y=-slp*X;
%
% -------------------------------------------------------------------------
% CASE H=0.0185m
% -------------------------------------------------------------------------
fig2=figure;
set(fig2,'name',str,...
         'position',[100 150 1400 900]);
annotation(fig2,'textbox',...
    'String',[str '  ' str1],...
    'FontWeight','bold',...
    'FontSize',16,...
    'FitHeightToText','off',...
    'LineStyle','none',...
    'Position',[0.29 0.94 0.46 0.05]);
% plot snapshot t=30 -lab result- H=0.0185m 
subplot(5,1,1)
Norm=max(a_30)-min(a_30);
a_30n=interp1(x_30,a_30,a1(:,1));
ina=find(~isnan(a_30n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a_30n(ina)-a1(ina,2)).^2) / N )/Norm*100.0;
%-------------------------------------------------------------------------
plot(X,Y,'b-','linewidth',3)
hold on
plot(x_30, a_30,'.','MarkerSize',14,'color',cc)
plot(a1(:,1),a1(:,2),'k--','linewidth',3)
ax1=[-2 20 -0.02 0.04];
axis(ax1)
legend('SLOPE','EXPERIMENT',str3,'location','best');
text(0.82*ax1(2),0.88*ax1(4),'t=30','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-2:2:20],...
        'xticklabel',[],...
        'Ytick',[-0.04:0.02:0.10],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on
%-------------------------------------------------------------------------
% plot snapshot t=40 -lab result- H=0.0185m 
subplot(5,1,2)
Norm=max(a_40)-min(a_40);
a_40n=interp1(x_40,a_40,a1(:,1));
ina=find(~isnan(a_40n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a_40n(ina)-a1(ina,3)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x_40, a_40,'.','MarkerSize',14,'color',cc)
plot(a1(:,1),a1(:,3),'k--','linewidth',3)
ax1=[-2 20 -0.02 0.04];
axis(ax1)
text(0.82*ax1(2),0.88*ax1(4),'t=40','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-2:2:20],...
        'xticklabel',[],...
        'Ytick',[-0.04:0.02:0.10],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on

% plot snapshot t=50 -lab result- H=0.0185m 
subplot(5,1,3)
Norm=max(a_50)-min(a_50);
a_50n=interp1(x_50,a_50,a1(:,1));
ina=find(~isnan(a_50n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a_50n(ina)-a1(ina,4)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x_50, a_50,'.','MarkerSize',14,'color',cc)
plot(a1(:,1),a1(:,4),'k--','linewidth',3)
ax1=[-2 20 -0.02 0.08];
axis(ax1)
ylabel('\eta ()','fontsize',16,'rotation',90);
text(0.82*ax1(2),0.88*ax1(4),'t=50','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-2:2:20],...
        'xticklabel',[],...
        'Ytick',[-0.04:0.02:0.10],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on

% plot snapshot t=60 -lab result- H=0.0185m 
subplot(5,1,4)
Norm=max(a_60)-min(a_60);
a_60n=interp1(x_60,a_60,a1(:,1));
ina=find(~isnan(a_60n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a_60n(ina)-a1(ina,5)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x_60, a_60,'.','MarkerSize',14,'color',cc)
plot(a1(:,1),a1(:,5),'k--','linewidth',3)
ax1=[-2 20 -0.02 0.10];
axis(ax1)
text(0.82*ax1(2),0.88*ax1(4),'t=60','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-2:2:20],...
        'xticklabel',[],...
        'Ytick',[-0.04:0.02:0.10],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on

% plot snapshot t=70 -lab result- H=0.0185m 
subplot(5,1,5)
Norm=max(a_70)-min(a_70);
a_70n=interp1(x_70,a_70,a1(:,1));
ina=find(~isnan(a_70n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a_70n(ina)-a1(ina,6)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x_70, a_70,'.','MarkerSize',14,'color',cc)
plot(a1(:,1),a1(:,6),'k--','linewidth',3)
ax1=[-2 20 -0.04 0.06];
axis(ax1)
xlabel('x()','fontsize',16)
text(0.82*ax1(2),0.88*ax1(4),'t=70','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-2:2:20],...
        'Ytick',[-0.04:0.02:0.10],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on

% -------------------------------------------------------------------------
% CASE H=0.30m
% -------------------------------------------------------------------------
fig1=figure;
set(fig1,'name',str,...
         'position',[150 100 1400 900]);
annotation(fig1,'textbox',...
    'String',[str '  ' str2],...
    'FontWeight','bold',...
    'FontSize',16,...
    'FitHeightToText','off',...
    'LineStyle','none',...
    'Position',[0.29 0.94 0.46 0.05]);
%
% plot snapshot t=15 -lab result- H=0.30m 
subplot(4,1,1)
Norm=max(a1_15)-min(a1_15);
a1_15n=interp1(x1_15,a1_15,a2(:,1));
ina=find(~isnan(a1_15n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a1_15n(ina)-a2(ina,2)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x1_15, a1_15,'.','MarkerSize',14,'color',cc)
plot(a2(:,1),a2(:,2),'k--','linewidth',3)
ax1=[-10 20 -0.10 0.50];
axis(ax1)
legend('SLOPE','EXPERIMENT',str3,'location','best');
text(0.80*ax1(2),0.9*ax1(4),'t=15','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-10:5:20],...
        'xticklabel',[],...
        'Ytick',[-0.10:0.10:0.50],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on

% plot snapshot t=20 -lab result- H=0.30m 
subplot(4,1,2)
Norm=max(a1_20)-min(a1_20);
a1_20n=interp1(x1_20,a1_20,a2(:,1));
ina=find(~isnan(a1_20n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a1_20n(ina)-a2(ina,3)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x1_20, a1_20,'.','MarkerSize',14,'color',cc)
plot(a2(:,1),a2(:,3),'k--','linewidth',3)
ax1=[-10 20 -0.10 0.50];
axis(ax1)
text(0.80*ax1(2),0.9*ax1(4),'t=20','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-10:5:20],...
        'xticklabel',[],...
        'Ytick',[-0.10:0.10:0.50],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on


% plot snapshot t=25 -lab result- H=0.30m 
subplot(4,1,3)
Norm=max(a1_25)-min(a1_25);
a1_25n=interp1(x1_25,a1_25,a2(:,1));
ina=find(~isnan(a1_25n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a1_25n(ina)-a2(ina,4)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x1_25, a1_25,'.','MarkerSize',14,'color',cc)
plot(a2(:,1),a2(:,4),'k--','linewidth',3)
ax1=[-10 20 -0.10 0.50];
axis(ax1)
ylabel('\eta ()','fontsize',16,'rotation',90);
text(0.80*ax1(2),0.9*ax1(4),'t=25','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-10:5:20],...
        'xticklabel',[],...
        'Ytick',[-0.10:0.10:0.50],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on

% plot snapshot t=30 -lab result- H=0.30m 
subplot(4,1,4)
Norm=max(a1_30)-min(a1_30);
a1_30n=interp1(x1_30,a1_30,a2(:,1));
ina=find(~isnan(a1_30n));
% GET RMS error
N=length(ina);
err = sqrt( sum( (a1_30n(ina)-a2(ina,5)).^2) / N )/Norm*100.0;
plot(X,Y,'b-','linewidth',3)
hold on
plot(x1_30, a1_30,'.','MarkerSize',14,'color',cc)
plot(a2(:,1),a2(:,5),'k--','linewidth',3)
ax1=[-10 20 -0.10 0.50];
axis(ax1)
xlabel('x()','fontsize',16)
text(0.80*ax1(2),0.9*ax1(4),'t=30','fontsize',16)
if ~isnan(err)
text(0.70*ax1(2),0.65*ax1(4),[str3 ' - NRMSD=' num2str(err) '%'],'fontsize',16)
end
set(gca,'xtick',[-10:5:20],...
        'Ytick',[-0.10:0.10:0.50],...
        'Xminortick','on',...
        'Yminortick','on',...
        'fontsize',14)
grid on
 return

