%%plot vector
%Gundula Winter
%2013-02-12

%% clean up
clc
clear
close all

%% input
dx=0.01;
x=0:dx:14;
y=0;
t=0:0.01:6;
[xgrid,ygrid]=meshgrid(x,y);

%% profile
z0=zeros(1,length(x));
z0(x>=1.807) = -0.435;
z0(x<1.807) = -(-0.19*(x(x<=1.807)-1.297).^3 - 0.1024*(x(x<=1.807)-1.297).^2 + 0.2728*(x(x<=1.807)-1.297) + 0.3475);
z0(x<=1.297) = -x(x<=1.297)*tand(15);
z0(x<=0.03) = -0.009;
z0=repmat(z0,length(y),1);
% plot(x,z0)
% hold on
% line([1.297 1.297],[-.5 0]);
% line([1.807 1.807],[-.5 0]);
% xlim([0 2])

%% velocity
a=1*ones(1,length(t));
v=zeros(1,length(t));
pos=zeros(1,length(t));
v(t<=1.2)=(a(t<=1.2)).*(t(t<=1.2));
pos(t<=1.2)=0.5*(a(t<=1.2)).*(t(t<=1.2)).^2;
a(t>1.2)=-t(t>1.2);
v(t>1.2)=v(t==1.2)+a(t>1.2).*(t(t>1.2)-1.2);
pos(t>1.2)=pos(t==1.2)+v(t==1.2).*(t(t>1.2)-1.2)+0.5*(a(t>1.2)).*(t(t>1.2)-1.2).^2;
a(t>1.31)=-0.25;
v(t>1.31)=v(t==1.31)+a(t>1.31).*(t(t>1.31)-1.31);
pos(t>1.31)=pos(t==1.31)+v(t==1.31).*(t(t>1.31)-1.31)+0.5*(a(t>1.31)).*(t(t>1.31)-1.31).^2;
a(v<0)=[];
pos(v<0)=[];
t(v<0)=[];
v(v<0)=[];
% plot(t,a,'r')
% hold on
% plot(t,v)
% grid on
% figure; plot(t,pos); grid on


%% landslide files and animation
z_abs=zeros(length(t),length(y),length(x));
block=zeros(length(y),length(x));
for i=1:length(t)
    % slide body    
    x0=0.4665+pos(i);
    xb=0:dx:.5;
    %     block=0.982*0.026*(1-(xb./0.25).^4);
    block=0.5*(1+tanh(tan(pi*((xb./0.25)-0.5))));
    block(xb>0.25)=1-block(xb>.25);
    block=.026*block;
    if x0<1.807
        for j=1:length(xb)
           xbrot(j)=cosd(15)*xb(j)-sind(15)*block(j);
           blockrot(j)=sind(15)*xb(j)+cosd(15)*block(j);
        end
        blockrot=blockrot-xbrot*tand(15);
        block=interp1(xbrot,blockrot,xb,'linear',0);
    end
    block=interp1(xb,block,x-x0,'linear',0);
    block=repmat(block,length(y),1);
    z_abs(i,:,:)=z0+block;
end

fb=figure;
axes('position',[.15 .15 .3 .7])
plot(xb,block)
xlim([0 0.5])
axes('position',[.6 .15 .3 .7])
plot(xb(1:end-1),diff(block))
xlim([0 0.5])
set(fb,'paperunits','centimeters', 'papersize',[15 4],...
    'paperposition',[0, 0, 15, 4],'renderer','zbuffer')
print(fb,'-dpng','block_gladijs.png')

zstart=squeeze(z_abs(1,:,:));
zstart=zstart';
% z0=-0.1*ones(1,length(z0));
z=diff(z_abs);

% animation to check

