OPT = struct(...
    'id', '',...
    'idformat', '%s',...
    'irun', 42,...
    'irunformat', 'run%02i',...
    'model', '',...'xbeach',...
    'modelformat', '%s',...
    'maindir', fullfile(fileparts(mfilename('fullpath')), '..', 'simulations'),...
    'WL', 5,...
    'Hsig', 9,...
    'Tp', 12,...
    'D50', 225e-6,...
    'duration', 5*3600,...
    'initialdurationmarge', 400,...
    'radii', roundoff(1000 ./ (pi*(30:-6:6)/180), -2),...
    'dtheta', 5,...
    'thetarange', 180,...
    'circleportion', .5,...
    'dxmin', 5,...
    'dymin', 20,...
    'morfac', 5,...
    'profilefun', {{@SimpleProfile 'zmin', -20}});%{{@deal [0 50 65 115 355] [7.5 7.5 2.5 0 -20]}});

%% define cross-shore profile
[x z] = feval(OPT.profilefun{:});
[x ix] = sort(-x);
z = z(ix);

% derive required depth at seaward boundary
zsea = min(floor([min(OPT.WL(:)) - max(OPT.Tp(:))/.55 -20]));
xb_setpref('grid_finalise', {'zmin', zsea})

%
tmarge = roundoff(OPT.initialdurationmarge * OPT.morfac, -2, @ceil);
if OPT.morfac > 0
    duration = ceil(OPT.duration/OPT.morfac)*OPT.morfac;
else
    duration = OPT.duration;
end
tstop = duration + tmarge;

%
xb = xb_generate_model(...
    'bathy', {'x', x, 'z', z, 'xgrid', {'dxmin', OPT.dxmin}, 'finalise', {'actions', {'seaward_extend'}}, 'superfast', true},...
    'tide', {'front' OPT.WL},...
    'waves', {...
    'type', 'jonswap',...
    'Hm0', OPT.Hsig,...
    'Tp', OPT.Tp,...
    ...'s', 1000, ...
    'duration', min([tstop 100*OPT.morfac*OPT.Tp*4]),...
    'timestep', 2,...
    'fnyq', 1/OPT.Tp *4 ...
    },...
    ...'wavegrid', {'nbins' 1},...
    'settings', {...
    'morfac', OPT.morfac,...
    'posdwn', -1,...
    'tintg', 200,...
    'tintm', 1000,...
    'tstop', tstop,...
    'D50', OPT.D50,...
    'D90', OPT.D50*1.5,...
    'random' 0,...
...    'avalanching', 0,...
...    'facDc', 0,...
...    'depthscale', 1,...
    'nc', 20,...
    'Tsmin', 3,...
    'epsi', -1,...
    'outputformat', 'netcdf',...
    'globalvars', {'H' 'ee' 'zs' 'zb' 'sedero' 'ue' 've' 'u' 'v' 'Susg' 'Svsg' 'Subg' 'Svbg' 'ccg' 'dzav' 'ceqsg' 'thetamean' 'cgx' 'cgy' 'Fy' 'Fx' 'umean' 'vmean' 'Qb' 'D'},...
    'meanvars', {'H' 'ee' 'zs' 'ue' 've' 'u' 'v' 'Susg' 'Svsg' 'Subg' 'Svbg' 'ccg' 'ceqsg' 'thetamean' 'cgx' 'cgy' 'Fy' 'Fx' 'umean' 'vmean' 'Qb' 'D'},...
    'gridform', 'delft3d',...
    'rotate', 0,...
...    'left', 'wall',...
...    'right', 'wall',...
    });

% wl = xb_bc_extractwl(xb);
% xb = xb_set_start_time(xb, 'waterlevel', wl);
% tstart = xb_get(xb, 'tstart');
xb = xs_set(xb, 'tstart', tmarge, 'morstart', tmarge);

x = xs_get(xb, 'xfile.xfile');
z = xs_get(xb, 'depfile.depfile');
xb = xs_del(xb, 'xfile', 'yfile', 'swtable', 'nx', 'ny', 'vardx');

xb = xs_set(xb,...
    'xyfile', 'xy.grd',...
    'depfile', 'bed.dep');


mainangle = xs_get(xb,'bcfile.mainang');
remainer = mod(.5*OPT.thetarange-.5*OPT.dtheta, OPT.dtheta);
thetarange = mainangle + [-1 1]*.5*(OPT.thetarange-2*remainer);
xb = xs_set(xb,...
    'thetamin', thetarange(1),...
    'thetamax', thetarange(2),...
    'dtheta', OPT.dtheta);

basepath = fullfile(OPT.maindir, sprintf(OPT.idformat, OPT.id), sprintf(OPT.modelformat, OPT.model), sprintf(OPT.irunformat, OPT.irun));

for iradii = 1:length(OPT.radii)
    calcdir = fullfile(basepath, sprintf('radius%05i', OPT.radii(iradii)));
    mkpath(calcdir)
    
    ymax = 2 * pi * OPT.radii(iradii) * OPT.circleportion;
    [xgr ygr dep] = decc_curved_grid(...
        'x', x',...
        'y', [0 ymax],...
        'z', z',...
        'contour', OPT.WL,...
        'dymin', OPT.dymin,...
        'radius', OPT.radii(iradii));
    dep = dep';
    xori = min(xgr(:));
    yori = min(ygr(:));
    xb = xs_set(xb, 'xori', xori, 'yori', yori);
    xgri = xgr + xori;
    ygri = ygr + yori;
    
    xb_write_input(fullfile(calcdir, 'params.txt'), xb);%,...
        %'params_header', {'XBeach parameter settings input file' '' ['date:     ' datestr(now)] ['function: ' mfilename] '$Id$'});
    
    ok1 = wlgrid('write', fullfile(calcdir, 'xy.grd'), xgri, ygri);
    
    % add -999 row and column for quickin
    depw = ones(size(dep) + [1 1]) * -999;
    depw(1:end-1, 1:end-1) = dep;
%     depw2 = [dep -999*ones(size(dep,1),1); ...
%         -999*ones(1,size(dep,2)+1)];
    
    ok2 = wldep('write', fullfile(calcdir, 'bed.dep'), depw');
end