% profile; x positive seaward
[x, z] = deal([0 50 65 115 355], [7.5 7.5 2.5 0 -20]);

% change profile direction: x positive landward
[x, ix] = sort(-x);
z = z(ix);

% water level
WL = 0;
% create x-grid
[xgr, zgr] = xb_grid_xgrid(x, z,...
    'wl', WL);
% add 3 extra cells at the offshore boundary
xgr = [xgr(1) - diff(xgr(1:2)) * (3:-1:1) xgr];
zgr = [repmat(zgr(1), 1, 3) zgr];


dymin = 5; % desired dy at waterline
r = 115; % radius at waterline
ymax = 2*pi*r; % circumference at waterline
ny = ceil(ymax/dymin); % number of alongshore cells
y = linspace(0, ymax, ny); % alongshore dimension at waterline
y(end+1) = y(end) + diff(y(1:2)); % at one extra alongshore cell for overlap

% x coordinate at waterline
x_contour = findCrossings(xgr, zgr, minmax(xgr(:)')', ones(2,1)*WL);
% amplitude
a = r + x_contour - xgr(:)';

% mean orientation
phi_mean = pi;
% orientation range
phi_range = diff(y([1 end]))/(2*pi*r) * 2 * pi * [-.5 .5] + phi_mean;
% difference in orientation between 2 adjacent transects
d_phi = atan(diff(y(1:2))/r);
% number of alongshore cells
ny = max( [ceil(abs(diff(phi_range)/d_phi)) 3] );
% evenly distribute the orientation angles
phi = linspace(phi_range(1), phi_range(2), ny);
% sort angles descending
phi = sort(phi', 'descend'); %#ok<UDIM>

% calculate x, y and z
x = cos(phi)*a;
y = sin(phi)*sign(r)*a;
z = repmat(zgr(:)', size(x,1), 1);

%% write grid in delft3d format
calcdir = fileparts(which(mfilename));
ok1 = wlgrid('write', fullfile(calcdir, 'xy.grd'), x, y);

% add -999 row and column for quickin
depw = ones(size(z) + 1) * -999;
depw(1:end-1, 1:end-1) = z;

ok2 = wldep('write', fullfile(calcdir, 'bed.dep'), depw');

%%
plot(x(:,1), y(:,1), x(1,:), y(1,:), x(end,:), y(end,:))
hold on
for i = 1:10:size(x,2); text(x(1,i),y(1,i), num2str(i)); end
for i = 1:10:size(x,1); text(x(i,1),y(i,1), num2str(i)); end

%%
ortho_grid = [xgr 0 fliplr(-xgr)];
[xo, yo] = meshgrid(ortho_grid-min(xgr), ortho_grid-min(xgr));
ao = sqrt((xo+min(xgr)).^2 + (yo+min(xgr)).^2);
zo = interp1(abs(xgr), zgr, ao, 'linear', 'extrap');

%% write grid in delft3d format
calcdir = fileparts(which(mfilename));
ok1 = wlgrid('write', fullfile(calcdir, 'xy_ortho.grd'), xo, yo);

% add -999 row and column for quickin
depwo = ones(size(zo) + 1) * -999;
depwo(1:end-1, 1:end-1) = zo;

ok2 = wldep('write', fullfile(calcdir, 'bed_ortho.dep'), depwo');
