How to draw evolution of concentration in 2D [solved]

Hi all,

I have coded a small program in order to compute the evolution of concentration of a plume of smoke and would like to plot/shade it as follow (I did it with Ferret but would like to know how to do it with Matlab) :

My matrix is defined along x & z axis with concentration(x,z) defined as a sinusoidal function

Thanks for all,

Florian

Answer 1

[Image Analyst]

Why not make save be a 3D array instead of a cell array? Why not just use imshow() or image() to display each 2D image of "save" right after you created it.

[Florian]

That's it ! imshow() is the function I was looking for (with some more parameters to add) Thank you so much for your help Best Florian

Answer 2

[Walter Roberson]

Perhaps surf() would be appropriate

[Florian]

Hi,

surf is not appropriated actually, it doesn't work so far ... I copy past my code, it may help :

Many thanks for your kind help

 deltat=0.2*deltax./umax;
vitu=umin+(umax-umin).*zzz(:)/zmax;
alphax=(deltat./deltax).*vitu(:);

 qqq(:,:)=zeros(imax,kmax); % we initialize the grid with zeros
qqq(imax/6,3*kmax/4)=1; % except in one point which corresponds to the source

 save=[];
        %========================================================================
 for it=1:nt % TIME LOOP
%=========================================================================

 temps=it.*deltat;
 qancien(:,:) = qqq(:,:); %copy of the old concentration

 % Advection
     for k=1:kmax
         vent = umin + (umax-umin).*k./kmax;
         alpha = vent*deltat./deltax;

 %Horizontal advection
         for i=2:imax-1
               qqq(i,k) = qancien(i,k).*(1-alpha)+qancien(i-1,k).*alpha;
         end
     end

     qqq(imax./6 , 3.*kmax./4)=1;
     qancien(:,:)=qqq(:,:); % we save the matrix for the vertical advection which will use it

%Advection verticale  
    for i=2:imax-1  
        for k=2:kmax-1
            qqq(i,k) = qancien(i,k)+deltat.*Kz.*(qancien(i,k+1)-2.*qancien(i,k)+ qancien(i,k-1))./(2.*deltaz.*deltax); 

           end
     end

     qqq(imax/6 , 3.*kmax./4)=1;

     save{it}(:,:)=qqq(:,:); %copy the qqq matrix at the time step 'it' in order to plot it later at different time steps

 end

[José-Luis]

And why is it not appropriate? You realize we can't run your code or see the error messages that you get?

[Florian]

There is no error message, it's just that I don't know how to "plot" the result, but I agree with you, I should have written the code. Please find it below :

 %========================================================================
% DECLARATIONS
%========================================================================
clear
close

 %========================================================================
% INITIALISATIONS
%========================================================================

 %------------------------------------------------------------------------
% Initialize constants
%------------------------------------------------------------------------
imax=300;
kmax=100;
nt=2001;
xmax=3000.;
zmax=100.;
umin=2.;
umax=10.;
x0=xmax./5.;
wmax=5.;
k1=30.;
k2=10.;
tau1=50.;
tau2=66.;
Kz=25.; %diffusion coefficient 
pi=2.*asin(1.) ;

 %------------------------------------------------------------------------
% Initialize grid
%------------------------------------------------------------------------
deltax=xmax./imax; % step in x
deltaz=zmax./kmax; % step in z

 for i=1:imax
   xxx(i)=i.*deltax; 
end

 for k=1:kmax
   zzz(k)=(k-0.5).*deltaz; 
end

 %------------------------------------------------------------------------
% Initialize tracer
%------------------------------------------------------------------------
for i=1:imax
   qinitial(i,:)=exp( -(25.*(xxx(i)-x0)./xmax).^2 ); 
end

 %------------------------------------------------------------------------
% alphax(z) = delta t * U(z) /deltax (sans dimension)
%------------------------------------------------------------------------
deltat=0.2*deltax./umax;
vitu=umin+(umax-umin).*zzz(:)/zmax;
alphax=(deltat./deltax).*vitu(:); 

 deltat=0.2*deltax./umax;
vitu=umin+(umax-umin).*zzz(:)/zmax;
alphax=(deltat./deltax).*vitu(:);

 qqq(:,:)=zeros(imax,kmax); % we initialize the grid with zeros
qqq(imax/6,3*kmax/4)=1; % except in one point which corresponds to the source

 save=[];
        %========================================================================
 for it=1:nt % TIME LOOP
%=========================================================================

 temps=it.*deltat;
 qancien(:,:) = qqq(:,:); %copy of the old concentration

 % Advection
     for k=1:kmax
         vent = umin + (umax-umin).*k./kmax;
         alpha = vent*deltat./deltax;

 %Horizontal advection
         for i=2:imax-1
               qqq(i,k) = qancien(i,k).*(1-alpha)+qancien(i-1,k).*alpha;
         end
     end

     qqq(imax./6 , 3.*kmax./4)=1;
     qancien(:,:)=qqq(:,:); % we save the matrix for the vertical advection which will use it

 %Advection verticale  
    for i=2:imax-1  
        for k=2:kmax-1
            qqq(i,k) = qancien(i,k)+deltat.*Kz.*(qancien(i,k+1)-2.*qancien(i,k)+ qancien(i,k-1))./(2.*deltaz.*deltax); 

           end
     end

     qqq(imax/6 , 3.*kmax./4)=1;

     save{it}(:,:)=qqq(:,:); %copy the qqq matrix at the time step 'it' in order to plot it later at different time steps

 end

[Image Analyst]

What is the thing you want to display/plot? qqq? save? qancien?

By the way, DON'T USE save AS THE NAME OF A VARIABLE. That is the name of a very important function, and you just destroyed it (not permanently, just in this m-file).

[Florian]

I'd like to display save by the way ... Considering that I've changed its name