Generating the right plot

12 次查看（过去 30 天）
Kyle Broder
Kyle Broder 2016-8-24
回答： Shruti Shivaramakrishnan 2016-8-31
I want to plot the orbit of a galaxy in the Miyanomata-Nagai potential using Euler's method. The plot generated should look something like this
https://www.google.com.au/url?sa=i&rct=j&q=&esrc=s&source=images&cd=&cad=rja&uact=8&ved=0ahUKEwic98-QttnOAhWGJpQKHc-FAMYQjRwIBw&url=https%3A%2F%2Finspirehep.net%2Frecord%2F1236152&psig=AFQjCNG2iylnwrn0MJPSH7xV4O9A5ArWxA&ust=1472106464926420
I keep getting a straight line however.
My code is given by function
Eulersystem_MNmodel()
parsec = 3.08*10^18;
r_1 = 8.5*1000.0*parsec; % This converts our value of r into cm.
z_1 = 0.0;
theta_1 = 0.0; %Initial Value for Theta.
U_1 = 100.0*10^5; %Initial value for U in cm/sec
V = 156.972*10^5; %Proposed value of the radial velocity in cm/sec
W_1 = 1;
grav = 6.6720*10^-8; %Universal gravitational constant
amsun = 1.989*10^33; %Proposed Angular momentum of the sun
amg = 1.5d11*amsun;
gm = grav*amg; %Note this is constant
nsteps = 50000; %The number of steps
deltat = 5.0*10^11; %Measured in seconds
angmom = r_1*V; %The angular momentum
angmom2 = angmom^2.0; %The square of the angular momentum
E = -gm/r_1 + U_1*U_1/2 + angmom2/(2.0*r_1*r_1); %Total Energy of the system
time = 0.0;
for i=1:nsteps
r_1 = r_1 + deltat*U_1;
U_1 = U_1 + deltat*(-gm*r_1/(r_1^2.0 + (1+sqrt(z_1^2.0+1)^2.0)^1.5));
z_1 = z_1 + deltat*W_1;
W_1 = W_1 + deltat*(gm*z_1*(1+sqrt(z_1^2.0+1))/(sqrt(z_1^2.0+1)*(r_1^2.0+(1+sqrt(z_1^2.0+1))^2.0)^1.5));
theta_1 = theta_1 + deltat*(angmom/(r_1^2.0));
E = -gm/r_1+U_1/2.0+angmom2/(2.0*r_1*r_1);
ecc = (1.0 + (2.0*E*angmom2)/(gm^2.0))^0.5;
time1(i) = time;
time = time+deltat;
r(i) = r_1;
z(i) = z_1;
end
figure()
plot(r,z)
The code runs fine, but I'm getting a straight line rather than the spiral seen in the linked picture.
  1 个评论
Image Analyst
Image Analyst 2016-8-27
编辑：Image Analyst 2016-8-27
The link does not show any image. Please insert the image into the body of your question with the green and brown frame icon.
Also, don't use time as a variable. It's a built in function. And r and z are determined in the first 3 lines of your for loop, so why are those other lines calculating E, ecc, and theta_1 there if you don't use them???

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Shruti Shivaramakrishnan
The figure posted seems to represent a mesh of lines and would mostly require multiple lines to be plotted. I notice that the code consists of only the single plot command which may not be successful in plotting a complete mesh. Could there be any other variants that need to be accounted for while plotting?

类别

在 Help Center 和 File Exchange 中查找有关 Graphics Performance 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by