1 view (last 30 days)

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.

Shruti Shivaramakrishnan
on 31 Aug 2016

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/300693-generating-the-right-plot#comment_387458

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/300693-generating-the-right-plot#comment_387458

Sign in to comment.