Generating the right plot

5 visualizzazioni (ultimi 30 giorni)
Kyle Broder
Kyle Broder il 24 Ago 2016
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 Commento
Image Analyst
Image Analyst il 27 Ago 2016
Modificato: Image Analyst il 27 Ago 2016
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???

Accedi per commentare.

Risposte (1)

Shruti Shivaramakrishnan
Shruti Shivaramakrishnan il 31 Ago 2016
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