% function makeprofile_mic_variable_dx
clear all
close all

addpath(genpath('E:\My software'));

pp0=cd;

% This function will work for beach profiles with depths defined by
% negative values and x values increasing shorewards
% it generates a beach profile grid with variable grid spacing by

% a) making sure that offshore n values don't exceed 0.85
% b) the ratio sqrt(depth)/dx remains constant in the whole domain


% Michalis Vousdoukas, 2008


% initial grid spacing (for interpolated profiles input your grid spacing) 
dx=2;

dx0=20;     % deep water dx (maximum)
dxs=0.5;      %minimum dx

% extension on the shore for len meters
len=10;


% maximum deep water parameters
per=8;
ho=4;

% slope for deep water extension of the profile
sl=1/30;

% maximum possible deep water depth
d0=-20;

% % input data (if your data are not interpolated-make sure depth values
% % are negative!)


% profname='profile.dat';
% 
% xz=load(profname)
% xdata1=xz(:,1);
% zdata1=xz(:,2);
% 
% xdata=[xdata1(1):dx:xdata1(end)];
% zdata=interp1(xdata1,zdata1,xgrid);

cd ..

load profs_2


ic=make_icol_scheme;

cd(pp0)

dx=2;

load(['start.dep'])

nst = start(1,:);

    
% my input for my example
% profname='profile.dat';
% 
% path='D:\Algarve\Field work\Data\Crop\';
% 
% nam='data_mix.mat';
% 
% load([path nam]);

% xdata=([1:length(nst)]-1)*dx;
% zdata=nst;

pn=3
xdata=prof(pn).x1;
zdata=flipud(prof(pn).z1);


% exclude nans
ind1=find(isnan(zdata)==0);

xdata=xdata(ind1);
zdata=zdata(ind1);

figure(1)

x1=abs(xdata-min(xdata))
plot(x1,zdata,'o-')
title('Initial profile')

title('Original profile')

xlabel('cross-shore distance')
ylabel('elevation')
grid on 
box on

% vertical increment
dy=-sl*dx;

% depths of the section added offshore
d1=[d0:-dy:min(zdata)];

% size
np=length(d1);

% cross-shore distance of the section added offshore
x2=flipud(min(x1)-dx*[1:np]');


% merging
d=[d1(:) ; zdata(:)];

d2=abs(d);

% estimate n parameter
for i=1:length(d)
    
    l(i)=dispeq(50,d2(i),per);

    kapa=2*pi/l(i);

    n(i)=0.5*(1+2*(kapa*d2(i)/sinh(2*kapa*d2(i))));
    
end

% criterion according to depth
ind=find(n<=0.85);

% find maximum depth so that offshore n won't exceed 0.85
if length(ind)>0
    
    mdep=max(d(ind))
    
else
    
    mdep=d0;
    
end

if mdep<-d0
    
    mdep=-d0;
    
end

% make flat for depth>mdep

d(find(d<=-mdep))=-mdep;

% extension on the shore for len meters
x3=[max(x1)+dx:dx:max(x1)+len]';

% merging with the existing profile
d=[d ; ones(length(x3),1).*max(d)];

nx=[x2(:) ;x1(:) ; x3(:)];

% plot
hold on

plot(nx,d,'r')

legend('raw','extended',2)

title('Constant dx grid extended')

% grid ratio to satisfy Courant
gratio_d=sqrt(abs(mdep))/dx0;

xx=nx-min(nx);

% first estimation of grid spacing
dxx=sqrt(abs(d))./gratio_d;

figure(2)
subplot(2,1,1)
plot(xx,d)

title('Final profile')

xlabel('cross-shore distance')
ylabel('elevation')
grid on 
box on

hold on

subplot(2,1,2)
plot(xx,dxx)

title('Grid spacing')

xlabel('cross-shore distance')
ylabel('grid spacing-dx')
grid on 
box on

hold on

% round to estimate grid spacing
dxx1=round(dxx);

% dry profile
dxx1(find(d>=0))=dxs;

plot(xx,dxx1,'r')


ind2=find(diff(dxx1)~=0);

plot(xx(ind2),dxx(ind2),'o');

legend('estimated','final','separators',2)

ind3=[0 ; ind2(:) ; length(dxx1)];

dxx2=[dxx1(ind2); min(dxx1)];

fx=[];

for i=1:length(ind3)-1
    
    if i==1
        
        temp=[xx(ind3(i)+1):dxx2(i):xx(ind3(i+1))];
        
    else
        
        temp2=max(temp);
        
        temp=[temp2+dxx2(i):dxx2(i):xx(ind3(i+1))];
        
    end
    
    fx=[fx; temp(:)];
      
end

indp=find(diff(fx)>0);

fx=fx(indp);

nd=interp1(xx,d,fx)

subplot(2,1,1)

hold on
plot(fx,nd,'g.')

legend('initial','final',2)

figure(3)

subplot(2,1,1)
plot(fx)
title('Cross-shore distance')

xlabel('grid points')
ylabel('x')
grid on 
box on

subplot(2,1,2)
plot(diff(fx))

title('Grid spacing final validation')

xlabel('grid points')
ylabel('diff(x)')
grid on 
box on

% final grid adding to the cross shore


[xgrid,ygrid]=meshgrid(fx,[0 10 20]);
zgrid=nd'%[nd nd nd]';

% ygrid=[0 10 20];

[ny,nx]=size(xgrid);

fi=fopen('test_griddata.txt','wt');
fprintf(fi,'nx      = %3i \n',nx-1);
fprintf(fi,'ny      = %3i \n',ny-1);
% fprintf(fi,'dx      = %6.1f \n',dx);
% fprintf(fi,'dy      = %6.1f \n',dx);
fprintf(fi,'xori    = %10.2f \n',0);
fprintf(fi,'yori    = %10.2f \n',0);
fprintf(fi,'alfa    = %10.2f \n',0);
fprintf(fi,'depfile = newdep.dep \n')
fprintf(fi,'posdwn  = -1 \n')
fprintf(fi,'vardx  = 1 \n')
fprintf(fi,'xfile    = test_x.txt  \n');
fprintf(fi,'yfile    = test_y.txt  \n');

fclose(fi);

fi=fopen('newdep.dep','wt');
for i=1:3
    fprintf(fi,'%7.3f ',zgrid);
    fprintf(fi,'\n');
end
fclose(fi);

fi=fopen('test_x.txt','wt');
for i=1:3
    fprintf(fi,'%7.3f ',xgrid(i,:));
    fprintf(fi,'\n');
end
fclose(fi);

fi=fopen('test_y.txt','wt');
for i=1:3
    fprintf(fi,'%7.3f ',ygrid(i,:));
    fprintf(fi,'\n');
end
fclose(fi);


% function l=dispeq(lo,d,per)
% 
% % solves disperion equation
% omega=2*pi/per;
% g=9.81;
% check=1.;
% 		
% l1=lo    ;
% 
% %     dispersion relation gia to ypologismo tou wave length aL
% while check>0.000001
%     kappa=2*pi/l1;
%     kd=kappa*d;
%     l2=g*2*pi*tanh(kd)/omega^2;	
%     check=abs(l2-l1);
%     l1=(l2+l1)/2;
% end
% l=l2;
%    
% end
% 
% end