3 views (last 30 days)

Show older comments

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.

Image Analyst
on 27 Aug 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???

Shruti Shivaramakrishnan
on 31 Aug 2016

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

Start Hunting!