% 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