Prev: Defining function in Mupad
Next: Problems with RMS and IFFT amplitude
From: Dave on 25 Jun 2010 20:58
Dear Matlab gurus,

I am trying to plot two data series on the same plot but with a different number of elements. In fact it is the same solution but at a higher resolution.


% MS2

% clear all variables from the memory
% clear all

% close all the figures
close all

% 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=60; % number of grid boxes around spiral phase
n=1200; % number of time-steps to run model
dt=0.00099725; % time-step (unit=hours)

b=1.0; % dissipative term coefficient

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

% declare arrays to store the values
U=zeros(1,N,3);
F=zeros(1,N,3);
H=zeros(1,N,3);

D=zeros(1,N,3);

Ubar=zeros(1,N,3);
Fbar=zeros(1,N,3);
Hbar=zeros(1,N,3);

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

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

rhobar=zeros(1,N);
ubar=zeros(1,N);
vbar=zeros(1,N);

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

% uniform initial values
rho=ones(1,N);
u(1,:)=u0;
v(1,:)=v0;

% define scheme
% begin time loop
for j=1:n
%------------------------------------------------------------------------forward corrector on ODD steps---------------------------------------------------------------------------
% begin space loop
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);
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);
H(1,i,1)=0; H(1,i,2)=2.*omega.*(v(1,i) - v0).*rho(1,i) + (2./(alpha.*r)).*rho(1,i).*A.*sin(((2.*eta(i))./(alpha.*r))); H(1,i,3)= -((kappa.^2)./(2.*omega)).*(u(1,i) - u0).*rho(1,i);
end

% i=1
D(:,1,:)=b.*(dt./deta).*((abs(u(:,2) - u(:,1))).*(U(:,2,:) - U(:,1,:)) - (abs(u(1,1) - u(1,N))).*(U(:,1,:) - U(:,N,:)));

for i=2:N-1
% define additional dissipative term
D(:,i,:)=b.*(dt./deta).*((abs(u(:,i+1) - u(:,i))).*(U(:,i+1,:) - U(:,i,:)) - (abs(u(1,i) - u(1,i-1))).*(U(:,i,:) - U(:,i-1,:)));
end

% i=N
D(:,N,:)=b.*(dt./deta).*((abs(u(:,1) - u(:,N))).*(U(:,1,:) - U(:,N,:)) - (abs(u(1,N) - u(1,N-1))).*(U(:,N,:) - U(:,N-1,:)));

for i=1:N-1
% perform backward predictor step (odd)
Ubar(:,i,:)=U(:,i,:) - (dt./deta).*(F(:,i+1,:) - F(:,i,:)) + dt.*H(:,i,:);
end

% i=N
Ubar(:,N,:)=U(:,N,:) - (dt./deta).*(F(:,N,:) - F(:,1,:)) + dt.*H(:,N,:);

for i=1:N
% return individual vector quantities
rhobar(1,i)=Ubar(1,i,1); ubar(1,i)=Ubar(1,i,2)./rhobar(1,i); vbar(1,i)=Ubar(1,i,3)./rhobar(1,i);
end

for i=1:N
% update vectors of conserved values
Fbar(1,i,1)=rhobar(1,i).*ubar(1,i); Fbar(1,i,2)=rhobar(1,i).*(ubar(1,i).^2 + c.^2); Fbar(1,i,3)=rhobar(1,i).*ubar(1,i).*vbar(1,i);
Hbar(1,i,1)=0; Hbar(1,i,2)=2.*omega.*(vbar(1,i) - v0).*rhobar(1,i) + (2./(alpha.*r)).*rhobar(1,i).*A.*sin(((2.*eta(i))./(alpha.*r))); Hbar(1,i,3)= -((kappa.^2)./(2.*omega)).*(ubar(1,i) - u0).*rhobar(1,i);
end

% i=1
Unew(1,1,:)=0.5.*(U(1,1,:) + Ubar(1,1,:) - (dt./deta).*(Fbar(1,1,:) - Fbar(1,N,:)) + dt.*Hbar(1,1,:)) + D(1,1,:);

for i=2:N
% perform forward corrector step with dissipative term
Unew(1,i,:)=0.5.*(U(1,i,:) + Ubar(1,i,:) - (dt./deta).*(Fbar(1,i,:) - Fbar(1,i-1,:)) + dt.*Hbar(1,i,:)) + D(1,i,:);
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

