Commit 49b4c29c authored by Jérémy Baudry's avatar Jérémy Baudry

clean matlab scripts

parent 39385008
File deleted
clear all;
close all;
fig=figure;
% fig.Units='inches';
% fig.Position=[6 6 9 3.2];
% fig.OuterPosition=[6 6 9 3.2];
%
% fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 0 9 3.2];
% fig.PaperPositionMode = 'manual';
% fig.PaperOrientation='landscape';
% fig.PaperSize=[9 3.2];
%
for n=3
titre1=['10^{-6}';'10^{-4}';'10^{-3}';'10^{-2}'];
titre2=cellstr(titre1);
dos1=['1e6 ';'1e4 ';'test ';'1e2 '];
dos2=cellstr(dos1);
dossier=[num2str(n),'mu1';num2str(n),'mu2';num2str(n),'mu3';dos2(n);num2str(n),'mu4'];
dossier=cellstr(dossier);
diffus=['1e-2_uniforme'];
mix=[1e-6 1e-4 1e-3 1e-2];
mu=[0.5 1 1.5 2 2.5]/86400;
col=['mo';'ro';'bo';'ko'];
col=cellstr(col);
for i=1:5
dos=char(dossier(i));
number(:,n)=sqrt(mu/mix(n));
lagbiomass=load(['~/projets/correction_memo/output/',dos,'/weight_',diffus,'.dat']);
eulbiomass=load(['~/projets/correction_memo/output/',dos,'/eulerian_',diffus,'.dat']);
biolag_tot=sum(lagbiomass,2)*0.5;
bioeul_tot=sum(eulbiomass,2)*0.5;
std=sqrt((biolag_tot-bioeul_tot).^2./bioeul_tot.^2)*100;
std(1:3)=0;
maxstd(i,n)=max(std);
maxbio(i,n)=max(biolag_tot);
plot(std)
hold on
end
end
clear all;
close all;
dossier=['temp1'];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
x=ncread(path,'x_axis');
t=ncread(path,'time');
om=ncread(path,'omega');
spectre=ncread(path,'Spectrum');
Dave=ncread(path,'Dave');
Dmax=ncread(path,'Dmax');
FSD=ncread(path,'Floe size distribution');
Fsize=ncread(path,'floe size');
thick=ncread(path,'Ice thickness');
conc=ncread(path,'Ice concentration');
Hs=ncread(path,'significant height');
IDT=ncread(path,'Ice Thickness Distribution');
hcat=ncread(path,'thickness categories');
fig=figure(1);
u = fig.Units;
fig.Units = 'inches';
set(fig,'outerposition',[3 3 12 10])
cm=colormap('jet');
cm1(1,:)=[1 1 1];
for i=2:10
cm1(i,:)=[1-0.08*i 1-0.06*i 1];
end
cmap(1:10,:)=cm1;
cmap(11:59,:)=cm(16:end,:);
w = waitforbuttonpress;
%figure
%filename = 'joint_distribution_anim.gif';
for j=1:length(t)
FSD1=reshape(FSD(j,j,:,:),length(Fsize),length(hcat));
IDT1=IDT(j,:);
for i=1:length(hcat)
FSD1(:,i)=FSD1(:,i)*IDT1(i);
nfloe(:,i)=((FSD1(:,i)*(5000)^2)./Fsize.^2) ;
end
nfloe=nfloe/sum(sum(nfloe));
FSD2(j,:)=sum(FSD1,2);
subplot(2,2,1)
hh=surf(hcat,Fsize,nfloe);
shading interp
%
zoom(1.5)
caxis([0 0.06])
colormap(cmap)
axis([0.1 10 0 400 0 0.1])
view(108,60)
xlabel('Ice thickness [m]')
ylabel('Floe size [m]')
zlabel('area fraction')
subplot(2,2,2)
nfloe2=((FSD2(j,:)*5000^2)./Fsize.^2')/sum((FSD2(j,:)*5000^2)./Fsize.^2');
plot((Fsize),nfloe2)
axis([0 400 0 1])
xlabel('Floe size [m]')
ylabel('normalized number of floe')
subplot(2,2,3:4)
plot(x,thick)
hold on
plot(x(j),thick(j),'ro','markerfacecolor','r')
hold off
xlabel('distance [km]')
ylabel('Average ice thickness [m]')
%pause(0.02)
w = waitforbuttonpress;
% drawnow
% frame = getframe(1);
% im = frame2im(frame);
% [A,map] = rgb2ind(im,256);
% if j == 1;
% imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
% else
% imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
% end
end
clear all;
close all;
x=ncread('simulation1.nc','x_axis');
t=ncread('simulation1.nc','time');
om=ncread('simulation1.nc','omega');
spectre=ncread('simulation1.nc','Spectrum');
Dave=ncread('simulation1.nc','Dave');
Dmax=ncread('simulation1.nc','Dmax');
FSD=ncread('simulation1.nc','Floe size distribution');
Fsize=ncread('simulation1.nc','floe size');
thick=ncread('simulation1.nc','Ice thickness');
conc=ncread('simulation1.nc','Ice concentration');
Hs=ncread('simulation1.nc','significant height');
IDT=ncread('simulation1.nc','Ice Thickness Distribution');
hcat=ncread('simulation1.nc','thickness categories');
f=om/(2*pi);
E=reshape(spectre(end,end,:),length(om),1);
Ei=reshape(spectre(30,:,:),length(x),length(om));
E1=reshape(spectre(1,1,:),length(om),1);
for i=1:length(f)
thick2D(:,i)=thick;
end
thick2D(:,1)=0;
thick2D(:,end)=0;
thick2D(end,:)=0;
fig=figure(1);
u = fig.Units;
fig.Units = 'inches';
set(fig,'outerposition',[3 3 12 7])
cm=colormap('jet');
cm1(1,:)=[1 1 1];
for i=2:10
cm1(i,:)=[1-0.08*i 1-0.06*i 1];
end
cmap(1:10,:)=cm1;
cmap(11:59,:)=cm(16:end,:);
filename = 'spectrum_attenuation_anim.gif';
% w = waitforbuttonpress;
for i=1:length(t)
figure(1)
h2=surf(f,x,thick2D*0.05);
view(-56,22)
hold on
sp=reshape(spectre(i,:,:),length(x),length(om));
h=surf(f,x,sp);
view(-56,22)
colormap(jet)
% caxis([-0.1 0.2])
shading interp
set(h2,'Facealpha',0.3,'facecolor',[0.9 0.9 0.9])
ll=light;
set(ll,'position',[0.25 -20 0.4])
lighting phong
axis([min(f) max(f) min(x) max(x) 0 max(E1)])
xlabel('Frequency [s^{-1}]')
zlabel('Energy')
hold off
ylabel('x [km]')
% pause(0.001)
%w = waitforbuttonpress;
drawnow
frame = getframe(1);
im = frame2im(frame);
[A,map] = rgb2ind(im,256);
if i == 1;
imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
else
imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
end
end
%
%
% for i=1:size(Hs,1)
%
% hs2(i)=Hs(i,i);
% end
%
% figure
% subplot(4,4,1:4)
% [ax,p1,p2] = plotyy(x,thick,x,conc,'plot','plot');
%
% ylabel(ax(1),'Ice thickness [m]')
% ylabel(ax(2),'Ice concentration')
% hold on
% line1=line([55 55], ylim);
% line2=line([150 150], ylim);
% line3=line([250 250], ylim);
% line4=line([400 400], ylim);
% set(line1,'color','k','linestyle','--')
% set(line2,'color','k','linestyle','--')
% set(line3,'color','k','linestyle','--')
% set(line4,'color','k','linestyle','--')
%
%
% subplot(4,4,5:8)
%
% plot(x,hs2)
%
% ylabel('Hs [m]')
%
%
% subplot(4,4,9)
% fsd1=reshape(FSD(:,find(x==55),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==55)),'-')
% xlabel('Floe size [m]')
% ylabel('area fraction')
% subplot(4,4,10)
% fsd1=reshape(FSD(:,find(x==150),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==150)),'-')
% xlabel('Floe size [m]')
%
% subplot(4,4,11)
% fsd1=reshape(FSD(:,find(x==250),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==250)),'-')
% xlabel('Floe size [m]')
%
% subplot(4,4,12)
% fsd1=reshape(FSD(:,find(x==400),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==400)),'-')
% xlabel('Floe size [m]')
%
% subplot(4,4,13:16)
% plot(x,Dave)
% xlabel('Distance [km]')
% ylabel('<D> [m]')
% saveas(gcf,'new_model.png')
%
clear all;
close all;
fig=figure;
fig.Units='inches';
fig.Position=[7.0312 6.6875 4.0312 3.3021];
fig.OuterPosition=[7.0312 6.6875 4.0312 4.3021];
fig.PaperUnits = 'inches';
fig.PaperPosition = [0 0 4.0312 3.3021];
fig.PaperPositionMode = 'manual';
fig.PaperOrientation='landscape';
fig.PaperSize=[4.0312 3.3021];
strain=[0:1e-6:3e-4];
prob=exp(-2*3e-5^2./strain.^2);
plot(strain(find(prob<0.37)),prob(find(prob<0.37)),'k--','linewidth',1.5)
hold on
plot(strain(find(prob>0.37)),prob(find(prob>0.37)),'k-','linewidth',1.5)
a=xlabel('$E_s$');
b=ylabel('$P_{\varepsilon}=P(E_w>\varepsilon_c)$');
set(a,'interpreter','latex','fontsize',11)
set(b,'interpreter','latex','fontsize',11)
scrit=3e-5;
probcrit=0.37;
plot(xlim,[probcrit probcrit])
plot([scrit scrit],ylim)
c=text(1.5e-4,0.4,'$P_c$');
d=text(3.5e-5,0.65,'$\varepsilon_c$');
set(c,'interpreter','latex','fontsize',11)
set(d,'interpreter','latex','rotation',90,'fontsize',11)
saveas(fig,'Pe_Es','pdf')
clear all;
close all;
dossier=['fig1_1'];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
Fsize=ncread(path,'floe size');
strain1=load('/home/bauj0001/projets/WIM/output/strain1.dat');
beta1=load('/home/bauj0001/projets/WIM/output/beta1.dat');
Q1=load('/home/bauj0001/projets/WIM/output/Q1.dat');
strain2=load('/home/bauj0001/projets/WIM/output/strain2.dat');
beta2=load('/home/bauj0001/projets/WIM/output/beta2.dat');
Q2=load('/home/bauj0001/projets/WIM/output/Q2.dat');
strain3=load('/home/bauj0001/projets/WIM/output/strain3.dat');
beta3=load('/home/bauj0001/projets/WIM/output/beta3.dat');
Q3=load('/home/bauj0001/projets/WIM/output/Q3.dat');
Lmin=pi*0.5*((5e6*2^3)/(3*10*(1-0.3^2)))^0.25;
Lmin=Lmin*2;
wavelength=Fsize*2;
fig=figure;
fig.Units='inches';
fig.Position=[7.4167 2.0833 2.9583 5.5312];
fig.OuterPosition=[7.4167 2.0833 2.9583 6.5312];
fig.PaperUnits = 'inches';
fig.PaperPosition = [0 0 2.9583 5.5312];
fig.PaperPositionMode = 'manual';
fig.PaperOrientation='landscape';
fig.PaperSize=[2.9583 5.5312];
subplot(3,1,1)
plot(wavelength(find(wavelength>Lmin)),strain1(find(wavelength>Lmin)),'color',[0.1 0.1 0.1],'linewidth',1.5)
hold on
plot(wavelength(find(wavelength>Lmin)),strain2(find(wavelength>Lmin)),'color',[0.5 0.5 0.5],'linewidth',1.5)
plot(wavelength(find(wavelength>Lmin)),strain3(find(wavelength>Lmin)),'color',[0.8 0.8 0.8],'linewidth',1.5)
a=ylabel('$E_s$');
b=xlabel('$\lambda_{\mathrm{ice}}$ [m]');
set(a,'interpreter','latex')
set(b,'interpreter','latex')
c=legend('Hs=0.1 m','Hs=0.5 m ','Hs=1 m');
set(c,'box','off')
hold on
plot([Lmin Lmin],ylim,'color',[0.85 0.325 0.098])
e=plot(wavelength(find(wavelength<Lmin)),strain1(find(wavelength<Lmin)),'color',[0.1 0.1 0.1],'linewidth',1.5);
set(e,'linestyle',':')
f=plot(wavelength(find(wavelength<Lmin)),strain2(find(wavelength<Lmin)),'color',[0.5 0.5 0.5],'linewidth',1.5);
set(f,'linestyle',':')
g=plot(wavelength(find(wavelength<Lmin)),strain3(find(wavelength<Lmin)),'color',[0.8 0.8 0.8],'linewidth',1.5);
set(g,'linestyle',':')
j=text(55,2e-3,'$\lambda_{MIN}$');
set(j,'interpreter','latex','rotation',90)
subplot(3,1,2)
plot(Fsize,Q1,'color',[0.1 0.1 0.1],'linewidth',1.5)
hold on
plot(Fsize,Q2,'color',[0.5 0.5 0.5],'linewidth',1.5)
plot(Fsize,Q3,'color',[0.8 0.8 0.8],'linewidth',1.5)
a=ylabel('$Q$');
b=xlabel('Floe size $l$ [m]');
set(a,'interpreter','latex')
set(b,'interpreter','latex')
subplot(3,1,3)
plot(Fsize,beta1(end,:),'color',[0.1 0.1 0.1],'linewidth',1.5)
hold on
plot(Fsize,beta2(end,:),'color',[0.5 0.5 0.5],'linewidth',1.5)
plot(Fsize,beta3(end,:),'color',[0.8 0.8 0.8],'linewidth',1.5)
a=ylabel('$\beta(l=400,l'')$');
b=xlabel('Floe size $l$ [m]');
set(a,'interpreter','latex')
set(b,'interpreter','latex')
saveas(fig,'test','pdf')
clear all;
close all;
dossier=['wavenumber'];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
x=ncread(path,'x_axis');
t=ncread(path,'time');
om=ncread(path,'omega');
spectre=ncread(path,'Spectrum');
Dave=ncread(path,'Dave');
Dmax=ncread(path,'Dmax');
FSD=ncread(path,'Floe size distribution');
Fsize=ncread(path,'floe size');
thick=ncread(path,'Ice thickness');
conc=ncread(path,'Ice concentration');
Hs=ncread(path,'significant height');
IDT=ncread(path,'Ice Thickness Distribution');
hcat=ncread(path,'thickness categories');
cm=colormap('jet');
cm1(1,:)=[1 1 1];
for i=2:10
cm1(i,:)=[1-0.08*i 1-0.06*i 1];
end
cmap(1:10,:)=cm1;
cmap(11:59,:)=cm(16:end,:);
alphmap=ones(63,1);
alphmap(1)=0;
xx=[55 70 80 100 200];
for j=1:5
FSD1=reshape(FSD(find(x==xx(j)),find(x==xx(j)),:,:),length(Fsize),length(hcat));
IDT1=IDT(find(x==xx(j)),:);
for i=1:length(hcat)
FSD1(:,i)=FSD1(:,i)*IDT1(i);
nfloe(:,i)=((FSD1(:,i)*(5000)^2)./Fsize.^2) ;
end
nfloe=nfloe/sum(sum(nfloe));
subplot(1,5,j)
pcolor(hcat,Fsize,nfloe)
shading interp
colormap(cmap)
hold on
alphamap(alphmap)
xlabel('thickness')
ylabel('floe size')
zlabel('Distance')
end
clear all;
close all;
fig=figure;
fig.Units='inches';
fig.Position=[7.0729 5.6146 5.8333 4.3750];
fig.OuterPosition=[7.0729 5.6146 5.8333 5.3750];
fig.PaperUnits = 'inches';
fig.PaperPosition = [0 0 5.8333 4.3750];
fig.PaperPositionMode = 'manual';
fig.PaperOrientation='landscape';
fig.PaperSize=[5.8333 4.3750];
for i=1:100
dossier=['thickhs',num2str(i-1)];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
Dave=ncread(path,'Dave');
x=ncread(path,'x_axis');
Lmiz(i)=length(find(Dave>0 & Dave<396))*(x(2)-x(1));
end
Lmiz=reshape(Lmiz,10,10);
for i=1:10
plot(hs,Lmiz(:,i))
hold on
end
clear all;
close all;
for i=1:1000
dossier=['thickhstp',num2str(i-1)];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
Dave=ncread(path,'Dave');
x=ncread(path,'x_axis');
Lmiz(i)=length(find(Dave>0 & Dave<396))*(x(2)-x(1));
i
end
hs=[0.2 0.4 0.8 1 1.5 2 2.5 3 3.5 4];
tp=[6 7 8 9 10 11 12 13 14 15];
a=0;
b=0;
for i=1:10
a=b+1;
b=b+100;
Lmiz1(:,:,i)=reshape(Lmiz(a:b),10,10);
end
c=[2 5 8];
fig=figure;
fig.Units='inches';
fig.Position=[3.6250 7.7812 8.0312 2.2083];
fig.OuterPosition=[3.6250 7.7812 8.0312 3.2083];
fig.PaperUnits = 'inches';
fig.PaperPosition = [0 0 8.0312 2.2083];
fig.PaperPositionMode = 'manual';
fig.PaperOrientation='Portrait';
fig.PaperSize=[8.0312 2.2083];
cm=colormap('jet');
cm1(1,:)=[1 1 1];
for i=2:10
cm1(i,:)=[1-0.08*i 1-0.06*i 1];
end
cmap(1:10,:)=cm1;
cmap(11:59,:)=cm(16:end,:);
titre=['$h_{\mathrm{ice}}=$ 1m ';'$h_{\mathrm{ice}}=$ 2.5m';'$h_{\mathrm{ice}}$= 3.5m'];
titre=cellstr(titre);
for i=1:3
subplot(1,3,i)
pcolor(tp,hs,Lmiz1(:,:,c(i)))
hold on
contour(tp,hs,Lmiz1(:,:,c(i)),[50 100 200 300 450],'color','k','showtext','on','linestyle',':','Linewidth',1.5)
shading interp
xlabel('$T_p$[s]','interpreter','Latex','Fontsize',10)
if i==1
ylabel('$H_s$[m]','interpreter','Latex','Fontsize',10)
end
p=gca;
if i==1
p.Position(1)=p.Position(1)-0.05;
end
if i>1
p.Position(1)=p.Position(1)-0.07;
end
p.Position(2)=p.Position(2)+0.08;
p.Position(4)=p.Position(4)-0.13;
colormap(cmap)
caxis([0 500])
t=title(char(titre(i)));
set(t,'interpreter','Latex','Fontsize',10)
end
col=colorbar;
set(col,'position',[0.86 p.Position(2) 0.03 p.Position(end)])
ylabel(col,'$L_{MIZ}$ [km]','interpreter','Latex','Fontsize',10)
clear all;
close all;
for i=1:100
dossier=['pc_ec',num2str(i-1)];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
Dave=ncread(path,'Dave');
x=ncread(path,'x_axis');
Lmiz(i)=length(find(Dave>0 & Dave<396))*(x(2)-x(1));
i
end
Lmiz2=reshape(Lmiz,10,10);
Ec=[3e-5 4e-5 5e-5 6e-5 7e-5 8e-5 9e-5 1e-4 2e-4 3e-4];