Generating the right plot

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.

[Image Analyst]

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???

Answer 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?