%------------------------------------------------------------------------backward corrector on EVEN steps--------------------------------------------------------------------------
% begin space loop
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);
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);
H(1,i,1)=0; H(1,i,2)=2.*omega.*(v(1,i) - v0).*rho(1,i) + (2./(alpha.*r)).*rho(1,i).*A.*sin(((2.*eta(i))./(alpha.*r))); H(1,i,3)= -((kappa.^2)./(2.*omega)).*(u(1,i) - u0).*rho(1,i);
end

% i=1
D(:,1,:)=b.*(dt./deta).*((abs(u(:,2) - u(:,1))).*(U(:,2,:) - U(:,1,:)) - (abs(u(1,1) - u(1,N))).*(U(:,1,:) - U(:,N,:)));

for i=2:N-1
% define additional dissipative term
D(:,i,:)=b.*(dt./deta).*((abs(u(:,i+1) - u(:,i))).*(U(:,i+1,:) - U(:,i,:)) - (abs(u(1,i) - u(1,i-1))).*(U(:,i,:) - U(:,i-1,:)));
end

% i=N
D(:,N,:)=b.*(dt./deta).*((abs(u(:,1) - u(:,N))).*(U(:,1,:) - U(:,N,:)) - (abs(u(1,N) - u(1,N-1))).*(U(:,N,:) - U(:,N-1,:)));

% i=1
Ubar(:,1,:)=U(:,1,:) - (dt./deta).*(F(:,1,:) - F(:,N,:)) + dt.*H(:,1,:);

for i=2:N
% perform forward predictor step (odd)
Ubar(:,i,:)=U(:,i,:) - (dt./deta).*(F(:,i,:) - F(:,i-1,:)) + dt.*H(:,i,:);
end

for i=1:N
% return individual vector quantities
rhobar(1,i)=Ubar(1,i,1); ubar(1,i)=Ubar(1,i,2)./rhobar(1,i); vbar(1,i)=Ubar(1,i,3)./rhobar(1,i);
end

for i=1:N
% update vectors of conserved values
Fbar(1,i,1)=rhobar(1,i).*ubar(1,i); Fbar(1,i,2)=rhobar(1,i).*(ubar(1,i).^2 + c.^2); Fbar(1,i,3)=rhobar(1,i).*ubar(1,i).*vbar(1,i);
Hbar(1,i,1)=0; Hbar(1,i,2)=2.*omega.*(vbar(1,i) - v0).*rhobar(1,i) + (2./(alpha.*r)).*rhobar(1,i).*A.*sin(((2.*eta(i))./(alpha.*r))); Hbar(1,i,3)= -((kappa.^2)./(2.*omega)).*(ubar(1,i) - u0).*rhobar(1,i);
end

for i=1:N-1
% perform backward corrector step with dissipative term
Unew(1,i,:)=0.5.*(U(1,i,:) + Ubar(1,i,:) - (dt./deta).*(Fbar(1,i+1,:) - Fbar(1,i,:)) + dt.*Hbar(1,i,:)) + D(1,i,:);
end

% i=N
Unew(1,N,:)=0.5.*(U(1,N,:) + Ubar(1,N,:) - (dt./deta).*(Fbar(1,1,:) - Fbar(1,N,:)) + dt.*Hbar(1,N,:)) + D(1,N,:);

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
end

% 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(eta, Analyticalv) 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',{' '})

SUBPLOT(3,1,2), plot(eta, u,'ko','MarkerSize',4)
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)
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'})

Sorry for the long code. The section I am trying to write is in the % plot results
From: dpb on 25 Jun 2010 21:09
Dave wrote:
....

> I am trying to plot two data series on the same plot but with a
> different number of elements. In fact it is the same solution but at a
> higher resolution.
....

....[_WAY_ too much code elided]...

And the problem is????...

At least say what's the problem you've run into; better give a very
short sample code that demonstrates it instead of all the generation
code that has no bearing on the plot...

cs-sm'ers are generally busy folk--willing to try to help but make it
easy on 'em, if you really want somebody to spend some time.

--
 | 
Pages: 1
Prev: Defining function in Mupad
Next: Problems with RMS and IFFT amplitude