data_treatment.m 3.38 KB
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')
%