Archive for the 'gravity' Category

13
Nov
07

Galactic Gravity Simulations

Here is a simulation of a galaxy using Newton’s law of universal gravitation. Click here to see a 3D rotatable model of the Milky Way Galaxy.

(* runtime: 8 seconds *)
n = 275; G = 1; SeedRandom[0]; image = Table[0, {n}, {n}]; nstar = 500;dt = 0.001; theta = Pi/4; Trot = {{1, 0, 0}, {0, Cos[theta], Sin[theta]}, {0, -Sin[theta], 1}};
stars = Table[r = 0.5 Random[]; theta = 2Pi Random[]; {Trot.{r Cos[theta], r Sin[theta],0.1(2Random[] - 1)Exp[-4r]} + {0, 0, 1}, Sqrt[G/r]{-Sin[theta], Cos[theta], 0}}, {nstar}];
Do[image *= 0.9; stars = Map[({r, v} = #; v -= G r dt/(r.r)^1.5; r += v dt; {x, y, z} = r; {j, i} = Round[n({x, y}/z + 0.5)]; If[0 < i <= n && 0 < j <= n, image[[i, j]] += 0.5/z]; {r, v}) &, stars]; ListDensityPlot[image, Mesh -> False, Frame -> False,PlotRange -> {0, 1}], {50}];

Links

01
Dec
05

Halo Orbit

In 1772, Joseph Lagrange showed that a small satelite moving in the Earth-Moon system can have 5 balance points, called Lagrange points. Here is a small halo orbit around Lagrange point L4. This is interesting because L4 is off to the side, so it looks like the satelite is orbiting nothing! The velocity components have been represented by the third axis and color.
(* runtime = 6 seconds *)
Clear[x, y, t, p, v]; mu = 0.0123001; x1 = -mu; x2 = 1 - mu; r1 = Sqrt[(x - x1)^2 + y^2]; r2 = Sqrt[(x - x2)^2 + y^2];
{x0, y0} = {x, y} /. Last[NSolve[{-x == -x2(x - x1)/r1^3 + x1(x - x2)/r2^3, -y == -(x2/r1^3 - x1/r2^3) y}, {x, y}]];
JacobianMatrix[p_, q_] := Outer[D, p, q];
vlist = Eigenvectors[JacobianMatrix[{xdot, -x2(x - x1)/r1^3 + x1(x - x2)/r2^3 + 2ydot + x, ydot, -(x2/r1^3 - x1/r2^3)y - 2xdot + y}, {x, xdot, y, ydot}] /. {x -> x0, xdot -> 0, y -> y0, ydot -> 0}];
r1 = Sqrt[(x[t] - x1)^2 + y[t]^2]; r2 = Sqrt[(x[t] - x2)^2 + y[t]^2];
tmax = 355; p0 = {x0, y0, 0, 0} + Re[0.5 vlist[[1]] + vlist[[3]]]/3.84400*^5;
p[t_] = {x[t], y[t]} /. NDSolve[{x''[t] - 2y'[t] - x[t] == -x2(x[t] - x1)/r1^3 + x1(x[t] - x2)/r2^3, y''[t] + 2x'[t] - y[t] == -(x2/r1^3 - x1/r2^3)y[t], x[0] == p0[[1]], y[0] == p0[[2]], x'[0] == p0[[3]], y'[0] == p0[[4]]}, {x[t], y[t]}, {t, 0, tmax}, MaxSteps -> 100000, AccuracyGoal -> 10][[1]];
v[t_] = D[p[t], t];
ParametricPlot3D[Append[p[t], v[t][[1]]], {t, 0, tmax}, Compiled -> False, PlotPoints -> 1000, BoxRatios -> {1, 1, 1}, ViewPoint -> {0, 10, 0}]

Links

05
Feb
05

Gravitational Lensing of a Black Hole

Gravitational Lensing of a Black Hole : According to Einstein’s Theory of Relativity, the intense gravity of a black hole can bend light into a circle called an Einstein Ring. See also Event horizon, Schwarzschild radius, photon sphere. The following code is only an approximation:
(* runtime: 30 seconds *)
Clear[stars]; SeedRandom[535]; stars[{x_, y_}] = Sum[Exp[-500((x - Random[])^2 + (y - Random[])^2)/Random[]^2], {10}];
DensityPlot[Module[{r = x^2 + y^2}, If[r 275, PlotRange -> {0, 1}, Mesh -> False, Frame -> False]

Links




Welcome !

You will find here some of my favorite hobbies and interests, especially science and art.

I hope you enjoy it!

Subscribe to the RSS feed to stay informed when I publish something new here.

I would love to hear from you! Please feel free to send me an email : bugman123-at-gmail-dot-com

Archives

Blog Stats

  • 544,672 hits