Generating the right plot

조회 수: 5 (최근 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
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
Shruti Shivaramakrishnan 2016년 8월 31일
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?

Community Treasure Hunt

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

Start Hunting!

Translated by