commit af986a42977a692fb925feec5cc9ef8b32dcf5a0
Author: Joshua Simmons <z3288267@ad.unsw.edu.au>
Date: Tue Feb 20 10:47:20 2018 +1100
init commit
diff --git a/RUN_sbeach_model.m b/RUN_sbeach_model.m
new file mode 100644
index 0000000..51f653a
--- /dev/null
+++ b/RUN_sbeach_model.m
@@ -0,0 +1,120 @@
+% Written by Joshua Simmons 11/02/2015 and updated 17/11/2017
+% to convert from fielddata storage struct to SBEACH toolbox.
+
+clear
+clc
+
+%% Settings
+%profiledata locations
+basedatapath = '.\example\input_data\';
+setupoutputpath = '.\temp\';
+infooutpath = '.\example\setup_info\';
+
+sbeachexepath = '.\sbeach_exe\final_exe\SBEACH.exe';
+
+%meta
+site = 'narrabeen';
+stormid = '2015';
+
+%write file to
+writetofile = 1; %boolean 1 = write to file, 0 = just SBEACHout.
+
+% other parameters;
+sbts = 2; %time step length minutes
+%% Code - format the data
+%load profile data - account for different data names - lin/poly
+pfdata = load([basedatapath 'morphology\processed\morphology_' site '_' stormid '.mat']);
+fns1 = fieldnames(pfdata);
+pfdata = pfdata.(fns1{1});
+
+%load wavedata
+wavedata = load([basedatapath 'fielddata\wave\processed\wavedata_' site '_' stormid '.mat']);
+fns2 = fieldnames(wavedata);
+wavedata = wavedata.(fns2{1});
+
+%load tidedata
+tidedata = load([basedatapath 'fielddata\tide\processed\tidedata_' site '_' stormid '.mat']);
+fns3 = fieldnames(tidedata);
+tidedata = tidedata.(fns3{1});
+
+% load the profileinfo file
+profileinfo = load([infooutpath 'profileinfo_' site '.mat']);
+fns2 = fieldnames(profileinfo);
+profileinfo = profileinfo.(fns2{1});
+
+%check for dimension match
+if length(pfdata.data) ~= length(wavedata.ns.data)
+ error('Wave data information does not match the morphology data')
+end
+
+%get the depth contour
+sbdepth = wavedata.ns.contour(1:2);
+
+pathrun = cell(length(pfdata.data),1);
+
+%do the modelling
+fprintf('Setting up %d SBEACH models...', length(pfdata.data))
+for ii = 1:length(pfdata.data)
+ if ~isnan(pfdata.data(ii).initial.x)
+ fprintf('%d..',ii)
+ pfid = ['PF' num2str(ii)];
+
+ %get profile info
+ pfx = pfdata.data(ii).initial.x;
+ pfz = pfdata.data(ii).initial.z;
+
+ %angle of pf from true north
+ pfangle = profileinfo.data(ii).TN; %angle from true north
+
+ outpath = [setupoutputpath 'sb_' pfid '\'];
+% outpath = [setupoutputpath 'sb_' site '_' stormid '_' pfid '\'];
+
+ sb_model = fsbeach_setupmodel(...
+ 'bathy', { 'x', pfx, 'z', pfz}, ...
+ 'waves', { 'time', wavedata.ns.data(ii).datenums, 'Hsig', wavedata.ns.data(ii).Hsig, 'Tp1', wavedata.ns.data(ii).Tp1, ...
+ 'ang', wavedata.ns.data(ii).Wdir, 'pfang', pfangle}, ...
+ 'tide', { 'time', tidedata.time,'elev', tidedata.tide}, ...
+ 'settings', { 'DT', sbts, 'K', 2.5e-6, 'BMAX', 28.5, 'D50', 0.455 ...
+ }, ...
+ 'path', outpath, ...
+ 'name', ['sb_' pfid] ...
+ );
+
+ % variable grid can be specified in settings with e.g.:
+ % 'NGRID', 3, 'DXV', [1 5 10], 'NDXV', [115 30 28]
+
+ % hard bottom can be specified in bathy with e.g.:
+ % 'hbz', pfz-1
+
+ measured(ii).x = pfdata.data(ii).measured.x;
+ measured(ii).z = pfdata.data(ii).measured.z;
+
+ pathrun{ii} = sb_model.path;
+ end
+end
+
+fprintf('\n')
+fprintf('Running %d SBEACH models...', length(pathrun))
+%now run sbeach EXE on the files setup
+for ii = 1:length(pathrun)
+ fprintf('%d..',ii)
+ [fail, cmdout] = fsbeach_runmodel(pathrun{ii},sbeachexepath);
+ if fail
+ error(cmdout);
+ end
+end
+
+fprintf('\n')
+fprintf('Collecting SBEACH results...\n')
+%now collect the sbeach results
+sbout = fsbeach_collectresults(pathrun);
+
+fprintf('Getting measured data...\n')
+%attached the measured data
+sbout.measured = measured;
+
+fprintf('Processing SBEACH results...\n')
+%evaulate results
+sbout = fsbeach_processresults(sbout);
+
+fprintf('Done!\n')
\ No newline at end of file
diff --git a/SBEACH_GLUE_setup.m b/SBEACH_GLUE_setup.m
new file mode 100644
index 0000000..42e4b7b
--- /dev/null
+++ b/SBEACH_GLUE_setup.m
@@ -0,0 +1,169 @@
+% Written by Joshua Simmons 11/02/2015 and updated on 11/01/2016
+% to convert from fielddata storage struct to SBEACH compatible format.
+
+clear
+clc
+
+%% Settings
+%all the initial files location
+%specify folder with Mitch's bathymetry and field data setup for xbeach.
+setupoutputpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\setup_folders\';
+setupbasepath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\';
+
+basebatchrunpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_results\src\';
+
+% meta
+site = 'deewhy';
+
+% stormid = '2007';
+% stormid = '2011';
+% stormid = '2014';
+stormid = '2015';
+
+%profile info - CHECK
+pfno = 5;
+
+% model parameters
+%parameter vals - leave as '' if params to be generated
+paramload = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\params\bigparams_K_EPS_BMAX.mat';
+% paramload = '';
+
+%or specify stats
+paramsname = {'K','EPS','BMAX'}; %K divided by 3600 below to bring into seconds, Cs for R2;
+% paramsrange = {[0.5e-6 2.5e-6],[0.001 0.003],[15 30]}; % give range as [min max] then fed into rand()
+paramsrange = {[1.75e-6 1.75e-6],[0.002 0.002],[28.5 28.5]}; % give range as [min max] then fed into rand()
+nruns = 10000; %number of parameter conmbinations
+
+nrunst = 1;
+nrunend = 10000;
+
+%% Code - format the data
+%generate param combos
+if isempty(paramload)
+ for aa = 1:length(paramsname)
+ tmpmin = paramsrange{aa}(1);
+ tmpmax = paramsrange{aa}(2);
+ params(:,aa) = (tmpmax-tmpmin).*rand(nruns,1) + tmpmin;
+ end
+ bigparams.name = paramsname;
+ bigparams.value = params;
+else
+ bigparams = load(paramload);
+ fnsparams = fieldnames(bigparams);
+ bigparams = bigparams.(fnsparams{1,1});
+ params = bigparams.value;
+ if sum(strcmp(bigparams.name,paramsname)) ~= length(paramsname)
+ error('Parameter names in loaded file do not match those specified')
+ end
+end
+
+%load pf info
+profileinfo = load([setupbasepath 'setup_info\profileinfo_' site '.mat']);
+fns2 = fieldnames(profileinfo);
+profileinfo = profileinfo.(fns2{1});
+%load cfg info
+cfginfo = load([setupbasepath 'setup_info\cfginfo_' site '_' stormid '_PF' num2str(pfno) '.mat']);
+fns3 = fieldnames(cfginfo);
+cfginfo = cfginfo.(fns3{1});
+
+%get the names of the files copying across
+dns1 = dir([setupoutputpath 'sb_' site '_' stormid '_PF' num2str(pfno) '\']);
+dns1dirs = {dns1.name}; % Get a list of the subdirectories
+dns1ind2 = ~ismember(dns1dirs,{'.','..'}); % Find index of subdirectories
+dns1 = dns1(dns1ind2);
+
+allfilelocs = cell(nruns,1);
+% not done VVVV
+%do the modelling
+for ii = nrunst:nrunend
+ CFGname = ['PF' num2str(pfno) '_Run_' num2str(ii) '.CFG'];
+ tmpprintdir = [basebatchrunpath site '_' stormid '_PF' num2str(pfno) '\PF' num2str(pfno) '_Run_' num2str(ii) '\'];
+ if ~exist(tmpprintdir,'dir')
+ mkdir(tmpprintdir)
+ end
+ %now copy the setup files across
+ copyfile([setupoutputpath 'sb_' site '_' stormid '_PF' num2str(pfno) '\*'],tmpprintdir);
+
+ % now change the contents of CFG file
+ fin = fopen([setupbasepath 'base.CFG']) ;
+ fout = fopen([tmpprintdir 'temp.CFG'], 'w');
+ tline = fgetl(fin);
+ while ischar(tline)
+ %replace:
+ %filename
+ if (~isempty(strfind(tline, '%RPLCNAME%')))
+ tline = strrep(tline, '%RPLCNAME%', ['_Run_' num2str(ii)]);
+ end
+ %profile ID
+ if (~isempty(strfind(tline, '%RPLCPFID%')))
+ tline = strrep(tline, '%RPLCPFID%', ['PF' num2str(pfno)]);
+ end
+ %run number
+ if (~isempty(strfind(tline, '%RPLCRNUM%')))
+ tline = strrep(tline, '%RPLCRNUM%', num2str(ii));
+ end
+ %calculation cells
+ %number
+ if (~isempty(strfind(tline, '%RPLCNDX%')))
+ tline = strrep(tline, '%RPLCNDX%', num2str(cfginfo.NDX));
+ end
+ %landward boundary
+ if (~isempty(strfind(tline, '%RPLCXSTART%')))
+ tline = strrep(tline, '%RPLCXSTART%', num2str(cfginfo.XSTART));
+ end
+ %spacing between
+ if (~isempty(strfind(tline, '%RPLCDXC%')))
+ tline = strrep(tline, '%RPLCDXC%', num2str(cfginfo.DXC));
+ end
+
+ %calculation time
+ %number
+ if (~isempty(strfind(tline, '%RPLCNDT%')))
+ tline = strrep(tline, '%RPLCNDT%', num2str(cfginfo.NDT));
+ end
+ %steps
+ if (~isempty(strfind(tline, '%RPLCDT%')))
+ tline = strrep(tline, '%RPLCDT%', num2str(cfginfo.DT));
+ end
+
+ %steps
+ if (~isempty(strfind(tline, '%RPLCD50%')))
+ tline = strrep(tline, '%RPLCD50%', num2str(profileinfo.data(pfno).d50*1000));
+ end
+
+ %check for params
+ for aa = 1:length(paramsname)
+ if (~isempty(strfind(tline, ['%RPLC' paramsname{aa} '%'])))
+ tline = strrep(tline, ['%RPLC' paramsname{aa} '%'], num2str(params(ii,aa)));
+ end
+ end
+ %now print the final string
+ fprintf(fout,'%s\r\n',tline);
+ tline = fgetl(fin);
+ end
+ fclose(fin);
+ fclose(fout);
+
+ %rename CFG to original name
+ try
+ movefile([tmpprintdir 'temp.CFG'],[tmpprintdir CFGname])
+ catch
+ fileattrib([tmpprintdir 'temp.CFG'],'+w')
+ movefile([tmpprintdir 'temp.CFG'],[tmpprintdir CFGname])
+ end
+
+ %change file name
+ for jj = 1:length(dns1)
+ movefile([tmpprintdir dns1(jj).name],[tmpprintdir 'PF' num2str(pfno) '_Run_' num2str(ii) dns1(jj).name(end-3:end)])
+ end
+
+ %save location of file
+ allfilelocs{ii,1} = ['PF' num2str(pfno) '_Run_' num2str(ii) '\PF' num2str(pfno) '_Run_' num2str(ii)];
+end
+
+%write allfilelocs to .txt
+fid = fopen([basebatchrunpath site '_' stormid '_PF' num2str(pfno) '\SBEACHrun_loc.txt'],'w');
+for kk = 1:length(allfilelocs)
+ fprintf(fid,'%s\r\n',allfilelocs{kk});
+end
+fclose(fid);
\ No newline at end of file
diff --git a/SBEACH_storm_setup.m b/SBEACH_storm_setup.m
new file mode 100644
index 0000000..f02f9e3
--- /dev/null
+++ b/SBEACH_storm_setup.m
@@ -0,0 +1,140 @@
+% Written by Joshua Simmons 11/02/2015 and updated 20/01/2016
+% to convert from fielddata storage struct to SBEACH compatible format.
+
+clear
+clc
+
+%% Settings
+%profiledata locations
+basedatapath = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\';
+setupoutputpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\setup_folders\';
+infooutpath = 'D:\7 Erosion Models\a_nb_model_comparison\sbeach\model_setup\setup_info\';
+
+%meta
+site = 'monavale';
+
+% stormid = '2007';
+% stormid = '2011';
+% stormid = '2014';
+stormid = '2015';
+
+%write file to
+writetofile = 1; %boolean 1 = write to file, 0 = just SBEACHout.
+
+% other parameters;
+sbts = 5; %time step length minutes
+%% Code - format the data
+%load profile data - account for different data names - lin/poly
+pfdata = load([basedatapath 'morphology\processed\morphology_' site '_' stormid '.mat']);
+fns1 = fieldnames(pfdata);
+pfdata = pfdata.(fns1{1});
+
+%load wavedata
+wavedata = load([basedatapath 'fielddata\wave\processed\wavedata_' site '_' stormid '.mat']);
+fns2 = fieldnames(wavedata);
+wavedata = wavedata.(fns2{1});
+
+%load tidedata
+tidedata = load([basedatapath 'fielddata\tide\processed\tidedata_' site '_' stormid '.mat']);
+fns3 = fieldnames(tidedata);
+tidedata = tidedata.(fns3{1});
+
+% load the profileinfo file
+profileinfo = load([infooutpath 'profileinfo_' site '.mat']);
+fns2 = fieldnames(profileinfo);
+profileinfo = profileinfo.(fns2{1});
+
+%check for dimension match
+if length(pfdata.data) ~= length(wavedata.ns.data)
+ error('Wave data information does not match the morphology data')
+end
+
+%check start time
+if tidedata.time(1) ~= wavedata.os.datenums(1)
+ error('Start times for wave and tide data dont match')
+end
+if tidedata.time(end) ~= wavedata.os.datenums(end)
+ error('End times for wave and tide data dont match. Make sure start/end times are hourly (for wave data).')
+end
+
+%check and allocate time parameters;
+sbtidetime = round((tidedata.time - tidedata.time(1))*24*60);
+%check time steps
+diffsbtidetime = diff(sbtidetime);
+if any(diffsbtidetime~=diffsbtidetime(1))
+ error('Inconsistent time steps in tide data');
+end
+
+%now check the wave data lines up (os only required the rest if copied)
+sbwavetime = round((wavedata.os.datenums - wavedata.os.datenums(1))*24*60);
+diffsbwavetime = diff(sbwavetime);
+if any(diffsbwavetime~=diffsbwavetime(1))
+ error('Inconsistent time steps in wave data');
+end
+
+%get the depth contour
+sbdepth = wavedata.ns.contour(1:2);
+
+%do the modelling
+for ii = 1:length(pfdata.data)
+ if ~isnan(pfdata.data(ii).initial.x)
+ pfid = ['PF' num2str(ii)];
+ %get profile info
+ prifile = [pfdata.data(ii).initial.x pfdata.data(ii).initial.z]; %initial pf
+
+ %angle of pf from true north
+ pfangle = profileinfo.data(ii).TN; %angle from true north
+
+ %create the other information required for files
+ % for .WAV and .ANG file - wave height/period and angle files
+ wavfile = [wavedata.ns.data(ii).Hsig wavedata.ns.data(ii).Tp1];
+ angfile = -1*(wavedata.ns.data(ii).Wdir-pfangle); %correct to SBEACH wave angle convention as per report 3 pg 19
+ wavets = sbwavetime(2)-sbwavetime(1); %wave time step minutes
+
+ % for .ELV file - water elevation file
+ elvfile = tidedata.tide;
+ tidets = sbtidetime(2)-sbtidetime(1); %tide time step
+
+ totaltimewav = sbwavetime(end);
+ totaltimetide = sbtidetime(end);
+ if totaltimewav - totaltimetide ~= 0
+ disp('Warning: difference between wave and tide total times')
+ end
+ totaltime = totaltimetide;
+ totalmodeltimesteps = totaltime/sbts;
+
+ disp([pfid ':'])
+ disp(['A.10 - Total Time steps: ' num2str(totalmodeltimesteps)])
+ disp(['B.6 - DTWAV: ' num2str(wavets)])
+ disp(['B.11 - DTANG: ' num2str(wavets)])
+ disp(['B.12 - DMEAS: ' sbdepth])
+ disp(['B.20 - DTELEV: ' num2str(tidets)])
+
+ outpath = [setupoutputpath 'sb_' site '_' stormid '_' pfid '\'];
+
+ cfginfo.NDT = totalmodeltimesteps;
+ cfginfo.DT = sbts;
+ cfginfo.DXC = pfdata.data(ii).initial.x(2)-pfdata.data(ii).initial.x(1);
+ cfginfo.NDX = length(pfdata.data(ii).initial.x)-1;
+ [cfginfo.XSTART,~] = min(pfdata.data(ii).initial.x);
+ cfginfo.DTWAV = wavets;
+ cfginfo.DTANG = wavets;
+ cfginfo.DMEAS = str2double(sbdepth);
+ cfginfo.DTELEV = tidets;
+
+
+ %write to files
+ if writetofile
+ if ~exist(outpath, 'dir')
+ mkdir(outpath)
+ else
+ delete([outpath '\*'])
+ end
+ writetofilepri(prifile, outpath, [site '_' stormid '_' pfid])
+ writetofilewav(wavfile, outpath, [site '_' stormid '_' pfid])
+ writetofileang(angfile, outpath, [site '_' stormid '_' pfid])
+ writetofileelv(elvfile, outpath, [site '_' stormid '_' pfid])
+ save([infooutpath 'cfginfo_' site '_' stormid '_' pfid '.mat'],'cfginfo')
+ end
+ end
+end
\ No newline at end of file
diff --git a/calc_TNorientation.m b/calc_TNorientation.m
new file mode 100644
index 0000000..92917b1
--- /dev/null
+++ b/calc_TNorientation.m
@@ -0,0 +1,38 @@
+% Written by Joshua Simmons 20/01/16 to calculate angle for TN for SBeach
+% modelling purposes.
+
+% NOTE: this code assumes that the furthest offshore point (-10m) is the
+% x(end) and the furthest inshore point is in x(1) (in line with the xbeach
+% co-ordinate system)
+
+clear
+clc
+
+% load('D:\7 Erosion Models\a_nb_model_comparison\xbeach\other_data\narra_xyz_for_alfa.mat')
+%
+% for ii = 1:length(pfxy)
+% x = [pfxy(ii).x(1) pfxy(ii).x(end)];
+% y = [pfxy(ii).y(1) pfxy(ii).y(end)];
+%
+% if x(2)-x(1)>0
+% intang = atand((y(2)-y(1))/(x(2)-x(1)));
+% ang = 90 - intang;
+% elseif x(2)-x(1)<0
+% intang = atand((y(2)-y(1))/(x(2)-x(1)));
+% ang = 270 - intang;
+% end
+% TN(ii,1) = ang;
+% end
+
+%or else get information from site specific file
+addpath('J:\Coastal\')
+coastalinit
+
+% pfs = {'PF1','PF2','PF4','PF6','PF8'};
+% pfs = {'PF1','PF2','PF3'};
+pfs = {'PF1','PF2','PF3','PF4','PF5'};
+
+for ii = 1:length(pfs)
+ [~,rot_angle] = get_rotation_specs('DEEWH',pfs{ii});
+ alfa(ii,1) = -rot_angle*180/pi;
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/stormdates_narra.mat b/example/input_data/fielddata/stormdates_narra.mat
new file mode 100644
index 0000000..8009d09
Binary files /dev/null and b/example/input_data/fielddata/stormdates_narra.mat differ
diff --git a/example/input_data/fielddata/stormdates_other.mat b/example/input_data/fielddata/stormdates_other.mat
new file mode 100644
index 0000000..3a7a562
Binary files /dev/null and b/example/input_data/fielddata/stormdates_other.mat differ
diff --git a/example/input_data/fielddata/tide/code/tidedataget.m b/example/input_data/fielddata/tide/code/tidedataget.m
new file mode 100644
index 0000000..fba5416
--- /dev/null
+++ b/example/input_data/fielddata/tide/code/tidedataget.m
@@ -0,0 +1,44 @@
+% Written by Joshua Simmons 04/01/2016
+% Description: used to tranform offshore to nearshore mike21 using the code
+% OS2NS_MIKE21.m
+
+clear
+clc
+
+%% Settings
+%tide gauge data
+guagedataloc = 'J:\Coastal\Data\Tide\Sydney\Processed\SydTideData.mat';
+
+site = 'Narrabeen'; %embayment of interest
+
+%storm details
+stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_narra.mat';
+storm = {'2007', '2011', '2014', '2015'};
+
+% save?
+savebool = 1;
+% datasave location
+saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\tide\processed\';
+
+%% Code
+%load buoy data,profile data and lookup table
+guagedata = load(guagedataloc);
+
+stormdates = load(stormdatesloc);
+
+% narrow down to date range
+for ii = 1:length(storm)
+ startdate = stormdates.(['storm_' storm{ii}]).startdate;
+ enddate = stormdates.(['storm_' storm{ii}]).enddate;
+
+ validentries = find(guagedata.SydTideData.time>=startdate & guagedata.SydTideData.time<=enddate);
+
+ tidedata.tide = guagedata.SydTideData.tide(validentries);
+ tidedata.time = guagedata.SydTideData.time(validentries);
+
+ if savebool
+ save([saveloc 'tidedata_' lower(site) '_' storm{ii} '.mat'],'tidedata')
+ end
+
+ clear tidedata startdate enddate validentries
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/tide/code/tidedataget_other.m b/example/input_data/fielddata/tide/code/tidedataget_other.m
new file mode 100644
index 0000000..bc645a8
--- /dev/null
+++ b/example/input_data/fielddata/tide/code/tidedataget_other.m
@@ -0,0 +1,44 @@
+% Written by Joshua Simmons 04/01/2016
+% Description: used to tranform offshore to nearshore mike21 using the code
+% OS2NS_MIKE21.m
+
+clear
+clc
+
+%% Settings
+%tide gauge data
+guagedataloc = 'J:\Coastal\Data\Tide\Sydney\Processed\SydTideData.mat';
+
+site = 'Monavale'; %embayment of interest
+
+%storm details
+stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_other.mat';
+storm = {'2014', '2015'};
+
+% save?
+savebool = 1;
+% datasave location
+saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\tide\processed\';
+
+%% Code
+%load buoy data,profile data and lookup table
+guagedata = load(guagedataloc);
+
+stormdates = load(stormdatesloc);
+
+% narrow down to date range
+for ii = 1:length(storm)
+ startdate = stormdates.(['storm_' storm{ii}]).startdate;
+ enddate = stormdates.(['storm_' storm{ii}]).enddate;
+
+ validentries = find(guagedata.SydTideData.time>=startdate & guagedata.SydTideData.time<=enddate);
+
+ tidedata.tide = guagedata.SydTideData.tide(validentries);
+ tidedata.time = guagedata.SydTideData.time(validentries);
+
+ if savebool
+ save([saveloc 'tidedata_' lower(site) '_' storm{ii} '.mat'],'tidedata')
+ end
+
+ clear tidedata startdate enddate validentries
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/tide/processed/tidedata_narrabeen_2015.mat b/example/input_data/fielddata/tide/processed/tidedata_narrabeen_2015.mat
new file mode 100644
index 0000000..9e71718
Binary files /dev/null and b/example/input_data/fielddata/tide/processed/tidedata_narrabeen_2015.mat differ
diff --git a/example/input_data/fielddata/wave/code/fmikeos2ns.m b/example/input_data/fielddata/wave/code/fmikeos2ns.m
new file mode 100644
index 0000000..bf1042f
--- /dev/null
+++ b/example/input_data/fielddata/wave/code/fmikeos2ns.m
@@ -0,0 +1,105 @@
+function nswaves = fmikeos2ns(oswaves, LookUps, site, contour, location)
+ % Written by Joshua Simmons (WRL) 03/11/2014
+ % to transform waves from offshore (OS) to nearshore (NS) using MIKE21 model
+ % oswaves - OS wave data - oswaves.Hsig(Significant wave height)/
+ % oswaves.Tp1(Peak period)/oswaves.Wdir(Incident wave direction)
+ % LookUps - MIKE21 wave model output in file "OEH_ARC_PointOutput.mat"
+ % site - embayment of interest - i.e. 'Narrabeen'
+ % contour - 10 or 15 (m) contour for NS output
+ % location - either pure S co-ord (see Mitchell Harley's xyz2spz
+ % transformation) or location.E and location.N - refers to output
+ % location for NS waves (gets nearest MIKE21 point)
+ %
+ %V1.1 - updated by Joshua Simmons 21/01/2015 to patch Wamberal/Wanda
+ %error.
+
+ %find site and contour indices
+ siteidx = getnameidx(LookUps(:,2), site);
+
+ %apply fix for Wamberal and Wanda
+ % from TM on 21/01/2015: "I've found an error in the .mat file for the
+ % look-up tables; in the pdf I say that all output for all locations runs
+ % north to south - for Wamberal and Wanda it seems I've output them south to north.
+ % So if you are using output from these two sites, if you haven't realised already,
+ % alongshore gradients will look inverted."
+ for aa = 1:2;
+ LookUps{6,1}{aa,1} = flipud(LookUps{6,1}{aa,1}); LookUps{6,1}{aa,1}(:,1) = flipud(LookUps{6,1}{aa,1}(:,1));
+ LookUps{7,1}{aa,1} = flipud(LookUps{7,1}{aa,1}); LookUps{7,1}{aa,1}(:,1) = flipud(LookUps{7,1}{aa,1}(:,1));
+ end
+
+ if contour == 10
+ contouridx = 1;
+ elseif contour == 15;
+ contouridx = 2;
+ else
+ error('Contour information not available, select either 10 or 15m contour')
+ end
+
+ %extract scenarios for easy lookup
+ scenariosidx = getnameidx(LookUps(:,2), 'Offshore Wave Scenarios');
+ scenarios = LookUps{scenariosidx,1};
+
+ %get nearest location available for output
+ Eouts = LookUps{siteidx,1}{contouridx,1}(:,2);
+ Nouts = LookUps{siteidx,1}{contouridx,1}(:,3);
+
+ if strcmp(site,'Narrabeen')
+ %compares closes alongshore spacing (s) for log-spiral narrabeen
+ xyz_data.x = cell2mat(Eouts);
+ xyz_data.y = cell2mat(Nouts);
+ xyz_data.z = zeros(size(Eouts));
+ spz = xyz2spz(xyz_data,'NARRA');
+ [~,locationidx] = min(abs(location-spz.s));
+ else
+ %gets nearest output point
+ dists = sqrt((location.E-cell2mat(Eouts)).^2+(location.N-cell2mat(Nouts)).^2);
+ [~,locationidx] = min(abs(dists));
+ end
+
+ fcount = 1;
+ %loop through oswaves and lookup nswaves
+ for ii = 1:length(oswaves.Hsig)
+ tmpHs = oswaves.Hsig(ii);
+ tmpTp = oswaves.Tp1(ii);
+ tmpWdir = oswaves.Wdir(ii);
+
+ %find closest scenario match
+ [~,Hsidx] = min(abs(tmpHs-scenarios(:,1)));
+ tmpHsround = scenarios(Hsidx,1);
+ [~,Tpidx] = min(abs(tmpTp-scenarios(:,2)));
+ tmpTpround = scenarios(Tpidx,2);
+ [~,Wdiridx] = min(abs(tmpWdir-scenarios(:,3)));
+ tmpWdirround = scenarios(Wdiridx,3);
+
+ iter1ind = find(scenarios(:,3) == tmpWdirround);
+ iter2ind = find(scenarios(iter1ind,2) == tmpTpround);
+ iter3ind = find(scenarios(iter1ind(iter2ind),1) == tmpHsround);
+
+ iter4ind = iter1ind(iter2ind(iter3ind));
+
+ %now get ns lookup information if available
+ if isempty(iter4ind)
+ %get the closest wave height to scenario chosen
+ [~,iter3ind] = min(abs(tmpHs-scenarios(iter1ind(iter2ind),1)));
+ iter4ind = iter1ind(iter2ind(iter3ind));
+ diffHs = scenarios(iter4ind,1)-tmpHs;
+ disp(['No match found for oswaves(' num2str(ii) ') wave height off by ' num2str(diffHs) 'm']);
+ end
+
+ nstmp = LookUps{siteidx,1}{contouridx,1}{locationidx,4}(iter4ind,:);
+ %store the values from the lookup table
+ nswaves.Hsig(ii,1) = nstmp(1);
+ nswaves.Hmax(ii,1) = nstmp(2);
+ nswaves.Tp1(ii,1) = nstmp(3);
+ nswaves.T01(ii,1) = nstmp(4);
+ nswaves.Wdir(ii,1) = nstmp(5);
+ nswaves.power(ii,1) = nstmp(6);
+ end
+ try
+ nswaves.date = oswaves.date;
+ nswaves.time = oswaves.time;
+ nswaves.datenums = oswaves.datenums;
+ catch
+ nswaves.datenums = oswaves.datenums;
+ end
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21.m b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21.m
new file mode 100644
index 0000000..5540332
--- /dev/null
+++ b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21.m
@@ -0,0 +1,95 @@
+% Written by Joshua Simmons 04/01/2016
+% Description: used to tranform offshore to nearshore mike21 using the code
+% OS2NS_MIKE21.m
+
+clear
+clc
+
+%% Settings
+%buoy data
+buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
+% the non CSIRO filtered data is the only data available after 2014 so fill
+% from the recent data with gaps.
+backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
+
+%location of mike21 runcode
+mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
+%location of MALT
+maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
+%location of profiles for which wave data is required.
+profiledataloc = [437 850 1750 2606 3420]; % s co-ordinate
+
+site = 'Narrabeen'; %embayment of interest
+contour = 10; %m - 10m or 15m contour?
+
+%storm details
+stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_narra.mat';
+storm = {'2007', '2011', '2014', '2015'};
+
+% save?
+savebool = 1;
+% datasave location
+saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
+
+%% Code
+% for MIKE21 code and MALT
+addpath(mike21loc);
+addpath(maltloc);
+%load buoy data,profile data and lookup table
+buoydata = load(buoydataloc);
+backupbuoydata = load(backupbuoydata);
+load('OEH_ARC_PointOutput.mat');
+stormdates = load(stormdatesloc);
+
+% narrow down to date range
+for ii = 1:length(storm)
+ startdate = stormdates.(['storm_' storm{ii}]).startdate;
+ enddate = stormdates.(['storm_' storm{ii}]).enddate;
+
+ validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
+ % if nothing then revert to backup
+ if isempty(validentries) || validentries(end) == length(buoydata.dates)
+ validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
+ oswaves.Hsig = backupbuoydata.Hsig(validentries);
+ oswaves.Tp1 = backupbuoydata.Tp(validentries);
+ oswaves.Wdir = backupbuoydata.Wdir(validentries);
+ oswaves.datenums = backupbuoydata.dates(validentries);
+ else
+ oswaves.Hsig = buoydata.Hsig(validentries);
+ oswaves.Tp1 = buoydata.Tp(validentries);
+ oswaves.Wdir = buoydata.Wdir(validentries);
+ oswaves.datenums = buoydata.dates(validentries);
+ end
+
+ fns1 = fieldnames(oswaves);
+ for kk = 1:length(fns1)
+ if any(isnan(oswaves.(fns1{kk})))
+ error('NaN values in wave data')
+ end
+ end
+
+ tic
+ parfor jj = 1:length(profiledataloc)
+ location = profiledataloc(jj);
+ nswaves = OS2NS_MIKE21(oswaves, LookUps, site, contour, location);
+
+ % bug fix to ensure no inaccurate wave period estimates
+ nswaves.Tp1 = oswaves.Tp1;
+
+ tmp(jj) = nswaves;
+ end
+ wavedata.ns.data = tmp;
+ toc
+
+ wavedata.os = oswaves;
+
+ wavedata.ns.model = 'MIKE21';
+ wavedata.ns.contour = [num2str(contour) 'm'];
+ wavedata.ns.location = profiledataloc;
+
+ if savebool
+ save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
+ end
+
+ clear wavedata oswaves startdate enddate validentries tmp
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_deewhy.m b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_deewhy.m
new file mode 100644
index 0000000..17ca5d3
--- /dev/null
+++ b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_deewhy.m
@@ -0,0 +1,99 @@
+% Written by Joshua Simmons 04/01/2016
+% Description: used to tranform offshore to nearshore mike21 using the code
+% OS2NS_MIKE21.m
+
+clear
+clc
+
+%% Settings
+%buoy data
+buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
+% the non CSIRO filtered data is the only data available after 2014 so fill
+% from the recent data with gaps.
+backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
+
+%location of mike21 runcode
+mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
+%location of MALT
+maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
+%location of profiles for which wave data is required.
+% profiledataloc = [6275797.660794820 6275692.739810870 6275536.815570840]; % Northings co-ordinate - Bilgola (at 10m contour)
+% profiledataloc = [6272148.185 6272338.391 6272389.496 6271924.185 6271823.724]; % Northings co-ordinate - Mona Vale (at 10m contour)
+profiledataloc = [6264042.586 6264127.494 6263995.016 6263745.956 6263879.058]; % Northings co-ordinate - Dee Why (at 10m contour)
+profiledataloc2 = [343633.9263 343038.1544 342940.1201 342757.2994 342848.2807]; % Eastings co-ordinate - Dee Why (at 10m contour)
+
+site = 'DeeWhy'; %embayment of interest
+contour = 10; %m - 10m or 15m contour?
+
+%storm details
+stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_other.mat';
+storm = {'2014', '2015'};
+
+% save?
+savebool = 1;
+% datasave location
+saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
+
+%% Code
+% for MIKE21 code and MALT
+addpath(mike21loc);
+addpath(maltloc);
+%load buoy data,profile data and lookup table
+buoydata = load(buoydataloc);
+backupbuoydata = load(backupbuoydata);
+load('OEH_ARC_PointOutput.mat');
+stormdates = load(stormdatesloc);
+
+% narrow down to date range
+for ii = 1:length(storm)
+ startdate = stormdates.(['storm_' storm{ii}]).startdate;
+ enddate = stormdates.(['storm_' storm{ii}]).enddate;
+
+ validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
+ % if nothing then revert to backup
+ if isempty(validentries) || validentries(end) == length(buoydata.dates)
+ validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
+ oswaves.Hsig = backupbuoydata.Hsig(validentries);
+ oswaves.Tp1 = backupbuoydata.Tp(validentries);
+ oswaves.Wdir = backupbuoydata.Wdir(validentries);
+ oswaves.datenums = backupbuoydata.dates(validentries);
+ else
+ oswaves.Hsig = buoydata.Hsig(validentries);
+ oswaves.Tp1 = buoydata.Tp(validentries);
+ oswaves.Wdir = buoydata.Wdir(validentries);
+ oswaves.datenums = buoydata.dates(validentries);
+ end
+
+ fns1 = fieldnames(oswaves);
+ for kk = 1:length(fns1)
+ if any(isnan(oswaves.(fns1{kk})))
+ error('NaN values in wave data')
+ end
+ end
+
+ tic
+ for jj = 1:length(profiledataloc)
+ location.N = profiledataloc(jj);
+ location.E = profiledataloc2(jj);
+ nswaves = fmikeos2ns(oswaves, LookUps, site, contour, location);
+
+ % bug fix to ensure no inaccurate wave period estimates
+ nswaves.Tp1 = oswaves.Tp1;
+
+ tmp(jj) = nswaves;
+ end
+ wavedata.ns.data = tmp;
+ toc
+
+ wavedata.os = oswaves;
+
+ wavedata.ns.model = 'MIKE21';
+ wavedata.ns.contour = [num2str(contour) 'm'];
+ wavedata.ns.location = profiledataloc;
+
+ if savebool
+ save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
+ end
+
+ clear wavedata oswaves startdate enddate validentries tmp
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_end_effects_check.m b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_end_effects_check.m
new file mode 100644
index 0000000..5dae8fa
--- /dev/null
+++ b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_end_effects_check.m
@@ -0,0 +1,105 @@
+% Written by Joshua Simmons 04/01/2016
+% Description: used to tranform offshore to nearshore mike21 using the code
+% OS2NS_MIKE21.m
+
+clear
+clc
+
+%% Settings
+%buoy data
+buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
+% the non CSIRO filtered data is the only data available after 2014 so fill
+% from the recent data with gaps.
+backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
+
+%location of mike21 runcode
+mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
+%location of MALT
+maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
+%location of profiles for which wave data is required.
+profiledataloc = [90:20:3870]; % s co-ordinate
+
+site = 'Narrabeen'; %embayment of interest
+contour = 10; %m - 10m or 15m contour?
+
+%storm details
+stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_narra.mat';
+storm = {'2007', '2011', '2014', '2015'};
+
+% save?
+savebool = 0;
+% datasave location
+saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
+
+%% Code
+% for MIKE21 code and MALT
+addpath(mike21loc);
+addpath(maltloc);
+%load buoy data,profile data and lookup table
+buoydata = load(buoydataloc);
+backupbuoydata = load(backupbuoydata);
+load('OEH_ARC_PointOutput.mat');
+stormdates = load(stormdatesloc);
+
+% narrow down to date range
+for ii = 1:length(storm)
+ startdate = stormdates.(['storm_' storm{ii}]).startdate;
+ enddate = stormdates.(['storm_' storm{ii}]).enddate;
+
+ validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
+ % if nothing then revert to backup
+ if isempty(validentries) || validentries(end) == length(buoydata.dates)
+ validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
+ oswaves.Hsig = backupbuoydata.Hsig(validentries);
+ oswaves.Tp1 = backupbuoydata.Tp(validentries);
+ oswaves.Wdir = backupbuoydata.Wdir(validentries);
+ oswaves.datenums = backupbuoydata.dates(validentries);
+ else
+ oswaves.Hsig = buoydata.Hsig(validentries);
+ oswaves.Tp1 = buoydata.Tp(validentries);
+ oswaves.Wdir = buoydata.Wdir(validentries);
+ oswaves.datenums = buoydata.dates(validentries);
+ end
+
+ fns1 = fieldnames(oswaves);
+ for kk = 1:length(fns1)
+ if any(isnan(oswaves.(fns1{kk})))
+ error('NaN values in wave data')
+ end
+ end
+
+ tic
+ parfor jj = 1:length(profiledataloc)
+ location = profiledataloc(jj);
+ nswaves = OS2NS_MIKE21(oswaves, LookUps, site, contour, location);
+
+ % bug fix to ensure no inaccurate wave period estimates
+ nswaves.Tp1 = oswaves.Tp1;
+
+% tmp(jj) = nswaves;
+ Hsig(:,jj) = nswaves.Hsig(:,1);
+ Wdir(:,jj) = nswaves.Wdir(:,1);
+ power(:,jj) = nswaves.power(:,1);
+ end
+ dataout(ii).s = profiledataloc;
+ dataout(ii).Hsig = Hsig;
+ dataout(ii).Wdir = Wdir;
+ dataout(ii).power = power;
+ dataout(ii).osWdir = oswaves.Wdir;
+ dataout(ii).osHsig = oswaves.Hsig;
+% wavedata.ns.data = tmp;
+ toc
+
+ wavedata.os = oswaves;
+
+ wavedata.ns.model = 'MIKE21';
+ wavedata.ns.contour = [num2str(contour) 'm'];
+ wavedata.ns.location = profiledataloc;
+
+ if savebool
+ save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
+ end
+
+ clear wavedata oswaves startdate enddate validentries
+ clear power s Hsig Wdir
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_other.m b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_other.m
new file mode 100644
index 0000000..f3927e1
--- /dev/null
+++ b/example/input_data/fielddata/wave/code/run_OS_to_NS_mike21_other.m
@@ -0,0 +1,99 @@
+% Written by Joshua Simmons 04/01/2016
+% Description: used to tranform offshore to nearshore mike21 using the code
+% OS2NS_MIKE21.m
+
+clear
+clc
+
+%% Settings
+%buoy data
+buoydataloc = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed_CSIROGapFill.mat';
+% the non CSIRO filtered data is the only data available after 2014 so fill
+% from the recent data with gaps.
+backupbuoydata = 'J:\Coastal\Data\Wave\Sydney\Processed\SydneyProcessed.mat';
+
+%location of mike21 runcode
+mike21loc = 'J:\Coastal\Tools\Nearshore Wave Transformations\MIKE21\';
+%location of MALT
+maltloc = 'J:\Coastal\Tools\MALT Logspiral Transformation\';
+%location of profiles for which wave data is required.
+% profiledataloc = [6275797.660794820 6275692.739810870 6275536.815570840]; % Northings co-ordinate - Bilgola (at 10m contour)
+% profiledataloc = [6272148.185 6272338.391 6272389.496 6271924.185 6271823.724]; % Northings co-ordinate - Mona Vale (at 10m contour)
+profiledataloc = [6264042.586 6264127.494 6263995.016 6263745.956 6263879.058]; % Northings co-ordinate - Dee Why (at 10m contour)
+
+
+
+site = 'DeeWhy'; %embayment of interest
+contour = 10; %m - 10m or 15m contour?
+
+%storm details
+stormdatesloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\stormdates_other.mat';
+storm = {'2014', '2015'};
+
+% save?
+savebool = 1;
+% datasave location
+saveloc = 'D:\7 Erosion Models\a_nb_model_comparison\input_data\fielddata\wave\processed\';
+
+%% Code
+% for MIKE21 code and MALT
+addpath(mike21loc);
+addpath(maltloc);
+%load buoy data,profile data and lookup table
+buoydata = load(buoydataloc);
+backupbuoydata = load(backupbuoydata);
+load('OEH_ARC_PointOutput.mat');
+stormdates = load(stormdatesloc);
+
+% narrow down to date range
+for ii = 1:length(storm)
+ startdate = stormdates.(['storm_' storm{ii}]).startdate;
+ enddate = stormdates.(['storm_' storm{ii}]).enddate;
+
+ validentries = find(buoydata.dates>=startdate & buoydata.dates<=enddate);
+ % if nothing then revert to backup
+ if isempty(validentries) || validentries(end) == length(buoydata.dates)
+ validentries = find(backupbuoydata.dates>=startdate & backupbuoydata.dates<=enddate);
+ oswaves.Hsig = backupbuoydata.Hsig(validentries);
+ oswaves.Tp1 = backupbuoydata.Tp(validentries);
+ oswaves.Wdir = backupbuoydata.Wdir(validentries);
+ oswaves.datenums = backupbuoydata.dates(validentries);
+ else
+ oswaves.Hsig = buoydata.Hsig(validentries);
+ oswaves.Tp1 = buoydata.Tp(validentries);
+ oswaves.Wdir = buoydata.Wdir(validentries);
+ oswaves.datenums = buoydata.dates(validentries);
+ end
+
+ fns1 = fieldnames(oswaves);
+ for kk = 1:length(fns1)
+ if any(isnan(oswaves.(fns1{kk})))
+ error('NaN values in wave data')
+ end
+ end
+
+ tic
+ for jj = 1:length(profiledataloc)
+ location = profiledataloc(jj);
+ nswaves = OS2NS_MIKE21(oswaves, LookUps, site, contour, location);
+
+ % bug fix to ensure no inaccurate wave period estimates
+ nswaves.Tp1 = oswaves.Tp1;
+
+ tmp(jj) = nswaves;
+ end
+ wavedata.ns.data = tmp;
+ toc
+
+ wavedata.os = oswaves;
+
+ wavedata.ns.model = 'MIKE21';
+ wavedata.ns.contour = [num2str(contour) 'm'];
+ wavedata.ns.location = profiledataloc;
+
+ if savebool
+ save([saveloc 'wavedata_' lower(site) '_' storm{ii} '.mat'],'wavedata')
+ end
+
+ clear wavedata oswaves startdate enddate validentries tmp
+end
\ No newline at end of file
diff --git a/example/input_data/fielddata/wave/processed/wavedata_narrabeen_2015.mat b/example/input_data/fielddata/wave/processed/wavedata_narrabeen_2015.mat
new file mode 100644
index 0000000..15b3848
Binary files /dev/null and b/example/input_data/fielddata/wave/processed/wavedata_narrabeen_2015.mat differ
diff --git a/example/input_data/fielddata/wave/~$2007_survey_to_survey_wave_data.xlsx b/example/input_data/fielddata/wave/~$2007_survey_to_survey_wave_data.xlsx
new file mode 100644
index 0000000..fd0b6a9
Binary files /dev/null and b/example/input_data/fielddata/wave/~$2007_survey_to_survey_wave_data.xlsx differ
diff --git a/example/input_data/morphology/processed/morphology_narrabeen_2015.mat b/example/input_data/morphology/processed/morphology_narrabeen_2015.mat
new file mode 100644
index 0000000..235d6a9
Binary files /dev/null and b/example/input_data/morphology/processed/morphology_narrabeen_2015.mat differ
diff --git a/example/setup_info/profileinfo_narrabeen.mat b/example/setup_info/profileinfo_narrabeen.mat
new file mode 100644
index 0000000..9203935
Binary files /dev/null and b/example/setup_info/profileinfo_narrabeen.mat differ
diff --git a/formatpfdata.m b/formatpfdata.m
new file mode 100644
index 0000000..c31dff5
--- /dev/null
+++ b/formatpfdata.m
@@ -0,0 +1,21 @@
+function [pfdataout] = formatpfdata(xin,zin,interpspace)
+ % Written by Joshua Simmons 11/02/2015 - to format profile data as
+ % required (linear interploation etc etc
+ % Input:
+ % xin = cross shore
+ % zin = elevation
+ % interppspace = regular spacing of the profile to be linearly
+ % interpolated to. Leave as "NaN" if raw input points to
+ % be used.
+ % Output:
+ % structure with .x,.z formatted as required.
+ if isnan(interpspace)
+ pfdataout.x = xin;
+ pfdataout.z = zin;
+ else
+ minx = min(xin);
+ maxx = max(xin);
+ pfdataout.x = [minx:interpspace:maxx];
+ pfdataout.z = interp1(xin,zin,pfdataout.x);
+ end
+end
\ No newline at end of file
diff --git a/fsbeach_checktsdata.m b/fsbeach_checktsdata.m
new file mode 100644
index 0000000..4e538e6
--- /dev/null
+++ b/fsbeach_checktsdata.m
@@ -0,0 +1,36 @@
+% SBEACH toolbox - fsbeach_checktsdata
+% Written by Joshua Simmons 11/2017
+
+% returns all the SBEACH default values
+
+% Syntax:
+% OPT = fsbeach_checktsdata()
+
+% Input:
+% waves = wave data structure - output by fsbeach_genwaves
+% tides = tide data structure - output by fsbeach_gentides
+
+function fsbeach_checktsdata(waves,tide)
+ % check the tide and wave data matches
+ if tide.time(1) ~= waves.time(1)
+ error('Start times for wave and tide data dont match')
+ end
+ if tide.time(end) ~= waves.time(end)
+ error('End times for wave and tide data dont match. Make sure start/end times are hourly (for wave data).')
+ end
+
+ %check and allocate time parameters;
+ tt = round((tide.time - tide.time(1))*24*60);
+ %check time steps
+ dtt = diff(tt);
+ if any(dtt~=dtt(1))
+ error('Inconsistent time steps in tide data');
+ end
+
+ %now check the wave data lines up (os only required the rest if copied)
+ wt = round((waves.time - waves.time(1))*24*60);
+ dwt = diff(wt);
+ if any(dwt~=dwt(1))
+ error('Inconsistent time steps in wave data');
+ end
+end
\ No newline at end of file
diff --git a/fsbeach_collectresults.m b/fsbeach_collectresults.m
new file mode 100644
index 0000000..27fff60
--- /dev/null
+++ b/fsbeach_collectresults.m
@@ -0,0 +1,34 @@
+% SBEACH toolbox - fsbeach_collectresults
+% Written by Joshua Simmons 11/2017
+
+% Syntax:
+% sbout = fsbeach_collectresults(sblocs)
+%
+% Input:
+% sblocs = Cell array of locations where SBEACH has been run
+
+% Output:
+% sbout = structure with sbeach data;
+
+function sbout = fsbeach_collectresults(sblocs)
+ for ii = 1:length(sblocs)
+ fid = fopen([sblocs{ii} '.Excel.ASSUPPLIED'],'r');
+ formatSpec = '%f %f %f';
+ sizeA = [3 Inf];
+ tempresults = fscanf(fid,formatSpec,sizeA)';
+ fclose(fid);
+ clear fid
+
+ %get computed profile
+ sbout.computed(ii).z(:,1) = tempresults(:,3);
+ sbout.computed(ii).x = tempresults(:,1);
+
+ %input bathy
+ sbout.input(ii).x = tempresults(:,1);
+ sbout.input(ii).z = tempresults(:,2);
+
+ %now get wave elevations etc
+ fid = fopen([sblocs{ii} '.XVR'],'r');
+ sbout.xvr(ii) = fsbeach_readxvrfile(fid);
+ end
+end
\ No newline at end of file
diff --git a/fsbeach_genbathy.m b/fsbeach_genbathy.m
new file mode 100644
index 0000000..6227ae8
--- /dev/null
+++ b/fsbeach_genbathy.m
@@ -0,0 +1,103 @@
+% SBEACH toolbox - fsbeach_genbathy
+% Written by Joshua Simmons 11/2017
+
+% Syntax:
+% bathyout = fsbeach_genbathy(varargin)
+%
+% Input:
+% varargin = x: x-coordinates of bathymetry
+% z: z-coordinates of bathymetry
+% hbz: z-coordinates of hard bottom
+% ne: vector or matrix of the size of z containing
+% either booleans indicating if a cell is
+% non-erodable or a numeric value indicating the
+% thickness of the erodable layer on top of a
+% non-erodable layer
+%
+% Output:
+% bathyout = stucture of x and z
+
+function bathyout = fsbeach_genbathy(outpath,name,settings,varargin)
+ OPT = struct( ...
+ 'x', [0 100 200 300 400 500]', ...
+ 'z', [-10 -5 0 1 3 10]', ...
+ 'hbz', NaN, ...
+ 'ne', [] ...
+ );
+
+ OPT = setproperty(OPT, varargin{:});
+
+ %transpose them if needed
+ OPT = fsbeach_makecolumn(OPT,{'x','z','ne'});
+
+ %check they can be zipped
+ if length(OPT.x) ~= length(OPT.z)
+ error('x is of length %d, while z is of length %d...',length(OPT.x),length(OPT.z))
+ end
+
+ xin = OPT.x; %store x input
+
+ %check that constant grid is provided as required by SBEACH
+ diffx = diff(OPT.x);
+ dx = diffx(1);
+ if any(diffx~=dx)
+ dx = mode(diffx);
+ fprintf('\n')
+ warning('SBEACH requires constant spacing - interpolating to a regular grid with spacing of %d...',dx)
+ newx = [OPT.x(1):dx:OPT.x(end)]';
+ newz = interp1(OPT.x,OPT.z,newx);
+ OPT.x = newx;
+ OPT.z = newz;
+ clear newx newz
+ end
+
+ if settings.NGRID
+ tempdx = [0];
+ for ii = 1:settings.NGRID
+ tempdx = [tempdx; ones(settings.NDXV(ii),1).*settings.DXV(ii)];
+ end
+ newx = OPT.x(1) + cumsum(tempdx);
+
+ if newx(end) ~= OPT.x(end)
+ error(['Variable grid is wrong, please adjust DXV and NDXV so that'...
+ 'xstart = %d and xend = %d. Currently xend = %d.'],newx(1),OPT.x(end),newx(end))
+ end
+
+ newz = interp1(OPT.x,OPT.z,newx);
+ OPT.x = newx;
+ OPT.z = newz;
+ clear newx newz
+ end
+
+ if ~isnan(OPT.hbz)
+ OPT.hbz = interp1(xin,OPT.hbz,OPT.x);
+ end
+
+ %check they are in the right direction
+ if mean(diff(OPT.x)) < 0
+ OPT.x = flip(OPT.x);
+ OPT.z = flip(OPT.z);
+ if ~isnan(OPT.hbz)
+ OPT.hbz = flip(OPT.hbz);
+ end
+ end
+ if mean(diff(OPT.z)) > 0
+ fprintf('\n')
+ error('Please ensure x coordinate system is increasing offshore')
+ end
+
+ %zip
+ prifile = [OPT.x OPT.z]; %initial pf
+
+ %and write
+ writetofilepri(prifile, outpath, name)
+
+ if ~isnan(OPT.hbz)
+ hdbfile = [OPT.x OPT.hbz]; %hard bottom pf
+ writetofilehdb(hdbfile, outpath, name);
+ end
+
+ bathyout.x = OPT.x;
+ bathyout.z = OPT.z;
+ bathyout.hbz = OPT.hbz;
+end
\ No newline at end of file
diff --git a/fsbeach_gentide.m b/fsbeach_gentide.m
new file mode 100644
index 0000000..069b32c
--- /dev/null
+++ b/fsbeach_gentide.m
@@ -0,0 +1,33 @@
+% SBEACH toolbox - fsbeach_gentide
+% Written by Joshua Simmons 11/2017
+
+% Syntax:
+% tideout = fsbeach_gentide(varargin)
+%
+% Input:
+% varargin = time: array of measurement times for tide (minutes)
+% elev: array of waterlevels at seaward model border (m)
+%
+% Output:
+% tideout = stucture of time and elevation
+
+function tideout = fsbeach_gentide(outpath,name,varargin)
+ OPT = struct( ...
+ 'time', [0 100 200 300 400 500]', ...
+ 'elev', [1 1.2 1.1 1.5 1]' ...
+ );
+
+ OPT = setproperty(OPT, varargin{:});
+
+ %transpose them if needed
+ OPT = fsbeach_makecolumn(OPT,{'time','elev'});
+
+ %zip
+ elvfile = OPT.elev; %initial pf
+
+ %and write
+ writetofileelv(elvfile, outpath, name)
+
+ tideout.time = OPT.time;
+ tideout.elev = OPT.elev;
+end
\ No newline at end of file
diff --git a/fsbeach_genwaves.m b/fsbeach_genwaves.m
new file mode 100644
index 0000000..96f63ef
--- /dev/null
+++ b/fsbeach_genwaves.m
@@ -0,0 +1,62 @@
+% SBEACH toolbox - fsbeach_genwaves
+% Written by Joshua Simmons 11/2017
+
+% Syntax:
+% wavesout = fsbeach_genwaves(varargin)
+%
+% Input:
+% varargin = time: array of measurement times for w (minutes)
+%
+% Hsig: significant wave height (m)
+% Tp1: peak wave period (s)
+% ang: wave direction (degrees from True North)
+% pfang: Profile orientation - from the shore to the
+% most seaward point (degrees from True North)
+%
+% Output:
+% wavesout = structure array of wave inputs
+
+function wavesout = fsbeach_genwaves(outpath,name,varargin)
+ OPT = struct( ...
+ 'time', [0 100 200 300 400 500]', ...
+ 'Hsig', [0 100 200 300 400 500]', ...
+ 'Tp1', [-10 -5 0 1 3 10]', ...
+ 'ang', [-10 -5 0 1 3 10]', ...
+ 'pfang', NaN ...
+ );
+
+ OPT = setproperty(OPT, varargin{:});
+
+ %transpose them if needed
+ OPT = fsbeach_makecolumn(OPT,{'time','Hsig','Tp1','ang'});
+
+ %check they can be zipped
+ if length(OPT.Hsig) ~= length(OPT.Tp1)
+ error('Hsig is of length %d, while Tp1 is of length %d...',length(OPT.Hsig),length(OPT.Tp1))
+ end
+
+ %zip wave data
+ wavfile = [OPT.Hsig OPT.Tp1];
+
+ % correct for SBEACH wave angle convention
+ % See SBEACH Report 3 - User's manual (pg 19)
+ if isnan(OPT.pfang)
+ error('Please provide the orientation of your profile as pfang.')
+ end
+
+ angfile = OPT.pfang - OPT.ang;
+
+ %check
+ if any(abs(angfile)>90)
+ warning('Check the orientation of your profile relative to the waves.')
+ end
+
+ %and write
+ writetofilewav(wavfile, outpath, name)
+ writetofileang(angfile, outpath, name)
+
+ wavesout.time = OPT.time;
+ wavesout.Hsig = OPT.Hsig;
+ wavesout.Tp1 = OPT.Tp1;
+ wavesout.ang = angfile;
+end
\ No newline at end of file
diff --git a/fsbeach_getdefaults.m b/fsbeach_getdefaults.m
new file mode 100644
index 0000000..65bfc29
--- /dev/null
+++ b/fsbeach_getdefaults.m
@@ -0,0 +1,109 @@
+% SBEACH toolbox - fsbeach_getdefaults
+% Written by Joshua Simmons 11/2017
+
+% returns all the SBEACH default values
+
+% Syntax:
+% OPT = fsbeach_getdefaults()
+
+function OPT = fsbeach_getdefaults()
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % A - MODEL SETUP
+ OPT.TITLE = 'Test';
+ OPT.UNITS = 1;
+ OPT.NDX = NaN;
+ OPT.XSTART = NaN;
+ OPT.IDX = 0;
+ OPT.DXC = 1;
+ OPT.NGRID = 0;
+ OPT.DXV = 0;
+ OPT.NDXV = 0;
+ OPT.NDT = NaN;
+ OPT.DT = 5;
+ OPT.NWR = 0;
+ OPT.WRI = 0;
+ OPT.ICOMP = 0;
+ OPT.ELV1 = 0;
+ OPT.ELV2 = 0.25;
+ OPT.ELV3 = 0.5;
+ OPT.EDP1 = 0.5;
+ OPT.EDP2 = 1;
+ OPT.EDP3 = 1.5;
+ OPT.REFELV = 0;
+
+ % Calibration parameters
+ OPT.K = 1.75e-06;
+ OPT.EPS = 0.002;
+ OPT.LAMM = 0.5;
+
+ OPT.TEMPC = 20;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % B - WAVES/WATER ELEVATION/WIND
+ OPT.WVTYPE = 2;
+ OPT.IWAVE = 1;
+ OPT.HIN = 0;
+ OPT.T = 0;
+ OPT.DTWAV = NaN;
+ OPT.IANG = 1;
+ OPT.ZIN = 0;
+ OPT.DTANG = NaN;
+ OPT.DMEAS = NaN; %0 is deepwater
+ OPT.IRAND = 0;
+ OPT.ISEED = 4567;
+ OPT.RPERC = 5;
+ OPT.IELEV = 1;
+ OPT.TELEV = 0;
+ OPT.DTELV = NaN;
+ OPT.IWIND = 0;
+ OPT.W = 0;
+ OPT.ZWIND = 0;
+ OPT.DTWND = 60;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % C - BEACH
+ OPT.TPIN = 1;
+ OPT.SCHEM = 1;
+ OPT.XLAND = 0;
+ OPT.DLAND = [];
+ OPT.XLBDUNE = [];
+ OPT.DLBDUNE = [];
+ OPT.XLCDUNE = [];
+ OPT.DLCDUNE = [];
+ OPT.XSCDUNE = [];
+ OPT.DSCDUNE = [];
+ OPT.XBERMS = [];
+ OPT.DBERMS = [];
+ OPT.XBERME = [];
+ OPT.DBERME = [];
+ OPT.XFORS = [];
+ OPT.DFORS = [];
+
+ %Calibration params
+ OPT.DFS = 0.5;
+ OPT.D50 = 0.35;
+ OPT.BMAX = 17.2;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % D - BEACH FILL
+ OPT.IBCHFILL = 0;
+ OPT.XBFS = 0;
+ OPT.XBFE = 0;
+ OPT.NFILL = 0;
+ OPT.XF = [];
+ OPT.EFILL = [];
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % E - SEAWALL/REVETMENT
+ OPT.ISWALL = 0;
+ OPT.XSWALL = 0;
+ OPT.ISWFAIL = 0;
+ OPT.PEFAIL = 0;
+ OPT.WEFAIL = 0;
+ OPT.HFAIL = 0;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % F - HARD BOTTOM
+ OPT.IHBOT = 0;
+ OPT.SCF = 1;
+end
\ No newline at end of file
diff --git a/fsbeach_getstdsettings.m b/fsbeach_getstdsettings.m
new file mode 100644
index 0000000..a30ad81
--- /dev/null
+++ b/fsbeach_getstdsettings.m
@@ -0,0 +1,55 @@
+% SBEACH toolbox - fsbeach_getstdsettings
+% Written by Joshua Simmons 11/2017
+
+% returns all the SBEACH default values
+
+% Syntax:
+% OPT = fsbeach_getstdsettings()
+
+% Input:
+% bathy = bathy data structure - output by fsbeach_genbathy
+% waves = wave data structure - output by fsbeach_genwaves
+% tides = tide data structure - output by fsbeach_gentides
+%
+% Output:
+% OPT = settings
+
+function OPT = fsbeach_getstdsettings(bathy, waves, tide, OPT)
+ % preliminary work
+ diffx = diff(bathy.x);
+ dx = diffx(1);
+
+ wavestt = (waves.time(end)-waves.time(1))*24*60;
+ dwt = round((waves.time(2)-waves.time(1))*24*60);
+ if isstruct(tide)
+ tidett = (tide.time(end)-tide.time(1))*24*60; %total time
+ dtt = round((tide.time(2)-tide.time(1))*24*60); %time steps
+ end
+
+ % Section A
+ OPT.DXC = dx;
+ OPT.NDX = length(bathy.x); %number of cells
+ if isnan(OPT.XSTART)
+ OPT.XSTART = bathy.x(1); %landward calculation boundary
+ else
+ warning('XSTART is %d m from first profile point',OPT.XSTART - bathy.x(1))
+ end
+
+ OPT.NDT = round(wavestt/OPT.DT);
+
+ % Section B
+ OPT.DTWAV = dwt;
+ OPT.DTANG = dwt;
+ if isnan(OPT.DMEAS)
+ OPT.DMEAS = -min(bathy.z); %assume it is input at the boundary if not specified
+ end
+
+ if isstruct(tide)
+ OPT.DTELV = dtt;
+ end
+
+ % Section F - check for hard bottom
+ if ~isnan(bathy.hbz)
+ OPT.IHBOT = 1;
+ end
+end
\ No newline at end of file
diff --git a/fsbeach_makecolumn.m b/fsbeach_makecolumn.m
new file mode 100644
index 0000000..fb2e141
--- /dev/null
+++ b/fsbeach_makecolumn.m
@@ -0,0 +1,18 @@
+% SBEACH toolbox - fsbeach_makecolumn
+% Written by Joshua Simmons 11/2017
+
+% convert to column vectors
+% Syntax:
+% OPT = fsbeach_makecolumn(OPT,vars)
+%
+% Input:
+% OPT: toolbox OPT structure
+% vars: variable names to be checked in structure
+
+function OPT = fsbeach_makecolumn(OPT,vars)
+ for ii = 1:length(vars)
+ if ~iscolumn(OPT.(vars{ii}))
+ OPT.(vars{ii}) = OPT.(vars{ii})';
+ end
+ end
+end
\ No newline at end of file
diff --git a/fsbeach_modelskill.m b/fsbeach_modelskill.m
new file mode 100644
index 0000000..1f051c3
--- /dev/null
+++ b/fsbeach_modelskill.m
@@ -0,0 +1,211 @@
+function [bss, rmse, dVc, dVm, ddV] = fsbeach_modelskill(measured, computed, initial)
+%edited by Joshua Simmons 04/04/2016 to revise definitions of r2 etc
+%XB_SKILL Computes a variety skill scores
+%
+% Computes a variety skill scores: R^2, Sci, Relative bias, Brier Skill
+% Score. Special feature: within the XBeach testbed, the results are
+% stored to be able to show the development of the different skill scores
+% in time.
+%
+% Syntax:
+% [r2 sci relbias bss] = xb_skill(measured, computed, varargin)
+%
+% Input:
+% measured = Measured data where the first column contains independent
+% values and the second column contains dependent values
+% computed = Computed data where the first column contains independent
+% values and the second column contains dependent values
+% initial = Initial data where the first column contains independent
+% values and the second column contains dependent values
+% varargin = var: Name of the variable that is supplied
+%
+% Output:
+% bss = Brier Skill Score
+% rb20 = Regression through origin R2 value
+%
+% Example
+% [r2 sci relbias bss] = xb_skill(measured, computed)
+% [r2 sci relbias bss] = xb_skill(measured, computed, 'var', 'zb')
+%
+% See also xb_plot_skill
+
+%% Copyright notice
+% --------------------------------------------------------------------
+% Copyright (C) 2011 Deltares
+% Bas Hoonhout
+%
+% bas.hoonhout@deltares.nl
+%
+% Rotterdamseweg 185
+% 2629HD Delft
+%
+% This library is free software: you can redistribute it and/or
+% modify it under the terms of the GNU Lesser General Public
+% License as published by the Free Software Foundation, either
+% version 2.1 of the License, or (at your option) any later version.
+%
+% This library is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+% Lesser General Public License for more details.
+%
+% You should have received a copy of the GNU Lesser General Public
+% License along with this library. If not, see <http://www.gnu.org/licenses/>.
+% --------------------------------------------------------------------
+
+% This tool is part of <a href="http://OpenEarth.nl">OpenEarthTools</a>.
+% OpenEarthTools is an online collaboration to share and manage data and
+% programming tools in an open source, version controlled environment.
+% Sign up to recieve regular updates of this function, and to contribute
+% your own tools.
+
+%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
+% Created: 13 Apr 2011
+% Created with Matlab version: 7.9.0.529 (R2009b)
+
+% $Id: xb_skill.m 4725 2011-06-27 12:58:48Z hoonhout $
+% $Date: 2011-06-27 14:58:48 +0200 (ma, 27 jun 2011) $
+% $Author: hoonhout $
+% $Revision: 4725 $
+% $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/applications/xbeach/xb_analysis/xb_skill.m $
+% $Keywords: $
+
+%% read options
+
+if ~exist('initial', 'var'); initial = []; end;
+
+%% xbeach bugfix
+% bug - sometime xbeach has suprious points on the end (a flat end to
+% the profile) this is fake so remove it.
+if length(computed)>=10
+ if nanmean(computed(1:10,2)) < nanmean(computed(end-9:end,2))
+ %xbeach direction
+ spr_len = length(diff(computed(:,2))) - find(diff(computed(:,2))~=0,1,'last');
+ computed(end-spr_len+1:end,2) = NaN;
+ else
+ spr_len = find(diff(computed(:,2))~=0,1,'first');
+ computed(1:spr_len-1,2) = NaN;
+ end
+ if spr_len > 2
+ disp(['Warning check computed points, ' num2str(spr_len) ' points found to be spurious'] )
+ end
+end
+
+%% remove nans
+
+measured = measured(~any(isnan(measured), 2),:);
+computed = computed(~any(isnan(computed), 2),:);
+
+if ~isempty(initial)
+ initial = initial (~any(isnan(initial ), 2),:);
+end
+
+%% compute skills
+
+[bss rmse dVc dVm ddV] = deal(nan);
+
+if size(measured,1) > 1 && size(computed,1) > 1
+
+ x = unique([computed(:,1) ; measured(:,1)]);
+ x = x(~isnan(x));
+
+ zmt = interp1(measured(:,1), measured(:,2), x);
+ zct = interp1(computed(:,1), computed(:,2), x);
+
+ % determine active zone
+ if ~isempty(initial)
+ zit = interp1(initial(:,1), initial(:,2), x);
+
+ %average over 10 steps for small spikes
+ dzmi = smooth(abs(zmt-zit),10);
+ dzmi((isnan(zmt)|isnan(zit))) = 0;
+ dzcm = smooth(abs(zct-zmt),10);
+ dzcm((isnan(zct)|isnan(zmt))) = 0;
+
+ dz = max(dzmi,dzcm);
+ %has to be more change than a standard measurement error (10cm say)
+ zi1 = find(dz>0.1,1,'first');
+ zi2 = find(dz>0.1,1,'last');
+
+ x = x (zi1:zi2);
+ zmt = zmt(zi1:zi2);
+ zct = zct(zi1:zi2);
+ zit = zit(zi1:zi2);
+ end
+
+ %% clear all nan values in input (to avoid problems with interpolation)
+ %all equal at this point and assign values for BSS calc
+ x0 = x;
+ xc = x;
+ xm = x;
+ z0 = zit;
+ zc = zct;
+ zm = zmt;
+
+ x0(isnan(z0)) = [];
+ z0(isnan(z0)) = [];
+ xc(isnan(zc)) = [];
+ zc(isnan(zc)) = [];
+ xm(isnan(zm)) = [];
+ zm(isnan(zm)) = [];
+
+ %% determine new grid covered by all profiles
+ x_new = unique([x0' xc' xm']);
+
+ %% interpolate profiles onto new grid
+ [x0,id1] = unique(x0);
+ [xc,id2] = unique(xc);
+ [xm,id3] = unique(xm);
+ z0 = z0(id1);
+ zc = zc(id2);
+ zm = zm(id3);
+
+ z0_new = interp1(x0, z0, x_new);
+ zc_new = interp1(xc, zc, x_new);
+ zm_new = interp1(xm, zm, x_new);
+
+ % reduce vectors if any nans are left
+ id = ~(isnan(z0_new) | isnan(zc_new) | isnan(zm_new));
+ [x_new z0_new zc_new zm_new] = deal(x_new(id), z0_new(id), zc_new(id), zm_new(id));
+
+ %% Calculate volume change
+ try
+ Vc = volumefind(x_new,zc_new,0);
+ Vm = volumefind(x_new,zm_new,0);
+ V0 = volumefind(x_new,z0_new,0);
+
+ dVm = Vm-V0; %measured volume change
+ dVc = Vc-V0; %computed volume change
+ ddV = dVc-dVm; %predicted volume change - measured
+ nddV = ddV/abs(V0-Vm);
+ catch
+ dVm = NaN;
+ dVc = NaN;
+ ddV = NaN;
+ nddV = NaN;
+ end
+ %% derive weight
+ % weight is defined as half of diff left and half of diff right of each
+ % point; first and last point are considered to have only one side.
+ weight = diff([x_new(1) diff(x_new)/2 + x_new(1:end-1) x_new(end)]);
+ weight = weight / sum(weight);
+
+ %% calculate BSS only
+ % brier skill score
+ if ~isempty(initial)
+ mse_p = sum( (zm_new - zc_new).^2 .* weight );
+ mse_0 = sum( (zm_new - z0_new).^2 .* weight );
+ bss = 1. - (mse_p/mse_0);
+ else
+ bss = 1-(std(zc-zm))^2/var(zm);
+ end
+
+ %% Calculate others
+
+ rmse = sqrt(mean((zc_new-zm_new).^2));
+ rmsm = sqrt(mean(zm_new.^2));
+ sci = rmse/max(rmsm,abs(mean(zm_new)));
+ relbias = mean(zc_new-zm_new)/max(rmsm,abs(mean(zm_new)));
+
+end
+
diff --git a/fsbeach_processresults.m b/fsbeach_processresults.m
new file mode 100644
index 0000000..08457e9
--- /dev/null
+++ b/fsbeach_processresults.m
@@ -0,0 +1,71 @@
+% SBEACH toolbox - fsbeach_processresults
+% Written by Joshua Simmons 11/2017
+
+% Syntax:
+% sbout = fsbeach_processresults(sblocs)
+%
+% Input:
+% sbout = structure with input, computed and measured
+
+% Output:
+% sbout = structure with sbeach data and statistics;
+
+function sbout = fsbeach_processresults(sbout)
+ for ii = 1:length(sbout.input)
+ xin = sbout.input(ii).x;
+ zin = sbout.input(ii).z;
+ xcomp = sbout.computed(ii).x;
+ zcomp = sbout.computed(ii).z;
+ xmeas = sbout.measured(ii).x;
+ zmeas = sbout.measured(ii).z;
+
+ try
+ [tmpbss, tmprmse, dVc, dVm, ddV] = fsbeach_modelskill([xmeas zmeas], [xcomp zcomp], [xin zin]);
+
+ %calculate BSS above MSL
+ %first get points above MSL - for measured and then apply that
+ %to all - cant take these independently as may only get say one
+ %point for modelled and then xb_skill interpolates.
+ %BSS deals well if based on measured as the -10 and x
+ %determinant - it only finds the interpolated values with
+ %info so will deal if initpf doesnt stretch to all of
+ %measured pf
+ keepme1 = find(zmeas>0);
+ if xmeas(keepme1(1))>xmeas(keepme1(end))
+ keepmo1 = find(xcomp>xmeas(keepme1(end)) & xcomp<xmeas(keepme1(1)));
+ keepin1 = find(xin>xmeas(keepme1(end)) & xin<xmeas(keepme1(1)));
+ elseif xmeas(keepme1(1))<xmeas(keepme1(end))
+ keepmo1 = find(xcomp<xmeas(keepme1(end)) & xcomp>xmeas(keepme1(1)));
+ keepin1 = find(xin<xmeas(keepme1(end)) & xin>xmeas(keepme1(1)));
+ end
+ [tmpabMSLbss, tmpabMSLrmse, ~, ~, ~] = fsbeach_modelskill([xmeas(keepme1) zmeas(keepme1)], [xcomp(keepmo1) zcomp(keepmo1)], [xin(keepin1) zin(keepin1)]);
+
+ %calculate BSS below MSL
+ keepme2 = find(zmeas<0);
+ if xmeas(keepme2(1))>xmeas(keepme2(end))
+ keepmo2 = find(xcomp>xmeas(keepme2(end)) & xcomp<xmeas(keepme2(1)));
+ keepin2 = find(xin>xmeas(keepme2(end)) & xin<xmeas(keepme2(1)));
+ elseif xmeas(keepme2(1))<xmeas(keepme2(end))
+ keepmo2 = find(xcomp<xmeas(keepme2(end)) & xcomp>xmeas(keepme2(1)));
+ keepin2 = find(xin<xmeas(keepme2(end)) & xin>xmeas(keepme2(1)));
+ end
+ [tmpbeMSLbss,tmpbeMSLrmse, ~, ~, ~] = fsbeach_modelskill([xmeas(keepme2) zmeas(keepme2)], [xcomp(keepmo2) zcomp(keepmo2)], [xin(keepin2) zin(keepin2)]);
+
+ %store results
+ sbout.stats(ii).bss = tmpbss;
+ sbout.stats(ii).rmse = tmprmse;
+ sbout.stats(ii).dVc = dVc;
+ sbout.stats(ii).dVm = dVm;
+ sbout.stats(ii).ddV = ddV;
+
+ sbout.stats(ii).aboveMSL.bss = tmpabMSLbss;
+ sbout.stats(ii).aboveMSL.rmse = tmpabMSLrmse;
+
+ sbout.stats(ii).belowMSL.bss = tmpbeMSLbss;
+ sbout.stats(ii).belowMSL.rmse = tmpbeMSLrmse;
+ catch
+ disp(['Run number ' num2str(ii) ' failed.'])
+ continue
+ end
+ end
+end
diff --git a/fsbeach_readxvrfile.m b/fsbeach_readxvrfile.m
new file mode 100644
index 0000000..79acae2
--- /dev/null
+++ b/fsbeach_readxvrfile.m
@@ -0,0 +1,43 @@
+% SBEACH toolbox - fsbeach_collectresults
+% Written by Joshua Simmons 02/2018
+
+% Syntax:
+% dataout = fsbeach_readxvrfile(fid)
+%
+% Input:
+% fid = file id of .XVR file from which data is to be extracted
+
+% Output:
+% dataout = structure with sbeach data;
+
+function dataout = fsbeach_readxvrfile(fid)
+ ln = fgets(fid); % first line
+ %header material
+ while ln(1) ~= '_'
+ ln = fgets(fid);
+ end
+
+ for ii = 1:8
+ tmpd = [];
+ while any(isletter(ln)) || ln(1) == '_'
+ ln = fgets(fid);
+ end
+ while ln(1) ~= '_' && ~any(isletter(ln)) && ischar(ln)
+ tmpd = [tmpd; sscanf(ln,'%f')];
+ ln = fgets(fid);
+ end
+ tmpdataout(ii).data = tmpd;
+ end
+
+ padlen = length(tmpdataout(1).data)-length(tmpdataout(2).data);
+
+ %max value
+ dataout.Hmax = padarray(tmpdataout(2).data,padlen,0,'pre');
+ dataout.TWLmax = padarray(tmpdataout(4).data,padlen,0,'pre');
+ dataout.depmax = padarray(tmpdataout(6).data,padlen,0,'pre');
+
+ %time step of max
+ dataout.tsHmax = padarray(tmpdataout(3).data,padlen,0,'pre');
+ dataout.tsTWLmax = padarray(tmpdataout(5).data,padlen,0,'pre');
+ dataout.tsdepmax = padarray(tmpdataout(7).data,padlen,0,'pre');
+end
\ No newline at end of file
diff --git a/fsbeach_runmodel.m b/fsbeach_runmodel.m
new file mode 100644
index 0000000..892d3b4
--- /dev/null
+++ b/fsbeach_runmodel.m
@@ -0,0 +1,18 @@
+% SBEACH toolbox - fsbeach_runmodel
+% Written by Joshua Simmons 11/2017
+
+% Syntax:
+% success = fsbeach_runmodel(cfgpath,exepath)
+%
+% Input:
+% cfgpath = path to SBEACH .CFG file (and other files setup by
+% fsbeach_setupmodel()
+% exepath = path to SBEACH exe
+
+% Output:
+% success = bool;
+
+function [fail, cmdout] = fsbeach_runmodel(cfgpath,exepath)
+ cmdin = ['"' exepath '"' ' -cfgFN ' '"' cfgpath '"'];
+ [fail, cmdout] = system(cmdin);
+end
\ No newline at end of file
diff --git a/fsbeach_runsetup.asv b/fsbeach_runsetup.asv
new file mode 100644
index 0000000..315765d
--- /dev/null
+++ b/fsbeach_runsetup.asv
@@ -0,0 +1,6 @@
+% SBEACH toolbox
+% Written by Joshua Simmons
+
+function out = fsbeach_runsetup()
+
+end
\ No newline at end of file
diff --git a/fsbeach_runsetup.m b/fsbeach_runsetup.m
new file mode 100644
index 0000000..fcca4c7
--- /dev/null
+++ b/fsbeach_runsetup.m
@@ -0,0 +1,6 @@
+% SBEACH toolbox
+% Written by Joshua Simmons
+
+function out = fsbeach_runsetup()
+
+end
\ No newline at end of file
diff --git a/fsbeach_setupmodel.m b/fsbeach_setupmodel.m
new file mode 100644
index 0000000..72f13a8
--- /dev/null
+++ b/fsbeach_setupmodel.m
@@ -0,0 +1,96 @@
+% SBEACH toolbox - fsbeach_setupmodel
+% Written by Joshua Simmons 11/2017
+
+% Note this code is specifically written to align with (and borrows code
+% from) the XBeach matlab toolbox so as to make the use of both models
+% streamlined.
+
+% Syntax:
+% OPT = fsbeach_setupmodel(varargin)
+%
+% Input:
+% varargin = bathy: cell array of name/value pairs of bathymetry
+% settings supplied to fsbeach_genbathy
+% waves: cell array of name/value pairs of waves
+% settings supplied to fsbeach_genwaves
+% tide: cell array of name/value pairs of tide
+% settings supplied to fsbeach_gentide
+% settings: cell array of name/value pairs of model
+% settings supplied to fsbeach_writecfg
+% path: destination directory of model setup, if
+% written to disk
+%
+% Output:
+% OPT = options
+
+function out = fsbeach_setupmodel(varargin)
+ OPT = struct( ...
+ 'bathy', {{}}, ...
+ 'waves', {{}}, ...
+ 'tide', {{}}, ...
+ 'settings', {{}}, ...
+ 'path', [pwd '\temp\'], ...
+ 'name', 'sbeach_test' ...
+ );
+
+ OPT = setproperty(OPT, varargin{:});
+
+ %check we are ready to write!
+ if ~isdir(OPT.path)
+ mkdir(OPT.path)
+ else
+ %delete any old SBEACH files
+ delete(fullfile(OPT.path,'*.PRI'));
+ delete(fullfile(OPT.path,'*.WAV'));
+ delete(fullfile(OPT.path,'*.ANG'));
+ delete(fullfile(OPT.path,'*.ELV'));
+ delete(fullfile(OPT.path,'*.HDB'));
+ delete(fullfile(OPT.path,'*.CFG'));
+ delete(fullfile(OPT.path,'*.AsSupplied'));
+ delete(fullfile(OPT.path,'*.LOG'));
+ delete(fullfile(OPT.path,'*.PRC'));
+ delete(fullfile(OPT.path,'*.RPTNew'));
+ delete(fullfile(OPT.path,'*.XVR'));
+ end
+
+ %%% Initial settings
+ %get defaults
+ settingsOPT = fsbeach_getdefaults();
+ %apply any setttings as needed
+ settingsOPT = setproperty(settingsOPT, OPT.settings);
+
+ %%% Bathy
+ %print to file
+ bathy = fsbeach_genbathy(OPT.path,OPT.name,settingsOPT,OPT.bathy{:});
+
+ %%% Waves
+ %print to file
+ waves = fsbeach_genwaves(OPT.path,OPT.name,OPT.waves{:});
+
+ %%% Tide
+ %print to file
+ if isempty(OPT.tide) && ~any(strcmp(OPT.settings,'IELEV'))
+ error('Please supply tide or use the IELEV setting to supply constant WL.')
+ elseif ~isempty(OPT.tide)
+ tide = fsbeach_gentide(OPT.path,OPT.name,OPT.tide{:});
+ %check timeseries data
+ fsbeach_checktsdata(waves,tide);
+ else
+ tide = NaN; % no tide supplied
+ end
+
+ %%% Final settings
+ % get standard settings from bathy, waves and tides
+ settingsOPT = fsbeach_getstdsettings(bathy, waves, tide, settingsOPT);
+
+ %write the configuration file
+ fsbeach_writecfg(OPT.path,OPT.name,settingsOPT);
+
+ %%% Output
+ %structure out
+ out.bathy = bathy;
+ out.waves = waves;
+ out.tides = tide;
+ out.settings = settingsOPT;
+ out.path = [OPT.path OPT.name];
+end
\ No newline at end of file
diff --git a/fsbeach_writecfg.m b/fsbeach_writecfg.m
new file mode 100644
index 0000000..19cc855
--- /dev/null
+++ b/fsbeach_writecfg.m
@@ -0,0 +1,18 @@
+% SBEACH toolbox - fsbeach_writecfg
+% Written by Joshua Simmons 11/2017
+
+% Syntax:
+% wavesout = fsbeach_writecfg(outpath,standard_settings,varargin)
+%
+% Input:
+% outpath: output path for CFG file
+% standard_settings: standard settings for bathy, waves, tide etc.
+%
+% varargin = sbeach settings
+%
+% Output:
+% cfg = all config settings
+
+function fsbeach_writecfg(outpath,name,OPT)
+ writetofilecfg(OPT, outpath, name);
+end
\ No newline at end of file
diff --git a/sbeach_exe/final_exe/SBEACH.exe b/sbeach_exe/final_exe/SBEACH.exe
new file mode 100644
index 0000000..c6eac25
Binary files /dev/null and b/sbeach_exe/final_exe/SBEACH.exe differ
diff --git a/setproperty.m b/setproperty.m
new file mode 100644
index 0000000..c4ec1bb
--- /dev/null
+++ b/setproperty.m
@@ -0,0 +1,471 @@
+function [OPT Set Default] = setproperty(OPT, inputCell, varargin)
+% SETPROPERTY generic routine to set values in PropertyName-PropertyValue pairs
+%
+% Routine to set properties based on PropertyName-PropertyValue
+% pairs (aka <keyword,value> pairs). Can be used in any function
+% where PropertyName-PropertyValue pairs are used.
+%
+%
+% [OPT Set Default] = setproperty(OPT, inputCell, varargin)
+%
+% input:
+% OPT = structure in which fieldnames are the keywords and the values
+% are the defaults
+% inputCell = must be a cell array containing single struct, or a set of
+% 'PropertyName', PropertyValue,... pairs. It is best practice
+% to pass the varargin of the calling function as the 2nd argument.
+% Property names must be strings. If they contain a dot (.) this is
+% interpreted a field separator. This can be used to assign
+% properties on a subfield level.
+% varargin = series of 'PropertyName', PropertyValue,... pairs to set
+% methods for setproperty itself.
+%
+% output:
+% OPT = structure, similar to the input argument OPT, with possibly
+% changed values in the fields
+% Set = structure, similar to OPT, values are true where OPT has been
+% set (and possibly changed)
+% Default = structure, similar to OPT, values are true where the values of
+% OPT are equal to the original OPT
+%
+% Example calls:
+%
+% [OPT Set Default] = setproperty(OPT, vararginOfCallingFunction)
+% [OPT Set Default] = setproperty(OPT, {'PropertyName', PropertyValue,...})
+% [OPT Set Default] = setproperty(OPT, {OPT2})
+%
+% Any number of leading structs with
+% any number of trailing <keyword,value> pairs:
+% [OPT Set Default] = setproperty(OPT1, OPT2, ..., OPTn)
+% [OPT Set Default] = setproperty(OPT1, OPT2, ..., OPTn,'PropertyName', PropertyValue,...)
+%
+% Different methods for dealing with class changes of variables, or
+% extra fields (properties that are not in the input structure) can
+% be defined as property-value pairs. Valid properties are
+% onExtraField and onClassChange:
+%
+% PROPERTY VALUE
+% onClassChange: ignore ignore (default)
+% warn throw warning
+% error throw error
+% onExtraField: silentIgnore silently ignore the field
+% warnIgnore ignore the field and throw warning
+% silentAppend silently append the field to OPT
+% warnAppend append the field to OPT and throw warning
+% error throw error (default)
+%
+% Example calls:
+%
+% [OPT Set Default] = setproperty(OPT, vararginOfCallingFunction,'onClassChange','warn')
+% [OPT Set Default] = setproperty(OPT, {'PropertyName', PropertyValue},'onExtraField','silentIgnore')
+%
+%
+% Example:
+%
+% +------------------------------------------->
+% function y = dosomething(x,'debug',1)
+% OPT.debug = 0;
+% OPT = setproperty(OPT, varargin);
+% y = x.^2;
+% if OPT.debug; plot(x,y);pause; end
+% +------------------------------------------->
+%
+% legacy syntax is also supported, but using legacy syntax prohibits the
+% setting of onClassChange and onExtraField methods:
+%
+% [OPT Set Default] = setproperty(OPT, varargin{:})
+% OPT = setproperty(OPT, 'PropertyName', PropertyValue,...)
+% OPT = setproperty(OPT, OPT2)
+%
+% input:
+% OPT = structure in which fieldnames are the keywords and the values are the defaults
+% varargin = series of PropertyName-PropertyValue pairs to set
+% OPT2 = is a structure with the same fields as OPT.
+%
+% Internally setproperty translates OPT2 into a set of
+% PropertyName-PropertyValue pairs (see example below) as in:
+% OPT2 = struct( 'propertyName1', 1,...
+% 'propertyName2', 2);
+% varcell = reshape([fieldnames(OPT2)'; struct2cell(OPT2)'], 1, 2*length(fieldnames(OPT2)));
+% OPT = setproperty(OPT, varcell{:});
+%
+% Change log:
+% 2011-09-30: full code rewrite to include:
+% - setpropertyInDeeperStruct functionality
+% - user defined handling of extra fields
+% - class change warning/error message
+%
+% See also: VARARGIN, STRUCT, MERGESTRUCTS
+
+%% Copyright notice
+% --------------------------------------------------------------------
+% Copyright (C) 2009 Delft University of Technology
+% C.(Kees) den Heijer
+%
+% C.denHeijer@TUDelft.nl
+%
+% Faculty of Civil Engineering and Geosciences
+% P.O. Box 5048
+% 2600 GA Delft
+% The Netherlands
+%
+% This library is free software; you can redistribute it and/or
+% modify it under the terms of the GNU Lesser General Public
+% License as published by the Free Software Foundation; either
+% version 2.1 of the License, or (at your option) any later version.
+%
+% This library is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+% Lesser General Public License for more details.
+%
+% You should have received a copy of the GNU Lesser General Public
+% License along with this library; if not, write to the Free Software
+% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+% USA
+% or http://www.gnu.org/licenses/licenses.html, http://www.gnu.org/, http://www.fsf.org/
+% --------------------------------------------------------------------
+
+% This tools is part of <a href="http://OpenEarth.Deltares.nl">OpenEarthTools</a>.
+% OpenEarthTools is an online collaboration to share and manage data and
+% programming tools in an open source, version controlled environment.
+% Sign up to recieve regular updates of this function, and to contribute
+% your own tools.
+
+%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
+% Created: 26 Feb 2009
+% Created with Matlab version: 7.4.0.287 (R2007a)
+
+% $Id: setproperty.m 8036 2013-02-05 21:06:43Z boer_g $
+% $Date: 2013-02-05 22:06:43 +0100 (Tue, 05 Feb 2013) $
+% $Author: boer_g $
+% $Revision: 8036 $
+% $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/matlab/setproperty.m $
+% $Keywords: $
+
+%% shortcut function if there is nothing to set (1)
+if nargin==1||isempty(inputCell)
+ if nargout > 1 % process Set
+ flds1 = fieldnames(OPT);
+ Set = [flds1,repmat({false},size(flds1))]';
+ Set = struct(Set{:});
+ end
+ if nargout > 2 % process Default
+ Default = [flds1,repmat({true},size(flds1))]';
+ Default = struct(Default{:});
+ end
+ return
+end
+
+%% check and parse inputCell (usually the varargin struct in the calling function)
+% determine mode from class of the second argument
+switch(class(inputCell))
+ case 'struct'
+ % recursively let setproperty peel all leading structs
+ % with optional trailing keyword,value> pairs
+ inputCell = setproperty(inputCell,varargin{:});
+ inputCell = struct2arg(inputCell);
+ varargin = {};
+ case 'char'
+ % legay syntax mode
+ inputCell = [{inputCell} varargin];
+ varargin = {};
+ case 'cell'
+ % recursively let setproperty peel all leading structs from cell
+ if isstruct(inputCell{1})
+ inputCell = {setproperty(inputCell{:})};
+ end
+ otherwise
+ error('SETPROPERTY:inputCell',...
+ 'Second input must be a cell, (or a char or struct for legacy syntax)')
+end
+
+
+%% shortcut function if there is nothing to set (2)
+if numel(inputCell) == 1 && isempty(inputCell{1})
+ if nargout > 1 % process Set
+ flds1 = fieldnames(OPT);
+ Set = [flds1,repmat({false},size(flds1))]';
+ Set = struct(Set{:});
+ end
+ if nargout > 2 % process Default
+ Default = [flds1,repmat({true},size(flds1))]';
+ Default = struct(Default{:});
+ end
+ return
+end
+
+if odd(length(inputCell))
+ %then length must be 1
+ if length(inputCell) == 1
+ %then the inputCell must be a structure
+ if ~isstruct(inputCell{1})
+ error('SETPROPERTY:inputCell',...
+ 'Second argument inputCell must be a cell containing a single struct or property/value pairs')
+ end
+ flds2 = fieldnames(inputCell{1})';
+ newVal = struct2cell(inputCell{1})';
+ else
+ error('SETPROPERTY:inputCell',...
+ 'Second argument inputCell must be a cell containing a single struct or property/value pairs')
+ end
+else
+ flds2 = inputCell(1:2:end);
+ if ~all(cellfun(@ischar,flds2))
+ error('SETPROPERTY:inputCell',...
+ 'Second argument inputCell, property names should be strings')
+ end
+ newVal = inputCell(2:2:end);
+end
+
+
+%% check and parse varargin
+
+if length(varargin) > 4
+ error('Too many input arguments to setpropoerty.');
+end
+% error(nargchk(0, 4, length(varargin), 'struct'))
+if odd(length(varargin))
+ error('SETPROPERTY:varargin',...
+ 'Set onClassChange and onExtraField with the varargin, as keyword value pairs');
+end
+
+onExtraField = 'error';
+onClassChange = 'ignore';
+if ~isempty(varargin)
+ for ii = 1:2:length(varargin)
+ switch varargin{ii}
+ case 'onClassChange'
+ if ischar(varargin{ii+1}) && ismember(varargin(ii+1),...
+ {'ignore','warn','error'})
+ onClassChange = varargin{ii+1};
+ else
+ error('SETPROPERTY:InvalidMethod',...
+ 'Valid methods for onClassChange are: ignore, warn and error');
+ end
+ case 'onExtraField'
+ if ischar(varargin{ii+1}) && ismember(varargin(ii+1),...
+ {'silentAppend','warnAppend','silentIgnore','warnIgnore','error'})
+ onExtraField = varargin{ii+1};
+ else
+ error('SETPROPERTY:InvalidMethod',...
+ 'Valid methods for onExtraField are: append, silentIgnore, warnIgnore, error');
+ end
+ otherwise
+ error('SETPROPERTY:InvalidKeyword', ...
+ 'Supported keywords are onClassChange and onExtraField');
+ end
+ end
+end
+
+switch onClassChange
+ case 'warn'
+ warnOnClassChange = true;
+ errorOnClassChange = false;
+ case 'error'
+ warnOnClassChange = false;
+ errorOnClassChange = true;
+ case 'ignore'
+ warnOnClassChange = false;
+ errorOnClassChange = false;
+end
+
+%% copy the original OPT only if Default has to be returned
+if nargout > 2
+ OPToriginal = OPT;
+end
+%% check field names, and distinguish between field to set and extra fields
+% flds1 contains field names of OPT
+flds1 = fieldnames(OPT);
+
+% fldsCell contains field names to assign, split per subfield
+fldsCell = regexp(flds2,'[^\.]*','match');
+fldsCellLen = cellfun(@length,fldsCell);
+% check if all field names are valid to assign
+allFieldNames = [fldsCell{:}];
+validFieldName = cellfun(@(x) (isvarname(x) | iskeyword(x)),allFieldNames); % a field like 'case' is allowed also, but is not a valid varname
+if ~all(validFieldName)
+ error('SETPROPERTY:invalidFieldName',...
+ '\nThe following field name is not valid: %s',allFieldNames{~validFieldName});
+end
+
+% identify fldsToSet for simple field names (without subfields)
+
+% fldsToSet = ismember(flds2,flds1);
+% As ismember is rather slow, this is replaced by more optimized code
+
+fldsToSet = false(size(flds2));
+for ii = find(fldsCellLen==1)
+ fldsToSet(ii) = any(strcmp(flds2(ii),flds1));
+end
+
+% identify fldsToSet for field names with subfields)
+for ii = find(fldsCellLen>1)
+ fldsToSet(ii) = isfield2(OPT,fldsCell{ii});
+end
+
+
+% all other fields are either to be set, or extra
+fldsExtra = ~fldsToSet;
+
+%% process fldsToSet
+if any(fldsToSet)
+ for ii = find(fldsToSet)
+ switch length(fldsCell{ii})
+ case 1
+ if warnOnClassChange || errorOnClassChange
+ class1 = class(OPT.(fldsCell{ii}{1}));
+ class2 = class(newVal{ii});
+ classChange(flds2{ii},class1,class2,warnOnClassChange,errorOnClassChange);
+ end
+ OPT.(fldsCell{ii}{1}) = newVal{ii};
+ case 2
+ if warnOnClassChange || errorOnClassChange
+ class1 = class(OPT.(fldsCell{ii}{1}).(fldsCell{ii}{2}));
+ class2 = class(newVal{ii});
+ classChange(flds2{ii},class1,class2,warnOnClassChange,errorOnClassChange);
+ end
+ OPT.(fldsCell{ii}{1}).(fldsCell{ii}{2}) = newVal{ii};
+ otherwise
+ if warnOnClassChange || errorOnClassChange
+ class1 = class(getfield(OPT,fldsCell{ii}{:}));
+ class2 = class(newVal{ii});
+ classChange(flds2{ii},class1,class2,warnOnClassChange,errorOnClassChange);
+ end
+ OPT = setfield(OPT,fldsCell{ii}{:},newVal{ii});
+ end
+ end
+end
+
+%% deal with fldsExtra
+if any(fldsExtra)
+ switch lower(onExtraField)
+ case {'silentappend','warnappend'}
+ % append extra fields to OPT
+ for ii = find(fldsExtra)
+ switch length(fldsCell{ii})
+ case 1
+ OPT.(flds2{ii}) = newVal{ii};
+ otherwise
+ OPT = setfield(OPT,fldsCell{ii}{:},newVal{ii});
+ end
+ end
+ if strcmpi(onExtraField,'warnAppend')
+ % throw a warning
+ warning('SETPROPERTY:ExtraField',...
+ '\nThe following field is Extra and appended to OPT: %s',flds2{fldsExtra});
+ end
+ case 'silentignore'
+ % do nothing, silently ignore extra fields
+ case 'warnignore'
+ warning('SETPROPERTY:ExtraField',...
+ '\nThe following field is Extra and thus ignored: %s',flds2{fldsExtra});
+ case 'error'
+ error('SETPROPERTY:ExtraField',...
+ '\nThe following field is Extra: %s',flds2{fldsExtra});
+ end
+end
+
+%% assign Output
+if nargout > 1 % process Set
+ Set = [flds1,repmat({false},size(flds1))]';
+ Set = struct(Set{:});
+ if any(fldsToSet)
+ for ii = find(fldsToSet)
+ switch length(fldsCell{ii})
+ case 1
+ Set.(fldsCell{ii}{1}) = true;
+ otherwise % this routine is rather slow, would be nice to add case 2 to cover the most common calls
+ for nn = 1:length(fldsCell{ii})-1
+ if ~isstruct(getfield(Set,fldsCell{ii}{1:nn})) %#ok<*GFLD>
+ substruct = [fieldnames(getfield(OPT,fldsCell{ii}{1:nn})),...
+ repmat({false},size(fieldnames(getfield(OPT,fldsCell{ii}{1:nn}))))]';
+ Set = setfield(Set,fldsCell{ii}{1:nn},[]); %#ok<*SFLD>
+ Set = setfield(Set,fldsCell{ii}{1:nn},struct(substruct{:}));
+ end
+ end
+ Set = setfield(Set,fldsCell{ii}{:},true);
+ end
+ end
+ end
+ if any(fldsExtra)
+ switch lower(onExtraField)
+ case {'silentappend','warnappend'}
+ for ii = find(fldsExtra)
+ switch length(fldsCell{ii})
+ case 1
+ Set.(flds2{ii}) = true;
+ otherwise
+ Set = setfield(Set,fldsCell{ii}{:},true);
+ end
+ end
+ end
+ end
+end
+
+if nargout > 2 % process Default
+ Default = [flds1,repmat({true},size(flds1))]';
+ Default = struct(Default{:});
+ if any(fldsToSet)
+ for ii = find(fldsToSet)
+ switch length(fldsCell{ii})
+ case 1
+ Default.(fldsCell{ii}{1}) = ...
+ isequalwithequalnans(OPToriginal.(fldsCell{ii}{1}),OPT.(fldsCell{ii}{1}));
+ otherwise % this routine is rather slow, would be nice to add case 2 to cover the most common calls
+ for nn = 1:length(fldsCell{ii})-1
+ if ~isstruct(getfield(Default,fldsCell{ii}{1:nn})) %#ok<*GFLD>
+ substruct = [fieldnames(getfield(OPT,fldsCell{ii}{1:nn})),...
+ repmat({true},size(fieldnames(getfield(OPT,fldsCell{ii}{1:nn}))))]';
+ Default = setfield(Default,fldsCell{ii}{1:nn},[]); %#ok<*SFLD>
+ Default = setfield(Default,fldsCell{ii}{1:nn},struct(substruct{:}));
+ end
+ end
+ Default = setfield(Default,fldsCell{ii}{:},...
+ isequalwithequalnans(getfield(OPToriginal,fldsCell{ii}{:}),getfield(OPT,fldsCell{ii}{:})));
+ end
+ end
+ end
+ if any(fldsExtra)
+ switch lower(onExtraField)
+ case {'silentappend','warnappend'}
+ for ii = find(fldsExtra)
+ switch length(fldsCell{ii})
+ case 1
+ Default.(flds2{ii}) = false;
+ otherwise
+ Default = setfield(Default,fldsCell{ii}{:},false);
+ end
+ end
+ end
+ end
+end
+
+function tf = isfield2(OPT,fldsCell)
+if isfield(OPT,fldsCell{1})
+ if length(fldsCell) > 1
+ tf = isfield2(OPT.(fldsCell{1}),fldsCell(2:end));
+ else
+ tf = true;
+ end
+else
+ tf = false;
+end
+
+function classChange(fld,class1,class2,warnOnClassChange,errorOnClassChange)
+if ~strcmp(class1,class2)
+ if warnOnClassChange
+ warning('SETPROPERTY:ClassChange', ...
+ 'Class change of field ''%s'' from %s to %s',...
+ fld,class1,class2);
+ elseif errorOnClassChange
+ error('SETPROPERTY:ClassChange', ...
+ 'Class change of field ''%s'' from %s to %s not allowed',...
+ fld,class1,class2);
+ end
+end
+
+function out = odd(in)
+%ODD test whether number if odd
+out = mod(in,2)==1;
+%% EOF
\ No newline at end of file
diff --git a/volumefind.m b/volumefind.m
new file mode 100644
index 0000000..bf3f1ff
--- /dev/null
+++ b/volumefind.m
@@ -0,0 +1,39 @@
+function [VolumeTotal, Beachwidth] = volumefind(x,z,shore)
+ % Calculates Volume and Beach Width
+ % Written by Matthew Phillips 2015
+
+ % Define Shoreline Contour and Fixed Landward Column Boundary
+ ShoreCont = shore; % m AHD
+ BackCol = 1; % column number
+
+ if ~iscolumn(z)
+ z = z';
+ end
+ if ~iscolumn(x)
+ x = x';
+ end
+
+ if z(end) > z(1)
+ z = flip(z);
+ x = flip(x);
+ end
+
+ if x(end) < x(1)
+ x = -x;
+ end
+
+ xstep = mean(diff(x));
+ shorecol = find(z>=ShoreCont,1,'last')+1;
+
+ % Width at MSL
+ Beachwidth = interp1(z(shorecol-1:shorecol,1),x(shorecol-1:shorecol,1),ShoreCont);
+
+ % Total Volume
+ Volx = [x(BackCol:shorecol-1); Beachwidth];
+ Volz = [z(BackCol:shorecol-1); ShoreCont];
+ VolumeTotal = trapz(Volx,Volz);
+end
+
+
+
+
diff --git a/writetofileang.m b/writetofileang.m
new file mode 100644
index 0000000..d8323c4
--- /dev/null
+++ b/writetofileang.m
@@ -0,0 +1,18 @@
+function writetofileang(angfile, outpath, pfstring)
+ fid = fopen([outpath '\' pfstring '.ANG'],'w');
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C> * SBEACH wave input file: %s.ANG *\r\n',pfstring);
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C> DIRECTION (DEGREES)\r\n');
+ fprintf(fid,'E>----------------------------------------------------------------\r\n');
+ for ii = 1:length(angfile)
+ if ii == length(angfile)
+ fprintf(fid,'%.3f\r\n', angfile(ii));
+ else
+ fprintf(fid,'%.3f\r\n', angfile(ii));
+ end
+ end
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/writetofilecfg.m b/writetofilecfg.m
new file mode 100644
index 0000000..620bda5
--- /dev/null
+++ b/writetofilecfg.m
@@ -0,0 +1,173 @@
+function writetofilecfg(OPT, outpath, pfstring)
+ fid = fopen([outpath '\' pfstring '.CFG'],'w');
+ fprintf(fid,' *******************************************************\r\n');
+ fprintf(fid,' * SBEACH model configuration file: %s.CFG *\r\n',pfstring);
+ fprintf(fid,' *******************************************************\r\n');
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % A - MODEL SETUP
+ fprintf(fid,'A---------------------------- MODEL SETUP -------------------------------A\r\n');
+ fprintf(fid,'A.1 RUN TITLE: TITLE\r\n');
+ fprintf(fid,' %s\r\n',pfstring);
+ fprintf(fid,'A.2 INPUT UNITS (SI=1, AMERICAN CUST.=2): UNITS\r\n');
+ fprintf(fid,' %d\r\n',OPT.UNITS);
+ fprintf(fid,'A.3 TOTAL NUMBER OF CALCULATION CELLS AND POSITION OF LANDWARD BOUNDARY\r\n');
+ fprintf(fid,' RELATIVE TO INITIAL PROFILE: NDX, XSTART\r\n');
+ fprintf(fid,' %d %d\r\n',OPT.NDX,OPT.XSTART);
+ fprintf(fid,'A.4 GRID TYPE (CONSTANT=0, VARIABLE=1): IDX\r\n');
+ fprintf(fid,' %d\r\n',OPT.IDX);
+ fprintf(fid,'A.5 COMMENT: IF GRID TYPE IS VARIABLE, CONTINUE TO A.8\r\n');
+ fprintf(fid,'A.6 CONSTANT GRID CELL WIDTH: DXC\r\n');
+ fprintf(fid,' %d\r\n',OPT.DXC);
+ fprintf(fid,'A.7 COMMENT: IF GRID TYPE IS CONSTANT CONTINUE TO A.10\r\n');
+ fprintf(fid,'A.8 NUMBER OF DIFFERENT GRID CELL REGIONS: NGRID\r\n');
+ fprintf(fid,' %d\r\n',OPT.NGRID);
+ fprintf(fid,'A.9 GRID CELL WIDTHS AND NUMBER OF CELLS IN EACH REGION FROM LANDWARD\r\n');
+ fprintf(fid,' TO SEAWARD BOUNDARY: (DXV(I), NDXV(I), I=1,NGRID)\r\n');
+ for ii = 1:length(OPT.NGRID)
+ fprintf(fid,' %d %d',OPT.DXV(ii),OPT.NDXV(ii));
+ end
+ fprintf(fid,'\r\n');
+ fprintf(fid,'A.10 NUMBER OF TIME STEPS AND VALUE OF TIME STEP IN MINUTES: NDT,DT\r\n');
+ fprintf(fid,' %d %d\r\n',OPT.NDT,OPT.DT);
+ fprintf(fid,'A.11 NUMBER OF TIME STEP(S) INTERMEDIATE OUTPUT IS WANTED: NWR\r\n');
+ fprintf(fid,' %d\r\n',OPT.NWR);
+ fprintf(fid,'A.12 TIME STEPS OF INTERMEDIATE OUTPUT: (WRI(I), I=1,NWR)\r\n');
+ for ii = 1:length(OPT.NWR)
+ fprintf(fid,' %d',OPT.WRI(ii));
+ end
+ fprintf(fid,'\r\n');
+ fprintf(fid,'A.13 IS A MEASURED PROFILE AVAILABLE FOR COMPARISON? (NO=0, YES=1): ICOMP\r\n');
+ fprintf(fid,' %d\r\n',OPT.ICOMP);
+ fprintf(fid,'A.14 THREE PROFILE ELEVATION CONTOURS (MAXIMUM HORIZONTAL RECESSION OF EACH\r\n');
+ fprintf(fid,' WILL BE DETERMINED): ELV1, ELV2, ELV3\r\n');
+ fprintf(fid,' %.2f %.2f %.2f\r\n',OPT.ELV1,OPT.ELV2,OPT.ELV3);
+ fprintf(fid,'A.15 THREE PROFILE EROSION DEPTHS AND REFERENCE ELEVATION (DISTANCE FROM\r\n');
+ fprintf(fid,' POSITION OF REFERENCE ELEVATION ON INITIAL PROFILE TO POSITION OF\r\n');
+ fprintf(fid,' LANDWARD MOST OCCURENCE OF EACH EROSION DEPTH WILL BE DETERMINED\r\n');
+ fprintf(fid,' EDP1, EDP2, EDP3, REFELV\r\n');
+ fprintf(fid,' %.2f %.2f %.2f %.2f\r\n',OPT.EDP1,OPT.EDP2,OPT.EDP3,OPT.REFELV);
+ fprintf(fid,'A.16 TRANSPORT RATE COEFFICIENT (m^4/N): K\r\n');
+ fprintf(fid,' %.3e\r\n',OPT.K);
+ fprintf(fid,'A.17 COEFFICIENT FOR SLOPE-DEPENDENT TERM (m^2/s): EPS\r\n');
+ fprintf(fid,' %.3f\r\n',OPT.EPS);
+ fprintf(fid,'A.18 TRANSPORT RATE DECAY COEFFICIENT MULTIPLIER: LAMM\r\n');
+ fprintf(fid,' %.3f\r\n',OPT.LAMM);
+ fprintf(fid,'A.19 WATER TEMPERATURE IN DEGREES C: TEMPC\r\n');
+ fprintf(fid,' %.2f\r\n',OPT.TEMPC);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % B - WAVES/WATER ELEVATION/WIND
+ fprintf(fid,'B---------------------- WAVES/WATER ELEVATION/WIND -----------------------B\r\n');
+ fprintf(fid,'B.1 WAVE TYPE (MONOCHROMATIC=1, IRREGULAR=2): WVTYPE\r\n');
+ fprintf(fid,' %d\r\n',OPT.WVTYPE);
+ fprintf(fid,'B.2 WAVE HEIGHT AND PERIOD INPUT (CONSTANT=0, VARIABLE=1): IWAVE\r\n');
+ fprintf(fid,' %d\r\n',OPT.IWAVE);
+ fprintf(fid,'B.3 COMMENT: IF WAVE HEIGHT AND PERIOD ARE VARIABLE, CONTINUE TO B.6\r\n');
+ fprintf(fid,'B.4 CONSTANT WAVE HEIGHT AND PERIOD: HIN, T\r\n');
+ fprintf(fid,' %.3f %.3f\r\n',OPT.HIN,OPT.T);
+ fprintf(fid,'B.5 COMMENT: IF WAVE HEIGHT AND PERIOD ARE CONSTANT, CONTINUE TO B.7\r\n');
+ fprintf(fid,'B.6 TIME STEP OF VARIABLE WAVE HEIGHT AND PERIOD INPUT IN MINUTES: DTWAV\r\n');
+ fprintf(fid,' %.2f\r\n',OPT.DTWAV);
+ fprintf(fid,'B.7 WAVE ANGLE INPUT (CONSTANT=0, VARIABLE=1): IANG\r\n');
+ fprintf(fid,' %d\r\n',OPT.IANG);
+ fprintf(fid,'B.8 COMMENT: IF WAVE ANGLE IS VARIABLE, CONTINUE TO B.11\r\n');
+ fprintf(fid,'B.9 CONSTANT WAVE ANGLE: ZIN\r\n');
+ fprintf(fid,' %.3f\r\n',OPT.ZIN);
+ fprintf(fid,'B.10 COMMENT: IF WAVE ANGLE IS CONSTANT, CONTINUE TO B.12\r\n');
+ fprintf(fid,'B.11 TIME STEP OF VARIABLE WAVE ANGLE INPUT IN MINUTES: DTANG\r\n');
+ fprintf(fid,' %.2f\r\n',OPT.DTANG);
+ fprintf(fid,'B.12 WATER DEPTH OF INPUT WAVES (DEEPWATER=0): DMEAS\r\n');
+ fprintf(fid,' %d\r\n',OPT.DMEAS);
+ fprintf(fid,'B.13 IS RANDOMIZATION OF WAVE HEIGHT DESIRED? (NO=0, YES=1): IRAND\r\n');
+ fprintf(fid,' %d\r\n',OPT.IRAND);
+ fprintf(fid,'B.14 COMMENT: IF RANDOMIZATION OF WAVE HEIGHT IS NOT DESIRED, CONTINUE TO B.16\r\n');
+ fprintf(fid,'B.15 SEED VALUE FOR RANDOMIZER AND PERCENT OF VARIABILITY: ISEED, RPERC\r\n');
+ fprintf(fid,' %d %.2f\r\n',OPT.ISEED,OPT.RPERC);
+ fprintf(fid,'B.16 TOTAL WATER ELEVATION INPUT (CONSTANT=0, VARIABLE=1): IELEV\r\n');
+ fprintf(fid,' %d\r\n',OPT.IELEV);
+ fprintf(fid,'B.17 COMMENT: IF WATER ELEVATION IS VARIABLE CONTINUE TO B.20\r\n');
+ fprintf(fid,'B.18 CONSTANT TOTAL WATER ELEVATION: TELEV\r\n');
+ fprintf(fid,' %.3f\r\n',OPT.TELEV);
+ fprintf(fid,'B.19 COMMENT: IF WATER ELEVATION IS CONSTANT, CONTINUE TO B.21\r\n');
+ fprintf(fid,'B.20 TIME STEP OF VARIABLE TOTAL WATER ELEVATION INPUT IN MINUTES: DTELV\r\n');
+ fprintf(fid,' %.2f\r\n',OPT.DTELV);
+ fprintf(fid,'B.21 WIND SPEED AND ANGLE INPUT (CONSTANT=0, VARIABLE=1): IWIND\r\n');
+ fprintf(fid,' %d\r\n',OPT.IWIND);
+ fprintf(fid,'B.22 COMMENT: IF WIND SPEED AND ANGLE ARE VARIABLE, CONTINUE TO B.25\r\n');
+ fprintf(fid,'B.23 CONSTANT WIND SPEED AND ANGLE: W,ZWIND\r\n');
+ fprintf(fid,' %.3f %.3f\r\n',OPT.W,OPT.ZWIND);
+ fprintf(fid,'B.24 COMMENT: IF WIND SPEED AND ANGLE ARE CONSTANT, CONTINUE TO C.\r\n');
+ fprintf(fid,'B.25 TIME STEP OF VARIABLE WIND SPEED AND ANGLE INPUT IN MINUTES: DTWND\r\n');
+ fprintf(fid,' %.2f\r\n',OPT.DTWND);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % C - BEACH
+ fprintf(fid,'C------------------------------- BEACH ----------------------------------C\r\n');
+ fprintf(fid,'C.1 TYPE OF INPUT PROFILE (ARBITRARY=1, SCHEMATIZED=2): TPIN\r\n');
+ fprintf(fid,' %d\r\n',OPT.TPIN);
+ fprintf(fid,'C.2 COMMENT: IF PROFILE TYPE IS ARBITRARY CONTINUE TO C.4\r\n');
+ fprintf(fid,'C.3 LOCATION AND ELEVATION OF LANDWARD BOUNDARY, LANDWARD BASE OF DUNE,\r\n');
+ fprintf(fid,' LANDWARD CREST OF DUNE, SEAWARD CREST OF DUNE, START OF BERM,\r\n');
+ fprintf(fid,' END OF BERM, AND FORESHORE: XLAND,DLAND,XLBDUNE,DLBDUNE,XLCDUNE,DLCDUNE,\r\n');
+ fprintf(fid,' XSCDUNE,DSCDUNE,XBERMS,DBERMS,XBERME,DBERME,XFORS,DFORS\r\n');
+ fprintf(fid,' %d %d %d %d %d %d\r\n',OPT.XLAND,OPT.DLAND,OPT.XLBDUNE,OPT.DLBDUNE,OPT.XLCDUNE,OPT.DLCDUNE);
+ fprintf(fid,' %d %d %d %d %d %d %d %d\r\n',OPT.XSCDUNE,OPT.DSCDUNE,OPT.XBERMS,OPT.DBERMS,OPT.XBERME,OPT.DBERME,OPT.XFORS,OPT.DFORS);
+ fprintf(fid,'C.4 DEPTH CORRESPONDING TO LANDWARD END OF SURF ZONE: DFS\r\n');
+ fprintf(fid,' %.3f\r\n',OPT.DFS);
+ fprintf(fid,'C.5 EFFECTIVE GRAIN SIZE DIAMETER IN MILLIMETERS: D50\r\n');
+ fprintf(fid,' %.3f\r\n',OPT.D50);
+ fprintf(fid,'C.6 MAXIMUM PROFILE SLOPE PRIOR TO AVALANCHING IN DEGREES: BMAX\r\n');
+ fprintf(fid,' %.3f\r\n',OPT.BMAX);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % D - BEACH FILL
+ fprintf(fid,'D------------------------------ BEACH FILL --------------------------------D\r\n');
+ fprintf(fid,'D.1 IS A BEACH FILL PRESENT? (NO=0, YES=1): IBCHFILL\r\n');
+ fprintf(fid,' %d\r\n',OPT.IBCHFILL);
+ fprintf(fid,'D.2 COMMENT: IF NO BEACH FILL, CONTINUE TO E.\r\n');
+ fprintf(fid,'D.3 POSITION OF START AND END OF BEACH FILL RELATIVE\r\n');
+ fprintf(fid,' TO INITIAL PROFILE: XBFS, XBFE\r\n');
+ fprintf(fid,' %.2f %.2f\r\n',OPT.XBFS,OPT.XBFE);
+ fprintf(fid,'D.4 NUMBER OF REPRESENTATIVE POINTS BETWEEN START\r\n');
+ fprintf(fid,'AND END OF BEACH FILL: NFILL\r\n');
+ fprintf(fid,' %d\r\n',OPT.NFILL);
+ fprintf(fid,'D.5 LOCATION AND ELEVATION OF REPRESENTATIVE POINTS RELATIVE TO THE\r\n');
+ fprintf(fid,'INITIAL PROFILE: (XF(I), EFILL(I), I=1,NFILL)\r\n');
+ for ii = 1:length(OPT.NFILL)
+ if isempty(OPT.XF)
+ fprintf(fid,'\r\n');
+ else
+ fprintf(fid,' %.2f %.2f',OPT.XF(ii),OPT.EFILL(ii));
+ end
+ end
+ fprintf(fid,'\r\n');
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % E - SEAWALL/REVETMENT
+ fprintf(fid,'E--------------------------- SEAWALL/REVETMENT ----------------------------E\r\n');
+ fprintf(fid,'E.1 IS A SEAWALL PRESENT? (NO=0, YES=1): ISWALL\r\n');
+ fprintf(fid,' %d\r\n',OPT.ISWALL);
+ fprintf(fid,'E.2 COMMENT: IF NO SEAWALL, CONTINUE TO F.\r\n');
+ fprintf(fid,'E.3 LOCATION OF SEAWALL RELATIVE TO INITIAL PROFILE: XSWALL\r\n');
+ fprintf(fid,' %.2F\r\n',OPT.XSWALL);
+ fprintf(fid,'E.4 IS SEAWALL ALLOWED TO FAIL? (NO=0, YES =1): ISWFAIL\r\n');
+ fprintf(fid,' %d\r\n',OPT.ISWFAIL);
+ fprintf(fid,'E.5 COMMENT: IF NO SEAWALL FAILURE, CONTINUE TO F.\r\n');
+ fprintf(fid,'E.6 PROFILE ELEVATION AT SEAWALL WHICH CAUSES FAILURE, TOTAL WATER ELEVATION\r\n');
+ fprintf(fid,' AT SEAWALL WHICH CAUSES FAILURE, AND WAVE HEIGHT AT SEAWALL WHICH CAUSES\r\n');
+ fprintf(fid,' FAILURE: PEFAIL, WEFAIL,HFAIL\r\n');
+ fprintf(fid,' %.2f %.2f %.2f\r\n',OPT.PEFAIL,OPT.WEFAIL,OPT.HFAIL);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % F - HARD BOTTOM
+ fprintf(fid,'F------------------------------ HARD BOTTOM -----------------------------F\r\n');
+ fprintf(fid,'F.1 IS A HARD BOTTOM PRESENT? (NO=0, YES=1): IHBOT\r\n');
+ fprintf(fid,' %d\r\n',OPT.IHBOT);
+ fprintf(fid,'F.2 COEFFICIENT FOR TRANSPORT GROWTH DOWNDRIFT EXPOSED HARD BOTTOM: SCF\r\n');
+ fprintf(fid,' %.2f\r\n',OPT.SCF);
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ fprintf(fid,'----------------------------------- END ------------------------------------\r\n');
+
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/writetofileelv.m b/writetofileelv.m
new file mode 100644
index 0000000..e03712d
--- /dev/null
+++ b/writetofileelv.m
@@ -0,0 +1,18 @@
+function writetofileelv(elvfile, outpath, pfstring)
+ fid = fopen([outpath '\' pfstring '.ELV'],'w');
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C> * SBEACH water elevation input file: %s.ELV *\r\n',pfstring);
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C> TOTAL WATER ELEVATION (m or ft)\r\n');
+ fprintf(fid,'E>----------------------------------------------------------------\r\n');
+ for ii = 1:length(elvfile)
+ if ii == length(elvfile)
+ fprintf(fid,'%.3f\r\n', elvfile(ii));
+ else
+ fprintf(fid,'%.3f\r\n', elvfile(ii));
+ end
+ end
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/writetofilehdb.m b/writetofilehdb.m
new file mode 100644
index 0000000..304e1c9
--- /dev/null
+++ b/writetofilehdb.m
@@ -0,0 +1,20 @@
+function writetofilehdb(hdbfile, outpath, namestr)
+ fid = fopen([outpath '\' namestr '.HDB'],'w');
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C> * SBEACH hard bottom input file: %s.HDB *\r\n',namestr);
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C> NUMBER OF POINTS\r\n');
+ fprintf(fid,'C> POSITION (m or ft), ELEVATION (m or ft)\r\n');
+ fprintf(fid,'E>----------------------------------------------------------------\r\n');
+ fprintf(fid,'%d\r\n',length(hdbfile));
+ for ii = 1:length(hdbfile)
+ if ii == length(hdbfile)
+ fprintf(fid,'%.3f %.3f\r\n', hdbfile(ii,1),hdbfile(ii,2));
+ else
+ fprintf(fid,'%.3f %.3f\r\n', hdbfile(ii,1),hdbfile(ii,2));
+ end
+ end
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/writetofilepri.m b/writetofilepri.m
new file mode 100644
index 0000000..69dff44
--- /dev/null
+++ b/writetofilepri.m
@@ -0,0 +1,20 @@
+function writetofilepri(prifile, outpath, namestr)
+ fid = fopen([outpath '\' namestr '.PRI'],'w');
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C> * SBEACH initial profile input file: %s.PRI *\r\n',namestr);
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C> NUMBER OF POINTS\r\n');
+ fprintf(fid,'C> POSITION (m or ft), ELEVATION (m or ft)\r\n');
+ fprintf(fid,'E>----------------------------------------------------------------\r\n');
+ fprintf(fid,'%d\r\n',length(prifile));
+ for ii = 1:length(prifile)
+ if ii == length(prifile)
+ fprintf(fid,'%.3f %.3f\r\n', prifile(ii,1),prifile(ii,2));
+ else
+ fprintf(fid,'%.3f %.3f\r\n', prifile(ii,1),prifile(ii,2));
+ end
+ end
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/writetofileprm.m b/writetofileprm.m
new file mode 100644
index 0000000..a5661a2
--- /dev/null
+++ b/writetofileprm.m
@@ -0,0 +1,20 @@
+function writetofileprm(prmfile, outpath, pfstring, rundescription)
+ fid = fopen([outpath pfstring '\' pfstring '_' rundescription '.PRM'],'w');
+ fprintf(fid,'C> *******************************************************\n');
+ fprintf(fid,'C> * SBEACH initial profile input file: %s.PRM *\n',[pfstring '_' rundescription]);
+ fprintf(fid,'C> *******************************************************\n');
+ fprintf(fid,'C>\n');
+ fprintf(fid,'C>\n');
+ fprintf(fid,'C> NUMBER OF POINTS\n');
+ fprintf(fid,'C> POSITION (m or ft), ELEVATION (m or ft)\n');
+ fprintf(fid,'E>----------------------------------------------------------------\n');
+ fprintf(fid,'%d\n',length(prmfile));
+ for ii = 1:length(prmfile)
+ if ii == length(prmfile)
+ fprintf(fid,'%4.3f %3.3f', prmfile(ii,1),prmfile(ii,2));
+ else
+ fprintf(fid,'%4.3f %3.3f\n', prmfile(ii,1),prmfile(ii,2));
+ end
+ end
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/writetofilewav.m b/writetofilewav.m
new file mode 100644
index 0000000..f7fe3a1
--- /dev/null
+++ b/writetofilewav.m
@@ -0,0 +1,18 @@
+function writetofilewav(wavfile, outpath, pfstring)
+ fid = fopen([outpath '\' pfstring '.WAV'],'w');
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C> * SBEACH wave input file: %s.WAV *\r\n',pfstring);
+ fprintf(fid,'C> *******************************************************\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C>\r\n');
+ fprintf(fid,'C> HEIGHT (m or ft), PERIOD (sec)\r\n');
+ fprintf(fid,'E>----------------------------------------------------------------\r\n');
+ for ii = 1:length(wavfile)
+ if ii == length(wavfile)
+ fprintf(fid,'%.3f %.3f\r\n', wavfile(ii,1),wavfile(ii,2));
+ else
+ fprintf(fid,'%.3f %.3f\r\n', wavfile(ii,1),wavfile(ii,2));
+ end
+ end
+ fclose(fid);
+end
\ No newline at end of file