Prev: remove duplicate rows in sparse matrix
Next: Number Plate Extraction Please HELP
From: Dave on 10 Jul 2010 10:23
Dear all,

Can anybody suggest a way to animate the plots at each j step for my code. I have tried to use the set(obj_handle,'YDataSource','varname') but cannot get it to work


% The Beam Scheme

% clear all variables from the memory
clear all

% close all the figures
close all

load('Analyticaldata.mat')

% define the physical constants
omega=25; % angular velocity at omega(r) (km /sec /kpc)
kappa=31.3; % epicyclic frequency (km /sec /kpc)
omegapattern=13.5; % angular pattern velocity (km /sec /kpc)
c=8.56; % speed of sound (km /sec)
r=10; % radius of circular flow (kpc)
alpha=0.11667; % angle of spiral inclination (radians)
A=72.92; % amplitude (km /sec)^2

u0=alpha*r*(omega - omegapattern); % equilibrium velocity u
v0= r*(omega - omegapattern); % equilibrium velocity v

% define the numerical constants
N=64; % number of grid boxes around spiral phase
n=1200; % number of time-steps to run model
dt=0.000964416667; % time-step (unit=hours)
totaltime=dt*n; % total time
epsilon=0.0894427191; % averaging term small non-vanishing bias of order O((deta)^3)

deta=(pi*alpha*r)/N; % step in eta-direction

% declare arrays to store the values
Fpos=zeros(1,N,3);
Fneg=zeros(1,N,3);
Phi=zeros(1,N,3);

U=zeros(1,N,3);
F=zeros(1,N,3);
H=zeros(1,N,3);

Unew=zeros(1,N,3);
Fnew=zeros(1,N,3);
Hnew=zeros(1,N,3);

rho=zeros(1,N);
u=zeros(1,N);
v=zeros(1,N);

% fill in the initial values of the arrays
eta=0:deta:(N-1)*deta;

% uniform initial values
rho(1,:)=datarho(1,:);
u(1,:)=datau(1,:);
v(1,:)=datav(1,:);

% define scheme

for j=1:n

% perform n+1/2 source term advance

% rho(t_n + tau, eta_i)=rho

for i=1:N
tau=0.5*dt;
u(1,i)= u0 + (u(1,i) - u0).*cos(kappa.*tau) + ((2.*omega)./(kappa)).*(v(1,i) - v0 + (A./(alpha.*r.*omega)).*sin((2.*eta(i))./(alpha.*r))).*sin(kappa.*tau);
v(1,i)= v0 - (A./(alpha.*r.*omega)).*sin((2.*eta(i))./(alpha.*r)) + (v(1,i) - v0 + (A./(alpha.*r.*omega)).*sin((2.*eta(i))./(alpha.*r))).*cos(kappa.*tau) - (kappa./(2.*omega)).*(u(1,i) - u0).*sin(kappa.*tau);
end

% Evaluate F(U)
for i=1:N
F(1,i,1)=rho(1,i).*u(1,i); F(1,i,2)=rho(1,i).*(u(1,i).^2 + c.^2); F(1,i,3)=rho(1,i).*u(1,i).*v(1,i);
end

% evaluate F+ at i and i+1
for i=1:N
if u(1,i) >= c*sqrt(3);
Fpos(1,i,:)=F(1,i,:);
elseif 0 <= u(1,i) < c*sqrt(3);
Fpos(1,i,1)=0.16666667.*rho(1,i).*(5.*u(1,i) + c*sqrt(3));
Fpos(1,i,2)=0.16666667.*rho(1,i).*(4.*(u(1,i).^2) + (u(1,i) + c*sqrt(3)).^2);
Fpos(1,i,3)=0.16666667.*rho(1,i).*v(1,i).*(5.*u(1,i) + c*sqrt(3));
elseif -c*sqrt(3) < u(1,i) < 0;
Fpos(1,i,1)=0.16666667.*rho(1,i).*(u(1,i) + c*sqrt(3));
Fpos(1,i,2)=0.16666667.*rho(1,i).*(u(1,i) + c*sqrt(3)).^2;
Fpos(1,i,3)=0.16666667.*rho(1,i).*v(1,i).*(u(1,i) + c*sqrt(3));
elseif u(1,i) <= -c*sqrt(3);
Fpos(1,i,:)=0;
end
end