% z_temp=zstart;%*0;
% fig=figure;
% for k=2:length(t)-1
%     z_temp=z_temp+(squeeze(z(k,:,:)))';
%     plot(x,z_temp)
% %     plot(x,(squeeze(z(k,:,:)))')
%     set(gca,'position',[.1 .1 .85 .3])
%     pause(0.005)
% end

% 
% 
% aviobj = avifile('Sue2007.avi','compression','Indeo5','keyframe',2.1429*50);
% z_temp=z0;
% fig=figure;
% for k=1:length(t)-1
%     z_temp=z_temp+(squeeze(z(k,:,:)))';
%     plot(x,z_temp)
%     set(gca,'position',[.1 .1 .85 .3])
% 
%     F = getframe(fig);
%     aviobj = addframe(aviobj,F);
% end
% close(fig);
% aviobj = close(aviobj);




%% write files for XBeach

mkdir('Sue_xbeach_gladijs')



xfile = fopen('Sue_xbeach_gladijs\x.grd', 'w');
for i=1:size(xgrid,1);
    fprintf(xfile,'%13E ',xgrid(i,:));
    fprintf(xfile,'\n');
end
fclose(xfile);

% y
yfile = fopen('Sue_xbeach_gladijs\y.grd', 'w');
for i=1:size(ygrid,1);
    fprintf(yfile,'%13E ',ygrid(i,:));
    fprintf(yfile,'\n');
end
fclose(yfile);

% z
zfile = fopen(['Sue_xbeach_gladijs\bed.dep'], 'w');
for i=1:size(z0,1);
    fprintf(zfile,'%13E ',zstart(i,:));
    fprintf(zfile,'\n');
end
fclose(zfile);

% landslide files
lsfile = fopen('Sue_xbeach_gladijs\landslide.txt', 'w');
for j=2:length(t)  % Do not start at 0 seconds because we do not want to update the bed before the frist time step
    % that would mean we also update the water level before the first time
    % step, but we want to statr with an even water level
    zfile = fopen(['Sue_xbeach_gladijs\',num2str(j-1,'%3d'),'.dep'], 'w');
    fprintf(lsfile,'%s %s \t',num2str(t(j),'%7.3f'),[num2str(j-1,'%3d'),'.dep']);
    fprintf(lsfile,'\n');
    for i=1:size(z,2);
        fprintf(zfile,'%13E ',z(j-1,i,:));
        fprintf(zfile,'\n');
    end
    fclose(zfile);
end
fclose(lsfile);

%% write files for Delft3D

% % x=repmat(x,3,1);
% y=repmat(y,3,1);
% [xgrid,ygrid]=meshgrid(x,y);
% z0=repmat(z0,3,1);
% z_d3d=zeros(size(z,1),size(xgrid,1),size(x,2));
% for i=1:size(z,1)
%     z_d3d(i,:,:)=repmat((squeeze(z(i,:,:)))',size(xgrid,1),1);
% end
% 
% d3d_path=('Sue_Delft3D');
% mkdir(d3d_path)
% 
% grdfile = ([d3d_path,'\Sue2007.grd']);
% OK=wlgrid('write','FileName',grdfile,'X',xgrid,'Y',ygrid);
% 
% d3d_grid=wlgrid('read',grdfile);
% 
% Frmt=repmat('%15.8f  ',[1 size(d3d_grid.X,1)]);
% k=8*12;
% Frmt((k-1):k:length(Frmt))='\';
% Frmt(k:k:length(Frmt))='n';
% Frmt(end-1:end)='\n';
% depfile = ([d3d_path,'\Sue2007.dep']);
% wldep('write',depfile,'',z0)
% % fid=fopen(depfile,'w');               
% % z0_d3d=[z0 -999*ones(size(z0,1),1); ...
% %                   -999*ones(1,size(z0,2)+1)]; 
% % fprintf(fid,Frmt,z0_d3d);        
% % fclose(fid);
% 
% % landslide files
% lsfile = fopen([d3d_path,'\landslide.txt'], 'w');
% for j=2:length(t)
%     zfile = ([d3d_path,'\',num2str(j-1,'%3d'),'.dep']);
%     fprintf(lsfile,'%s %s \t',num2str(t(j-1),'%7.3f'),[num2str(j-1,'%3d'),'.dep']);
%     fprintf(lsfile,'\n');
%     wldep('write',zfile,'',(squeeze(z_d3d(j-1,:,:))))
% end
% fclose(lsfile);

%% netdf
% nccreate('test.nc','x','Dimensions',{'nx',length(x)});
% ncwrite('test.nc','x',x);
% nccreate('test.nc','y','Dimensions',{'ny',length(y)});
% ncwrite('test.nc','y',y);
% nccreate('test.nc','t','Dimensions',{'nt',length(t)});
% ncwrite('test.nc','t',t);
% nccreate('test.nc','z','Dimensions',{'nt' 'ny' 'nx'},'ChunkSize', [length(t) length(y) length(x)]);
% ncwrite('test.nc','z',z(length(t),length(y),length(x)));


