Generating the right plot
1 回表示 (過去 30 日間)
古いコメントを表示
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
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
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?
0 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!