% evaluate F- = F+ - F(U)
for i=1:N
Fneg(1,i,:)=Fpos(1,i,:) - F(1,i,:);
end

for i=1:N-1
Phi(1,i,:)=Fpos(1,i,:) - Fneg(1,i+1,:);
end

% i=N
Phi(1,N,:)=Fpos(1,N,:) - Fneg(1,1,:);

for i=1:N
% define vectors of conserved quantities
U(1,i,1)=rho(1,i); U(1,i,2)=rho(1,i).*u(1,i); U(1,i,3)=rho(1,i).*v(1,i);
end

% i=1
Unew(1,1,:)=U(1,1,:) - (dt./deta).*(Phi(1,1,:) - Phi(1,N,:));
for i=2:N
Unew(1,i,:)=U(1,i,:) - (dt./deta).*(Phi(1,i,:) - Phi(1,i-1,:));
end

for i=1:N
% return individual vector quantities
rho(1,i)=Unew(1,i,1); u(1,i)=Unew(1,i,2)./rho(1,i); v(1,i)=Unew(1,i,3)./rho(1,i);
end

for i=1:N
tau=0.5*dt;
u(1,i)=u0 + (u(1,i) - u0).*cos(kappa.*tau) + ((2.*omega)./(kappa)).*(v(1,i) - v0 + (A./(alpha.*r.*omega)).*sin((2.*eta(i))./(alpha.*r))).*sin(kappa.*tau);
v(1,i)=v0 - (A./(alpha.*r.*omega)).*sin((2.*eta(i))./(alpha.*r)) + (v(1,i) - v0 + (A./(alpha.*r.*omega)).*sin((2.*eta(i))./(alpha.*r))).*cos(kappa.*tau) - (kappa./(2.*omega)).*(u(1,i) - u0).*sin(kappa.*tau);
end

end

% RMSE calculation
RMSE=zeros(1,4);

RMSE(1,1)=sqrt(sum((datarho-rho).^2)/(length(datarho)-1));
RMSE(1,2)=sqrt(sum((datau-u).^2)/(length(datau)-1));
RMSE(1,3)=sqrt(sum(((datarho.*datau)-(rho.*u)).^2)/(length((datarho.*datau))-1));
RMSE(1,4)=sqrt(sum((datav-v).^2)/(length(datav)-1));

% plot results

scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4) 400 scrsz(4)])

subplot(3,1,1), plot(eta, v, 'ko','MarkerSize',4)
hold on
plot(Analyticaleta, Analyticalv, 'k', 'linewidth', 1) % these two vectors are both 1x1000
xlim([0 3.6])
ylim([100 140])
ylabel('V')
set(gca,'xtick',0:0.3:3.6)
set(gca,'ytick',100:10:140)
set(gca,'XTickLabel',{' '})

title(['ITER=',num2str(2*n),' TIME=',num2str(totaltime)])

subplot(3,1,2), plot(eta, u,'ko','MarkerSize',4)
hold on
plot(Analyticaleta, Analyticalu, 'k', 'linewidth', 1)
xlim([0 3.6])
ylim([0 40])
ylabel('U')
set(gca,'xtick',0:0.3:3.6)
set(gca,'ytick',0:10:40)
set(gca,'XTickLabel',{' '})

subplot(3,1,3), plot(eta, rho,'ko','MarkerSize',4)
hold on
plot(Analyticaleta, Analyticalrho, 'k', 'linewidth', 1)
xlim([0 3.6])
xlabel('SPIRAL PHASE')
ylim([0 4])
ylabel('RHO')
set(gca,'xtick',0:0.3:3.6)
set(gca,'ytick',0:1:4)
set(gca,'XTickLabel',{'0',' ',' ','90',' ',' ','180',' ',' ','270',' ',' ','360'})
From: us on 10 Jul 2010 15:39
"Dave " <Davefulton(a)rocketmail.com> wrote in message <i19vo8$5jp$1(a)fred.mathworks.com>...
> Dear all,
>
> Can anybody suggest a way to animate the plots at each j step for my code. I have tried to use the set(obj_handle,'YDataSource','varname') but cannot get it to work

in your slew of code, which - by the way - is not helpful to demonstrate this problem, you do not show how you use the YDATASOURCE/REFRESH syntax...
produce a small(!) skeleton and come back to CSSM if you don't succeed...

us
 | 
Pages: 1
Prev: remove duplicate rows in sparse matrix
Next: Number Plate Extraction Please HELP