L = 50e3;
L0 = 5e3;
DeltaL0 = 1000;

dx = 10;
x = 0:dx:L;

h = 15;
lambda_d = 137;
A_d = 0.7;
H_d = 0.5;
D = 1;
y = [1,101];

% functions params
b = @(x) 1./(1+exp(-4*(x-L0)/DeltaL0)) + 1./(1+exp(-4*(L-x-L0)/DeltaL0))-1;
if A_d == 0 
    Lambda_d = @(x) sin(2*pi*x/lambda_d);
elseif abs(A_d) < 1
    Lambda_d = @(x) 2/(pi*A_d).* atan(sin(pi*A_d/2) .*sin(2*pi.*x./lambda_d)./(1+sin(pi*A_d/2).*cos(2*pi.*x./lambda_d)));
else
    error('Value of A_d is not allowed')
end
f = figure
f.Position = [100,100,500,500]
subplot(211)
fplot(b,[0,L])
subplot(212)
fplot(Lambda_d,[0,lambda_d*5])
mx_lambda = max(Lambda_d(x))

z_b = -h + b(x).*H_d/2.* (Lambda_d(x) ./ mx_lambda - D);

f = figure
f.Position = [100,100,500,500]
plot(x,z_b,'-o')

Z = meshgrid(-[z_b,h],ones(1,3));
OK = wlgrid('write','FileName',['grd_dx',int2str(dx),'.grd'],'X',x,'Y',y,'AutoEnclosure');

if H_d > 1
    Hdstring = ['_Hd', int2str(round(H_d*100))];
else
    Hdstring = ['_Hd0', int2str(round(H_d*100))];
end
lambdastring = ['_l',int2str(round(lambda_d))];
astring = ['_a0', int2str(round(A_d*100))];

OK = wldep('write',['dep_dx',int2str(dx), Hdstring, lambdastring, astring,'_D',int2str(D),'.dep'],'',Z');