function  FieldBM_CS_01(filename)
%
% -------------------------------------------------------------------------
%               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:  Rat Islands Tsunami
% 
%               TYPE:           Field Benchmarking (FieldBM)
%                               NUMBER: 01
%               TSUNAMI SOURCE: CoSeismic wave (CS)                       
%               SCRIPT NAME   : FieldBM_CS_01.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., E.N. Bernard, V.V. Titov, U. Kâno?lu,
%               and F.I. González (2007): Standards, criteria, and
%               procedures for NOAA evaluation of tsunami numerical models.
%               NOAA Tech. Memo. OAR PMEL-135, NOAA/Pacific Marine
%               Environmental Laboratory, Seattle, WA, 55 pp.
%
%                   Titov, V.V., F.I. González, E.N. Bernard, M.C. Eble,
%               H.O. Mofjeld, J.C. %Newman, and A.J. Venturato (2005):
%               Real-time tsunami forecasting: Challenges and solutions.
%               Nat. Hazards, 35(1), Special Issue, U.S. National Tsunami
%               Hazard Mitigation Program, 41–58.
% 
%                   NOAA/PMEL Center for Tsunami research
%               http://nctr.pmel.noaa.gov/benchmark/index.html
% 
%               INCLUDED FILE WITH THIS SCRIPT: 
%               CS_01_gage.mat (Hilo gage time serie, Nov. 17/2003.
%                               Rat Islands tsunami)
%
%               Note: File mentioned above is 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/
% 
% -------------------------------------------------------------------------
%   > FieldBM_CS_01(filename)
%
%   FieldBM_CS_01 reads input files 'filename' containing 
%   tide gage dataset recorded at Hilo Bay on Nov/17/2007 during the Rat
%   Islands Tsunami. This Benchmark case is described in detail in the
%   following URL adress 
%   NOAA/PMEL(http://nctr.pmel.noaa.gov/benchmark/index.html).
%
%   filename is a user numerical model result file containig dataset of
%   surface elevation at Hilo
%
%   Input files 'filename' needs to be named as follows:
%   XXXXXX_CS_01.dat where XXXXXX is
%   the name of the numerical model that is used,
%   eg. ADCIRC_CS_01.dat.
% 
%   The input file, 'filename' must have the following
%   configuration or structure:
%
%   1- File 'filename' is a matrix ASCII files of size
%      180x2(row x column)
%   2- Column 1 is time given in second, ranging form 
%      [5.015 to 7.998333333333333] with increment of
%      0.01666666667 hour (60 sec).
%   3- Columns 2 is user numerical results (surface elev.) in meter at
%      Hilo's gage
%   +-----------------+------------------+
%   | Time(hours)     |      Gage(m)     |
%   +-----------------+------------------+
%   5.015000000000000                   0
%   5.031666666666666   0.012355555555555
%   5.048333333333333  -0.009605714285714
%   5.065000000000000   0.016423490794919
%   5.081666666666667   0.008023499376441
%   5.098333333333334  -0.006937938878834
%           .                  . 
%           .                  .
%           .                  .
%   7.898333333333333   0.014060703969055
%   7.915000000000000   0.007230160026670
%   7.931666666666667  -0.008925599455011
%   7.948333333333333   0.006420027591456
%   7.965000000000000  -0.009392380952381
%   7.981666666666667  -0.005177777777778
%   7.998333333333333                   0
%   +-----------------+------------------+   
%   
%   NOTE:
%   This script is self-containing. It loads the field data automatically to
%   facilitate user's comparison, peer review, and determine model
%   performance.
%   The field data are saved on file CS_01_gage.mat.
%   Be sure that the field data file CS_01_gage resides
%   in the same directory or subdirectory as the Matlab script
%   FieldBM_CS_01.m
%   The user's input file 'filename' must be present in the same directory
%   or subdirectory of the script file FieldBM_CS_01.m    
% 
%   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. The NRMSD is plotted in time to visualize model
%   performance in a particular moment. As it is known tsunami models
%   usually predict fairly well the leading waves but perform poorly in
%   predicting subsequent waves. In this way users can quantify model
%   performance for the first two or three waves. 
%
%   Examples:
%   > FieldBM_CS_01
%     In this example the Matlab scrip just plots laboratory results.
%   > FieldBM_CS_01('ADCIRC_CS_01.dat')
%     In this example the Matlab script plots recorded Field dataset
%     against user's model results.
% -------------------------------------------------------------------------
help FieldBM_CS_01
if nargin < 1
   a=ones(180,2)*nan;
   filename=['       '];
else  
   load(filename);
   len=strfind(filename, '.');
   filename=filename(1:len-1);      
   eval(['a=' filename ';']);
   eval(['clear ' filename ';']);
   disp('size of model-data (Row,Column)' )
   disp(size(a))
end
% read Field time serie (Hilo gage)
ax=[5 8 -0.3 0.3];
load CS_01_gage
time=time/3600.0;
is=find(time>=ax(1) & time<=ax(2));
time1=time(is,1);
height1=height(is,1);
% Remove tide (detide)
h_smooth=smooth(smooth(height1,50),100);
height1=height1-h_smooth;
Norm=max(height1)-min(height1);
% BENCHMARK NAME
str='Rat Islands Tsunami';
% BENCHMARK CASE
% COLOR LAB. RESULT (for LINE)
cc=[0.4 0.6 0.8];  % light blue-gray
%
% -------------------------------------------------------------------------
% FIGURE
% -------------------------------------------------------------------------
fig1=figure;
set(fig1,'name',str,...
         'position',[100 150 1400 900]);
annotation(fig1,'textbox',...
    'String',[str],...
    'FontWeight','bold',...
    'FontSize',20,...
    'FitHeightToText','off',...
    'LineStyle','none',...
    'Position',[0.42 0.94 0.46 0.05]);

% GET RMS error
N=numel(height1);
err = sqrt( sum( (height1-a(:,2)).^2) / N )/Norm*100.0;
% Error
% -----
for ii=1:N, 
    err1(ii,1)=sqrt( sum( (height1(1:ii,1)-a(1:ii,2)).^2) / ii )/Norm*100.0;
end
%
if ~isnan(err),
   % Plot Hilo gage, numerical result and NRMSD
   [AX,H1,H2] = plotyy([time1 a(:,1)],[height1 a(:,2)],a(:,1),err1,'plot');
   %
   set(get(AX(1),'Ylabel'),'String','Amplitude (m)','fontsize',14) 
   set(get(AX(2),'Ylabel'),'String','NRMSD%','fontsize',14,'color',[1 0 0])
   set(H1(1),'LineStyle','-','color',cc,'linewidth',4)
   set(H1(2),'LineStyle','--','color',[0 0 0],'linewidth',2)
   set(H2,'LineStyle','-.','color',[1 0 0],'linewidth',2)
   set(AX(1),'Ylim',[ax(3) ax(4)],...
             'ytick',[ax(3):0.1:ax(4)],...
             'yticklabel',[ax(3):0.1:ax(4)],...
             'Yminortick','on',...
             'fontsize',14,'Ycolor',[0 0 0])
   set(AX(1),'xtick',[ax(1):ax(2)],...
             'xticklabel',[ax(1):ax(2)],...
             'Xminortick','on',...
             'Xgrid','on',...
             'Xlim',[ax(1) ax(2)],...
             'fontsize',14)
   set(AX(2),'xtick',[],'xticklabel',[],'fontsize',14,'Ycolor',[1 0 0])
   %    
   text(0.85*ax(2),0.85*ax(4),[filename(1:7) ' - NRMSD=' num2str(err) '%'],'fontsize',16)
   %
   grid on
else
   % plot Hilo gage
   plot(time1,height1,'linewidth',4,'color',cc)
   set(gca,'ytick',[ax(3):0.1:ax(4)],...
           'yticklabel',[ax(3):0.1:ax(4)],...
           'Ylim',[ax(3) ax(4)],...
           'Yminortick','on',...
           'Ygrid','on',...
           'xtick',[ax(1):ax(2)],...
           'xticklabel',[ax(1):ax(2)],...
           'Xminortick','on',...
           'Xgrid','on',...
           'Xlim',[ax(1) ax(2)],...
           'fontsize',14)  
   ylabel('Amplitude (m)','fontsize',14)
end
legend('HILO GAGE',filename(1:7),'NRMSD(t)%','location','best');
xlabel(' t(Hours after earthquake 06:43:07 UTC, November 17, 2003)','fontsize',14)
return

