 |
Home |
 |
Physics Simulations |
 |
Engineering |
 |
Fractals |
 |
Math Artwork |
 |
Rendered Artwork |
 |
Hand-made Artwork |
 |
Bug Collection |
 |
Programming |
 |
High Voltage |
 |
Holography |
 |
Physics Experiments |
 |
Legos |
 |
Claymation |
 |
The Bible |
 |
Links |
Here are some physics simulations and artwork along with some Mathematica code. See also my Engineering page.
Magnetic Field of a Solenoid - new version: POV-Ray 3.6.1, 5/17/06
old version: calculated in Mathematica 4.2, rendered in TrueSpace 4.3, 1/19/03
 This magnetic field was approximated by a superposition of 2D point sources using the Biot-Savart Law. This entry received 2nd place in the Art in Science Contest. You can also see an older version of this picture on Jeff Bryant’s Mathematica visualization site. Click here to download the POV-Ray code for this image. See also my magnetic field representations for a motor, Tesla coil, and horseshoe magnets. (* runtime: 12 seconds *) plist = Table[{(4 i - 26)/6, -(-1)^i}, {i, 1, 12}]; r[{xi_, yi_}] := Sqrt[(x - xi)^2 + (y - yi)^2]; DensityPlot[2Sqrt[((Plus @@ Map[#[[2]](x - #[[1]])/r[#]^2 &, plist])^2 + (Plus @@ Map[-#[[2]](y - #[[2]])/r[#]^2 &, plist])^2)] + Cos[18.8Plus @@ Map[#[[2]]/r[#] &, plist]] + 1, {x, -6, 6}, {y, -3, 3}, Mesh -> False, Frame -> False, PlotRange -> {0, 10}, PlotPoints -> {275, 138}, AspectRatio -> 1/2] |
Water Ripple Simulation - calculated in Mathematica 4.2, 12/2/04; rendered in POV-Ray 3.6.1, 7/29/05
 Here is a simple water ripple simulation showing single slit wave diffraction. The following Mathematica code solves the wave equation with damping using the finite difference method. You can read more about this algorithm on Hugo Elias’ website. (Note: technically this simulation should use Neumann boundary conditions but I decided it was simplier to demonstrate using Dirichlet boundary conditions). (* runtime: 18 seconds, c is the wave speed and b is a damping factor *) n = 64; c = 1; b = 5; dx = 1.0/(n - 1); Courant = Sqrt[2.0]/2;dt = Courant dx/c; z1 = z2 = Table[0.0, {n}, {n}]; Do[{z1, z2} = {z2, z1}; z1[[n/2, n/4]] = Sin[16Pi t]; Do[If[0.45 < i/n < 0.55 || ! (0.48 < j/n < 0.52), z2[[i, j]] = (2(1 - 2Courant^2)z1[[i, j]] + Courant^2(z1[[i - 1, j]] + z1[[i + 1, j]] + z1[[i, j - 1]] + z1[[i, j + 1]]))/(1 + b dt) - z2[[i, j]]], {i, 2, n - 1}, {j, 2, n - 1}]; ListPlot3D[z1, Mesh -> False, PlotRange -> {-1, 1}], {t, 0, 1, dt}]; Links Water Figures - beautiful high-speed camera splashes by Fotoopa LSSM - POV-Ray animations using this same technique, by Tim Nikias Wenclawiak POV-Ray Ripples - another one by Rune Johansen Water Ripples - Delphi program by Jan Horn Waves - C++ program by Maciej Matyka Fluid Animations - amazing animations by Ron Fedkiw, with Eran Guendelman, Andrew Selle, Frank Losasso, et al. Free Surface Fluid Simulations - impressive simulations using the Lattice-Boltzmann method with level sets, by Nils Thuerey, author of Blender’s fluid package Fluid v1.0 - C++ program for fluid surfaces by Maciej Matyka, uses Marker And Cell (MAC) method ball in water - nice 3D POV-Ray animation by Fidos Splash - analytical Mathematica plot by Dr. Jim Swift Smoke Animations - by Jos Stam, Henrik Jensen, Ron Fedkiw, et al. Fire Animations - by Duc Nguyen, Henrik Jensen, Ron Fedkiw |
Water Caustics - Mathematica 4.2, POV-Ray 3.1, 4/22/06
This water-like surface was generated in Mathematica using frequency filtered random noise, and then it was raytraced in POV-Ray and water caustics were added using Henrik Jensen’s photon mapping technique. The smooth_triangle command was used for phong normal interpolation of the water’s surface.Link: Water Caustics Generator - C++ program by Kjell Andersson, uses Heron’s Formula to calculate the caustics |
Leap-Frogging Bubble Rings - Mathematica 4.2, POV-Ray 3.6.1, 12/31/05
 Here is a simple technique for modeling vortices assuming inviscid incompressible potential flow (irrotational). This technique was inspired from Kerry Mitchell’s paper based on Ultra Fractal. You can also see this here on Dynamical Systems’ Picture Gallery. Here is some Mathematica code to plot the entrained fluid: (* runtime: 25 seconds, increase n for better resolution *) tmax = 0.85; rcore = 0.1; Klist = {1, -1, 1, -1}; zlist0 = {-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}; m = Length[zlist0]; v[K_, z_, z0_] := Module[{r2 = Abs[z - z0]^2}, I K(z0 - z)/r2(1 - Exp[-r2/rcore^2])]; zlist = NDSolve[Flatten[Table[{Subscript[z,i]'[t] == Sum[If[i == j, 0, v[Klist[[j]], Subscript[z, i][t], Subscript[z, j][t]]], {j, 1, m}], Subscript[z, i][0] ==zlist0[[i]]}, {i, 1, m}]], Table[Subscript[z, i][t], {i, 1, m}], {t, 0, tmax}][[1, All, 2]]; n = 23; image = Table[NDSolve[{z'[t] == Sum[v[Klist[[i]], z[t], zlist[[i]]], {i, 1, m}], z[tmax] == x + I y}, z[t], {t, 0, tmax}, MaxSteps -> 5000][[1, 1, 2]] /. t -> 0, {y, -1.125, 1.125, 2.25/n}, {x, -0.35, 1.9, 2.25/n}]; ListDensityPlot[Map[Sign[Im[#]]Arg[#] &, image, {2}], Mesh -> False, Frame -> False, AspectRatio -> Automatic] |
Here is some Mathematica code to numerically solve this using the 4th order Runge-Kutta method:(* runtime: 21 seconds, increase n for better resolution *) n = 23; tmax = 0.85; dt = tmax/100; rcore = 0.1; Klist = {1, -1, 1, -1}; zlist = {-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}; m = Length[zlist]; v[K_, z_, z0_] := Module[{r2 = Abs[z - z0]^2}, I K(z0 - z)/r2(1 - Exp[-r2/rcore^2])]; Runge[z_] := Module[{}, k1 = dt vtot[z]; k2 = dt vtot[z + k1/2]; k3 = dt vtot[z + k2/2]; k4 = dt vtot[z + k3]; (k1 + 2 k2 + 2 k3 + k4)/6]; vtot[zlist_] := Table[Sum[If[i == j, 0, v[Klist[[j]], zlist[[i]], zlist[[j]]]], {j, 1, m}], {i, 1, m}]; zlists = Transpose[Table[zlist += Runge[zlist], {t, 0, tmax, dt}]]; vtot[z_] := Plus @@ v[Klist, z, zlists[[All, t]]]; image = Table[z = x + I y; Do[z -= Runge[z], {t, 100, 1, -1}]; z, {y, -1.125, 1.125, 2.25/n}, {x, -0.35, 1.9, 2.25/n}]; ListDensityPlot[Map[Sign[Im[#]]Arg[#] &, image, {2}], Mesh -> False, Frame -> False, AspectRatio -> Automatic] Links My Educational Project - some basic potential flows using Matlab and Mathematica Vortex Movies - by Tee Tai Lim, I especially like the movies of leap-frogging and vortex ring collisions Bubble Ring Videos - blown by human, whale, dolphins. This is fun... try it some time! Leap-Frogging Bubble Rings - blown by scuba diver Bubble Rings - by David Whiteis Vortex rings blown out by Mt. Etna - same as a smoker puffing out an “O” ring Nuclear Bombs - mushroom clouds are giant vortex rings Birth of the Jellyfish - nice animation by Kerry Mitchell Smoke Ring - giant vortex ring at the Miramar Air Show Zero Blaster - fun toys that create fog vortex rings Ideal Flow Machine - Java applet by William Devenport Inviscid Flows - paper on incompressible, irrotational flows by Paul Hammerton |
2D Wave-Particle Simulations - Mathematica 4.2, 1/15/05; C++ version: 5/16/07
 The left animation shows the quantum mechanical probabilty distribution of an electron as it passes through two narrow slits. The electron itself is much smaller than its probability distribution cloud, which is dispersed over a large area, creating an interference pattern. As soon as the actual location of the electron itself is detected, then its probability distribution cloud will become small again. The color represents the complex phase. The following Mathematica code uses the Crank-Nicolson method with time-splitting to solve the Schrödinger wave equation: (* runtime: 10 seconds, V is the potential energy *) << LinearAlgebra`Tridiagonal`; n = 55; dx = 1.0/(n - 1); hbar = m = 1; dt = m dx^2/hbar; a = I hbar/dt; b = hbar^2/(4 m dx^2); blist = Table[b, {n - 1}]; psi = Table[Module[{p = {x - 0.25, y - 0.5}}, Exp[I {100, 0}.p - p.p/(2 0.14^2)]], {y, 0, 1, dx}, {x, 0, 1, dx}]; V = Table[If[0.48 < x < 0.52 && ! (0.4 < y < 0.475 || 0.525 < y < 0.6), 10000, 0], {y, 0, 1, dx}, {x, 0, 1, dx}]; Do[psi = Table[TridiagonalSolve[blist, a - 2b - V[[All, j]]/2, blist, (a + 2b + V[[All, j]]/2)psi[[All, j]] - b Table[psi[[i, Min[n, j + 1]]] + psi[[i, Max[1, j - 1]]], {i, 1, n}]], {j, 1, n}]; psi = Table[TridiagonalSolve[blist, a - 2b - V[[i,All]]/2, blist, (a + 2b + V[[i, All]]/2)psi[[All, i]] - b Table[psi[[j, Min[n, i + 1]]] + psi[[j, Max[1, i - 1]]], {j, 1, n}]], {i, 1, n}]; Show[Graphics[RasterArray[Map[(x = Min[Abs[#]^2, 1];Hue[Arg[#]/(2Pi), Min[1, 2(1 - x)], Min[1, 2x]]) &, psi + V, {2}]], AspectRatio -> 1]], {40}]; |
In case you’re wondering how TridiagonalSolve works:TridiagonalSolve[a_, b_, c_, r_] := Module[{n = Length[r], p, q}, p = q =Table[0, {n}]; {p[[1]], q[[1]]} = {c[[1]], r[[1]]}/b[[1]]; Do[{p[[i]], q[[i]]} = {If[i < n, c[[i]], 0], r[[i]] - a[[i - 1]]q[[i - 1]]}/(b[[i]] - a[[i - 1]]p[[i - 1]]), {i, 2, n}]; Do[q[[i]] -= p[[i]]q[[i + 1]], {i, n - 1, 1, -1}]; q]; Links Advanced Visual Quantum Mechanics - nice animations Fortran simulations - by Terry Robb Shroedinger simulator - C++ program by Maciej Matyka high energy quantum eigenfunctions in hyperbolic space - unexplained patterns by Dennis Hejhal |
Pool Simulation - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 7/9/05
 Here’s a fun one. The very instant when the cue ball breaks a tightly-packed triangle, the inner balls must rapidly bounce around until they spread out. (* runtime: 2 seconds *) Normalize[p_] := p/Sqrt[p.p]; Assoc[x_, data_] := Flatten[Select[data, #[[1]] == x &], 1]; r = 1.125; xmax = 50; ymax = 25; n = 16; balls = Prepend[Flatten[Table[{{xmax/2 + 1.05Sqrt[3]x, 1.05y}, {0, 0}}, {x, 4r, 0, -r}, {y, -x, x, 2r}], 1], {{-xmax/2 - 1.1r, 0}, {100, 0}}]; Step[dt0_] := Module[{dt = dt0, collisions = {}}, Do[{pa1, va} = balls[[i]]; pa2 = va dt0 + pa1; Do[If[j != i, {pb1, vb} = balls[[j]]; pb2 = vb dt0 + pb1; p1 = pa1 - pb1; p2 = pa2 - pb2; p = p1 - p2; a = p.p; b = -2p1.p; c = p1.p1 - (2r)^2; d = b^2 - 4a c; If[d >= 0, t = If[a == 0, If[b == 0, 1, -c/b], Min[(-b + {1, -1}Sqrt[d])/(2a)]]; If[0 <= t dt0 <= dt, pa = (1 - t)pa1 + t pa2; pb = (1 - t)pb1 + t pb2; normal = Normalize[pa - pb]; va -= (va - vb).normal normal; If[t dt0 < dt, dt =t dt0; collisions = {{i, va}}, collisions = Append[collisions, {i, va}]]]]], {j, 1, n}], {i, 1, n}]; balls = Table[{pa1, va1} = balls[[i]]; temp = Assoc[i, collisions]; {va1 dt + pa1, If[temp != {}, temp[[2]], va1]}, {i, 1, n}]; If[dt < dt0, Step[dt0 - dt]]]; Do[Show[Graphics[Map[Disk[#[[1]],r] &, balls], PlotRange -> {{-xmax, xmax}, {-ymax, ymax}}, AspectRatio -> 0.5]]; Step[0.1], {50}]; Link: Billiard Ball Simulation - includes angular momentum, by Jason Stewart |
3D Collisions - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 7/9/05
Here is a box of frictionless bouncy balls with gravity. The balls don’t rotate and never stop bouncing because there is no friction. This simulation uses the same basic code as the pool table simulation shown above.(* runtime : 52 seconds *) << Graphics`Shapes`; g = {0, 0, -9.80665}; xmax = 1; Normalize[p_] := p/Sqrt[p.p]; Distance[p1_, p2_] := Sqrt[(p2 - p1).(p2 - p1)]; Move[{p_, v_}, t_] := {g t^2/2 + v t + p,g t + v}; Interpolate[t_, x1_, x2_] := (1 - t)x1 + t x2; Step[dt_] := Module[{balls2 = Table[{pa1, va1} = balls[[i]]; {pa2, va2} = Move[{pa1, va1}, dt]; result = {pa2, va2}; tmin = Infinity; Do[If[j != i, {pb1, vb1} = balls[[j]]; {pb2, vb2} = Move[{pb1, vb1}, dt]; p1 = pa1 - pb1; p2 = pa2 - pb2; p = p1 - p2; a = p.p; b = -2p1.p;c = p1.p1 - (2r)^2; d = b^2 - 4a c; If[d >= 0, If[a == 0, t = If[b == 0, 1, -c/b], t = (-b + {1, -1}Sqrt[d])/(2a); t = If[0 <= Min[t] <= 1, Min[t], t[[If[0 <= t[[1]] <= 1, 1, 2]]]]]; If[0 <= t <= 1 &&t < tmin, tmin = t; pa = Interpolate[t, pa1, pa2]; pb =Interpolate[t, pb1, pb2]; va = Interpolate[t, va1, va2]; vb = Interpolate[t, vb1, vb2]; normal = Normalize[pa - pb]; va -= (va - vb).normal normal;result = Move[{pa, va}, (1 - t)dt]]]], {j, 1, n}];Do[If[va1[[dim]] != 0, Scan[Module[{t = (#(xmax - r) - pa1[[dim]])/(va1[[dim]]dt)}, If[0 <= t <= 1 && t < tmin, tmin = t;pa = Interpolate[t, pa1, pa2]; va = Interpolate[t, va1, va2]; va[[dim]] = -va[[dim]];result = Move[{pa, va}, (1 - t)dt]]] &, {-1, 1}]], {dim, 1, 3}]; result, {i, 1, n}]},collision = False; Do[pa2 = balls2[[i, 1]]; If[Min[pa2] < -(xmax - r) || Max[pa2] > xmax -r, collision = True]; Do[If[Distance[balls2[[j, 1]], pa2] < 2r, collision = True], {j, i + 1, n}], {i, 1, n}]; If[collision, Do[Step[dt/2], {2}], balls = balls2]]; r = 0.25; n = 20; balls = {}; SeedRandom[0]; Do[start = True; While[start || Or @@ Map[Distance[p, #[[1]]] < 2r &, balls], start = False; p = (xmax - r)Table[2 Random[] - 1, {3}]];balls = Append[balls, {p, {0, 0, 0}}], {n}]; Do[Step[0.05]; Show[Graphics3D[{EdgeForm[], Map[TranslateShape[Sphere[r, 10, 6], #[[1]]] &, balls]}, PlotRange -> xmax{{-1, 1}, {-1, 1}, {-1, 1}}]], {50}]; Collision Links Phun - real-time 2D sketched simulations, by Emil Ernerfeldt Collision Detection - similar program by Jason Stewart Rigid Body Animations - by Eran Guendelman, Ron Fedkiw, and Robert Bridson Rune’s Particle System - nice POV-Ray program by Rune Johansen Rigid Body Dynamics - Pixar notes by David Baraff Particle Fluid - Java applet using a simpler technique, by Chris Laurel |
Driven Cavity - calculated in Fortran 90, plotted in Mathematica 4.2, 11/2/05
Here is the flow inside a square box where the flow across the top of the box is moving to the right at Reynolds number 1000. This program uses the finite volume method to solve the Navier-Stokes equations assuming steady, incompressible, viscous, laminar flow. This was calculated on a 300×300 non-uniform grid and took 3 hours to run on my laptop. The animation shows the motion along the streamlines/pathlines. The following Mathematica code is not completely accurate at the boundaries, but it gives the basic idea:(* runtime: 20 minutes *) n = 50; dx = 1.0/(n - 1); RE = 1000; relax = 0.4; residu = residv = residp = 1; u = Table[If[j == 1 || i == 2 || i == n, 0, 1], {i, 1, n}, {j, 1, n}]; v = p = pstar = ap = an = as = aw = ae = sc = sp = du = dv = Table[0, {n}, {n}]; lisolv[i0_, j0_, phi0_] := Module[{phi = phi0, p = q =Table[0, {n}]}, Do[q[[j0 - 1]] = phi[[i, j0 - 1]]; Do[p[[j]] = an[[i, j]]/(ap[[i, j]] -as[[i, j]]p[[j - 1]]); q[[j]] = (as[[i, j]]q[[j - 1]] + ae[[i, j]] phi[[i + 1, j]] + aw[[i, j]]phi[[i - 1, j]] + sc[[i, j]])/(ap[[i, j]] - as[[i, j]]p[[j - 1]]), {j, j0, n - 1}]; Do[phi[[i,j]] = p[[j]]phi[[i, j + 1]] + q[[j]], {j, n - 1, j0, -1}], {i, i0, n - 1}]; phi]; calc[i0_, j0_] := Module[{}, {cn, cs, ce, cw} = 0.5 {v[[i,j + 1]] + v[[i0, j0 + 1]], v[[i, j]] + v[[i0, j0]], u[[i + 1, j]] + u[[i0 + 1, j0]], u[[i, j]] + u[[i0, j0]]}dx; cp = Max[0, cn - cs + ce - cw]; {an[[i, j]], as[[i, j]], ae[[i, j]], aw[[i, j]]} = Map[Max[0, 1 - 0.5 Abs[#]] &, RE{cn, cs, ce, cw}]/RE + Map[Max[0, #] &, {-cn, cs, -ce, cw}]; sp[[i, j]] = -cp; ap[[i, j]] = an[[i, j]] + as[[i, j]] + ae[[i, j]] + aw[[i, j]] - sp[[i, j]]]; While[Max[residu, residv, residp] > 0.001, residu = residv = residp = 0; Do[calc[i - 1, j]; sc[[i, j]] = cp u[[i, j]] + dx(p[[i - 1, j]] - p[[i, j]]); du[[i, j]] = relax dx/ap[[i, j]]; residu += Abs[an[[i, j]]u[[i, j + 1]] + as[[i, j]]u[[i, j - 1]] + ae[[i, j]]u[[i + 1, j]] + aw[[i, j]]u[[i - 1, j]] - ap[[i, j]]u[[i, j]] + sc[[i, j]]], {i, 3, n - 1}, {j, 2, n - 1}]; ap /= relax; sc += (1 - relax) ap u; u = lisolv[3, 2, u]; Do[calc[i, j - 1]; sc[[i, j]] = cp v[[i, j]] + dx(p[[i, j - 1]] - p[[i, j]]); dv[[i, j]] = relax dx/ap[[i, j]];residv += Abs[an[[i, j]]v[[i, j + 1]] + as[[i, j]]v[[i, j - 1]] + ae[[i, j]]v[[i + 1, j]] + aw[[i, j]]v[[i - 1, j]] - ap[[i, j]]v[[i, j]] + sc[[i, j]]], {i, 2, n - 1}, {j, 3, n - 1}]; ap /= relax; sc += (1 - relax)ap v; v = lisolv[2, 3, v];Do[an[[i, j]] = dx dv[[i, j + 1]]; as[[i, j]] = dx dv[[i, j]]; ae[[i, j]] = dx du[[i + 1, j]]; aw[[i, j]] = dx du[[i, j]]; sp[[i, j]] = 0; sc[[i, j]] = -((v[[i, j + 1]] - v[[i, j]])dx + (u[[i + 1, j]] - u[[i, j]]) dx); residp += Abs[sc[[i, j]]]; ap[[i, j]] = an[[i, j]] + as[[i, j]] + ae[[i, j]] + aw[[i, j]] - sp[[i, j]]; pstar[[i, j]] = 0, {i, 2, n - 1}, {j, 2, n - 1}]; Do[pstar = lisolv[2, 2, pstar], {2}]; Do[If[i != 2, u[[i, j]] += du[[i, j]](pstar[[i - 1, j]] - pstar[[i, j]])]; If[j != 2, v[[i, j]] += dv[[i, j]](pstar[[i, j - 1]] - pstar[[i, j]])]; p[[i, j]] += 0.3(pstar[[i, j]] - pstar[[n - 1, n - 1]]), {i, 2, n - 1}, {j, 2, n - 1}]]; |
 Here is some Mathematica code to plot streamlines. The blue lines are rotating clockwise and the red lines are rotating counter-clockwise. The plot on the left agrees fairly well with Ercan Erturk’s results. psi = Table[0, {n - 1}, {n - 1}]; Do[psi[[i, j]] = psi[[i, j - 1]] + dx(u[[i + 1, j - 1]] + u[[i + 1, j]])/2, {j, 2, n - 1}, {i, 1, n - 1}]; psi = ListInterpolation[psi, {{0, 1}, {0, 1}}]; ContourPlot[psi[x, y], {x, 0, 1}, {y, 0, 1}, PlotPoints -> 50, PlotRange -> All, ContourShading -> False,Contours -> {-0.08, -0.077, -0.07, -0.06, -0.045, -0.025, -0.01, -0.0025, 0, -8*^-6, 2*^-6, 3*^-5, 8*^-5, 1*^-4, 3*^-4, 6*^-4, 9*^-4}, ContourStyle -> Table[{Hue[2(1 - x)/3]}, {x, 0, 1, 1/16}]] |
 Here is some Mathematica code to plot vorticity contours. The plot on the left agrees fairly well with Ghia’s results. omega = ListInterpolation[Table[(v[[i, j]] - v[[i - 1, j]])/dx - (u[[i, j]] - u[[i, j - 1]])/dx, {i, 2, n - 1}, {j, 2, n - 1}], {{0, 1}, {0, 1}}]; ContourPlot[omega[x, y], {x, 0, 1}, {y, 0, 1},PlotPoints -> 50, PlotRange -> All, ContourShading -> False, Contours -> Range[-14, 7],ContourStyle -> Table[{Hue[2(1 - x)/3]}, {x, 0, 1, 1/21}]] Computational Fluid Dynamics (CFD) Links CFD Gallery - I especially like the simulations of turbulence Turbulent Boundary Layer - large simulation by Said Elghobashi More CFD Movies - nice simulations of fire Colorful Fluid Mixing Gallery Gallery of Fluid Dynamics - excellent web site on fluid dynamics by Mark Cramer Gallery of Fluid Motion - from the American Institute of Physics CFD code - simple Fortran codes by Abdusamad Salih FlowLab - educational resource by Fluent |
Driven Cavity - Mathematica 4.2, 11/16/06
|
Here is the same driven cavity using the finite difference method. This code uses the Alternating Direction Implicit (ADI) method to solve the vorticity transport equation assuming unsteady, incompressible, viscous flow. (* runtime: 1 minute *) << LinearAlgebra`Tridiagonal`; n = 20; RE = 1000; dx = 1.0/(n - 1); dt = 0.25; cxy = dt/(4dx); sxy = dt/(2RE dx^2); lambda = 1.5; free = Range[2, n - 1]; omega = v = psi = Table[0, {n}, {n}]; u = Table[If[j == n, 1, 0], {n}, {j, 1, n}]; Do[omega[[All, 1]] = 2psi[[All, 2]]/dx^2; omega[[All, n]] = 2 (psi[[All, n - 1]] + dx)/dx^2; omega[[1, All]] = 2 psi[[2, All]]/dx^2; omega[[n, All]] = 2 psi[[n - 1, All]]/dx^2; omega[[free, free]] = Transpose[Partition[TridiagonalSolve[Drop[Flatten[Table[If[i == n - 1, 0, -(cxy u[[i, j]] + sxy)], {j, 2, n - 1}, {i, 2, n - 1}]], -1], Flatten[Table[1 + 2sxy, {j, 2, n - 1}, {i, 2, n - 1}]], Drop[Flatten[Table[If[i == n - 1,0, cxy u[[i + 1, j]] - sxy], {j, 2,n - 1}, {i, 2, n - 1}]], -1], Flatten[Table[(cxy v[[i, j - 1]] + sxy)omega[[i, j - 1]] + (1 - 2sxy)omega[[i, j]] + (-cxy v[[i, j + 1]] + sxy)omega[[i, j + 1]] + If[i == 2, (cxy u[[i - 1, j]] + sxy)omega[[i - 1, j]],0] - If[i == n - 1, (cxy u[[i + 1, j]] - sxy) omega[[i + 1, j]], 0], {j, 2, n - 1}, {i, 2, n - 1}]]], n - 2]]; omega[[free, free]] = Partition[TridiagonalSolve[Drop[Flatten[Table[If[j == n - 1, 0, -(cxy v[[i, j]] + sxy)], {i, 2, n - 1}, {j, 2, n - 1}]], -1], Flatten[Table[1 + 2sxy, {i, 2, n - 1}, {j, 2, n - 1}]], Drop[Flatten[Table[If[j == n - 1, 0, cxy v[[i, j + 1]] - sxy], {i, 2, n - 1}, {j, 2, n - 1}]], -1], Flatten[Table[(cxy u[[i - 1, j]] + sxy)omega[[i - 1, j]] + (1 - 2sxy)omega[[i, j]] + (-cxy u[[i + 1, j]] +sxy)omega[[i + 1, j]] + If[j == 2, (cxy v[[i, j - 1]] + sxy)omega[[i, j - 1]], 0] - If[j == n - 1, (cxy v[[i, j + 1]] - sxy) omega[[i, j + 1]], 0], {i, 2, n - 1}, {j, 2, n - 1}]]], n - 2]; resid = 1; While[resid > 0.001, resid = 0; Do[resid += Abs[psi[[i, j]] - (psi[[i,j]] = lambda(psi[[i, j - 1]] + psi[[i - 1, j]] + psi[[i + 1, j]] + psi[[i, j + 1]] - dx^2omega[[i, j]])/4 + (1 - lambda)psi[[i, j]])], {i, 2, n - 1}, {j, 2, n - 1}]]; Do[u[[i, j]] = (psi[[i, j + 1]] - psi[[i, j - 1]])/(2dx); v[[i, j]] = -(psi[[i + 1, j]] - psi[[i - 1, j]])/(2dx), {i, 2, n - 1}, {j, 2, n - 1}]; ListContourPlot[Transpose[psi], PlotRange -> All, ContourShading -> False], {100}]; Link: Driven Cavity - simple Fortran program by Zheming Zheng that uses this method |
Raindrop Light Dispersion - Mathematica 4.2, 1/21/05; C++ version: 7/26/05
When light enters near the side of a raindrop, much of the light is internally reflected. This is how rainbows are formed. You can also see this animation on Jeff Bryant’s Mathematica visualization site.(* runtime: 12 seconds *) Spectrum = {{380, {0, 0, 0}}, {420, {1/3, 0, 1}}, {440, {0, 0, 1}}, {490, {0, 1, 1}}, {510, {0, 1,0}}, {580, {1, 1, 0}}, {645, {1, 0, 0}}, {700, {1, 0, 0}}, {780, {0, 0, 0}}}; Snell[theta_, n1_, n2_] := ArcSin[Min[1, Max[-1, Sin[theta]n1/n2]]]; FresnelR[t1_, t2_, n1_, n2_] := (((n1 Cos[t2] - n2 Cos[t1])/(n1 Cos[t2] + n2 Cos[t1]))^2 + ((n1 Cos[t1] - n2 Cos[t2])/(n1 Cos[t1] + n2 Cos[t2]))^2)/2; Angle[{x1_, y1_}, {x2_, y2_}] := ArcTan[x2 - x1, y2 - y1]; IntersCircle[{x_, y_}, theta_, {x0_,y0_}, r_] := Module[{t1 = Tan[theta], t2, a, b, x2}, t2 = ((x -x0) t1 - (y - y0))/r; a = t1^2 + 1; b = -2t1 t2; x2 = x0 + r(-b + If[Cos[theta] > 0, 1, -1]Sqrt[b^2 - 4a (t2^2 - 1)])/(2a); {x2, y + (x2 - x)t1}]; AddLine[{x1_, y1_}, {x2_, y2_}, color_] := Do[image[[Round[m(s y1 + (1 - s) y2)], Round[m(s x1 + (1 - s)x2)]]] += color, {s, 0,1, 1/(m Sqrt[(x2 - x1)^2 + (y2 - y1)^2])}]; m = 275; image =Table[{0, 0, 0}, {m}, {m}]; x0 = y0 = 0.5; p0 = {x0, y0}; r = 0.25; Do[n = 1.31Exp[12.4/lambda]; color = Table[Interpolation[Map[{#[[1]], #[[2, i]]} &, Spectrum], InterpolationOrder -> 1][lambda], {i, 1, 3}]; p1 = {x0 - Sqrt[r^2 - (y - y0)^2], y}; t = Angle[p0, p1]; t1 = t + Snell[Pi - t, 1,n] - Pi; While[Max[color] > 0.01, p2 = IntersCircle[p1, t1, p0, r]; AddLine[p1,p2, color]; t = Angle[p0, p2]; t2 =2t - t1 - Pi; t3 = t - Snell[t - t1, n, 1]; R1 = FresnelR[t - t1, t -t3, n, 1]; AddLine[p2, IntersCircle[p2, t3, p0, 0.49], (1 - R1) color]; p1 = p2; t1 = t2; color *= R1], {lambda, 380, 780, 25}, {y, y0 + r - 0.001, y0 + r - 0.0001, 0.0001}]; Show[Graphics[RasterArray[Map[RGBColor @@ Map[Min[1, #] &, #] &, image, {2}]], Epilog -> {{RGBColor[1, 1, 1], Circle[m p0, m r]}}, AspectRatio -> 1]] Click here to see some POV-Ray code for a similar-looking image using photon mapping. |
Diamond Light Dispersion - Mathematica 4.2, 1/29/05
Diamonds have a high index of refration (n = 2.4195) and coefficient of dispersion (COD = 0.044). When light reaches a boundary, it can refract or reflect, in proportion to the Fresnel Reflection Coefficient. The angle of the refracted light is given by Snell’s Law: n1sin(θ1) = n2sin(θ2)Links Click here to see my rock and gemstone collection. Ray Traced Diamond - POV-Ray code by PM 2Ring, based on Marcel Tolkowsky’s “ideal” diamond Prism Refraction - Java applet 3D Dispersion - C++ renderings by Werner Benger |
Galactic Gravity Simulations - C++, 11/13/07
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 Milky Way and Andromeda Galaxy Collision - amazing simulations by John Dubinski Millennium Simulation - interesting video, see also Dark Matter Simulation Type Ia Supernova - simulations Simulation of Newtonian Gravity - Java applet by Mathew Tizard Hyperion Spiral Galaxy Demo - beautiful simulation by Philippe Guglielmetti Colliding Galaxies - fun Java applet, originally by Uli Siegmund galaxy - nice-looking Java applet Universe Holocube - 100 Megaparsec scale laser-etched crystal cube Milky Way Galaxy Holocube - scale laser-etched crystal cube Astronomy Picture of the Day - amazing astronomy pictures, some of my favorites are Crab Pulsar, Cat's Eye Nebula, Eta Carinae, Variable star V838 Monocerotis, and Hourglass Nebula Cosmographica - impressive astronomical art gallery by Don Dixon Breakthrough Propulsion Physics (BPP) - is warp drive possible? |
Von Kármán Vortex Street - C++, 1/11/05
 When fluid passes an object, it can leave a trail of vortices called a Von Kármán Vortex Street. This animation shows the vorticity where blue is clockwise and red is counter-clockwise. This simulation assumes unsteady, incompressible, viscid, laminar flow. For solenoidal flows, mass conservation can be achieved by taking the Fast Fourier Transform (FFT) of the velocity and then removing the radial component of the wave number vectors. The code was adapted from Jos Stam’s C Source Code and paper which is patented by Alias. (* Note: Something is wrong with this code. Please tell me if you find out what it is. runtime: 14 seconds *) n = 64; dt = 0.3; mu = 0.001; v = Table[{0, 0}, {n}, {n}]; Do[Do[If[i < 5, v[[i, j]] = {0.1, 0}]; If[(i - n/4)^2 + (j -n/2)^2 < 4^2, v[[i, j]] = {0, 0}], {i, 1, n}, {j, 1, n}]; ui = ListInterpolation[v[[All, All,1]]]; vi = ListInterpolation[v[[All, All, 2]]]; v = Table[{i2, j2} = {i, j} - n dt v[[i, j]]; {ui[i2, j2], vi[i2, j2]}, {i, 1, n}, {j, 1, n}]; v = Transpose[Map[Fourier[v[[All, All, #]]] &, {1,2}], {3, 1, 2}]; v = Table[x = Mod[i + n/2, n] - n/2; y = Mod[j + n/2, n] - n/2; k = x^2 + y^2; Exp[-k dt mu]If[k > 0, (v[[i, j]].{-y, x}/k){-y, x}, v[[i, j]]], {i, 1, n}, {j, 1, n}]; v = Transpose[Map[Re[InverseFourier[v[[All, All, #]]]] &, {1, 2}], {3, 1, 2}]; ListDensityPlot[Table[((v[[i + 1, j, 2]] - v[[i - 1, j, 2]]) - (v[[i, j + 1, 1]] - v[[i, j - 1, 1]]))/2, {j, 2, n - 1}, {i, 2, n - 1}], Mesh -> False,Frame -> False, ColorFunction -> (Hue[2#/3] &)], {t, 1, 25}]; Links smoke program - uses this method, by Gustav Taxén fire program - uses this method, Hugo Elias Fluid Simulations - nice vortex street simulations by Maciej Matyka Vortex street in the clouds - photo of a vortex street created by an island Airplane Vortices - nice photo of trailing airplane vortices Harmonic Balance CFD |
Kelvin-Helmholtz Instability Waves - Mathematica 4.2, C++, POV-Ray 3.6.1, 11/16/06
Kelvin-Helmholtz instabilities are formed when one fluid layer is on top of another fluid layer moving with a different velocity. These instabilities take the form of small waves that can eventually grow into vortex rollers. This is a purely potential flow (irrotational) phenomenon. The following Mathematica code was adapted from Zheming Zheng’s Fortran program. It uses vortex blobs to simulate smooth rollers and periodic point vortices to simulate the far field. A simple predictor-corrector scheme is used to integrate the differential equations:(* runtime: 9 minutes *) n = 40; dt = 0.05; L = 2Pi; gamma = L/n; rcore = 0.5; plist = Table[{x, 0.01 Sin[2Pi x/L]}, {x, 0, L(1 - 1/n), L/n}]; vcalc[plist_] := Table[Sum[{x, y} = plist[[j]] - plist[[i]]; gamma/(2Pi) Sum[If[i == j && k == 0, 0, r2 = (x + k L)^2 +y^2; {-y, x + k L}/r2(1 - Exp[-r2/rcore^2] - 1)], {k, -3, 3}] + If[i ==j, 0, gamma/(2L) {-Sinh[2Pi y/L], Sin[2Pi x/L]}/(Cosh[2Pi y/L] - Cos[2Pi x/L])], {j, 1, n}], {i, 1, n}]; Do[vlist = vcalc[plist]; vlist0 =vlist; vlist = vcalc[plist + dt vlist]; plist += dt(vlist0 + vlist)/2; ListPlot[plist, PlotJoined -> True, PlotRange -> {-2, 2}, AspectRatio -> Automatic], {i, 1, 400}]; |
Here is some Mathematica code to plot streamlines assuming periodic point vortices:(* runtime: 2 seconds *) Clear[psi]; psi[x_, y_] := Plus @@ Map[-gamma/(4 Pi) Re[Log[Cos[2Pi(x - #[[1]])/L] - Cosh[2Pi (y - #[[2]])/L]]] &, plist]; ContourPlot[psi[x, y], {x, 0, L}, {y, -L/3, L/3}, PlotRange -> All, ContourShading -> False, PlotPoints -> 20, AspectRatio -> Automatic] Link: Kelvin-Helmholtz waves in the clouds |
Here is another variation showing vortex pairing and pathlines.Link: movie - famous experiment by Bernal |
Rayleigh Surface Waves - adapted from Daniel Russell's Mathematica code, C++, 11/2/07
Rayleigh Surface Waves are associated with ocean waves and earthquakes. The following Mathematica code was adapted from Daniel Russell's Rayleigh Surface Waves notebook:(* runtime: 0.2 second *) omega = 2 Pi; k = 7; nu = 0.33; n = 0.175; xi = Sqrt[(0.5 - nu)/(1 - nu)]; q = k Sqrt[1 - n^2 xi^2]; s = k Sqrt[1 - n^2]; Do[grid = Table[{x + k (Exp[q y] - 2 q s Exp[s y]/(s^2 + k^2)) Sin[omega t - k x], y + q (Exp[q y] - 2 k^2 Exp[s y]/(s^2 + k^2)) Cos[omega t - k x]}, {x, 0, 2, 0.1}, {y, -1, 0, 0.1}]; Show[Graphics[{Map[Line, grid], Map[Line, Transpose[grid]]}, AspectRatio -> Automatic, PlotRange -> {{-0.1, 2.1}, {-1.1, 0.25}}]], {t, 0, 1, 0.05}]; Links Animated Gif - by Daniel Russell Oceanography Waves - webpage by J Floor Anthoni |
Magnus Effect - C++, Mathematica 4.2, POV-Ray 3.6.1, 10/24/07
A Curveball is a spinning ball that accelerates sideways due to the Magnus Effect. Here is some Mathematica code to plot the streamlines around a lifting cylinder, assuming inviscid incompressible potential flow (irrotational):(* runtime: 0.3 second *) U = 1; gamma = 4; a = 1; w[z_] := U(z + a^2/z) + I gamma Log[z]/(2Pi); ContourPlot[Im[w[x + I y]], {x, -2, 2}, {y, -2, 2}, Contours -> Table[psi, {psi, -3, 3, 0.25}], PlotPoints -> 100, ContourShading -> False, AspectRatio -> Automatic] Here is some Mathematica code to plot the entrained fluid using the 4th order Runge-Kutta method: (* runtime: 17 seconds, increase n for better resolution *) Clear[v]; n = 40; U = 1; gamma = 4; a = 1; tmax = 5.0; dt = tmax/100; v[z_] := If[Abs[z] < 1, 0, Conjugate[U (1 - (a/z)^2) + I gamma/(2Pi z)]]; image = Table[z = x + I y; Do[k1 = dt v[z]; k2 = dt v[z + k1/2]; k3 = dt v[z + k2/2]; k4 =dt v[z + k3]; z -= (k1 + 2 k2 + 2 k3 + k4)/6, {t, tmax, 0, -dt}]; z, {y, -2, 2, 4.0/n}, {x, -2, 2, 4.0/n}]; ListDensityPlot[Map[Abs[Mod[#, 1]] &, image, {2}], Mesh -> False, Frame -> False, AspectRatio -> 1] Link: My Educational Project - some basic potential flows using Matlab and Mathematica |
Horseshoe Magnets - new version: POV-Ray 3.6.1, 7/15/06
old version: calculated in Mathematica 4.2, rendered in TrueSpace 4.3, 3/9/03
 These magnetic fields were made using the same technique as the solenoid picture. |
LinksLevitron - magnetic levitating top Molecular Magnetism - magnetic levitation of a frog in a Bitter solenoid Gravitator - magnetic levitating sphere Geodynamo Lab - I recently had the opportunity to tour this lab at the University of Maryland Electromagnetic Field in Microwave Oven - interesting discussion by Wenjun Deng (my roommate) |
Halo Orbit - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 12/1/05
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 Interplanetary Superhighway - by Shane Ross Stable 3-Body Orbits - Java applet by Robert Vanderbei, be sure to see the Ducati 3-Body Orbit |
Gravitational Lensing of a Black Hole - Mathematica 4.2, 2/5/05
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 < 0.1^2, 0, stars[Mod[0.5((1 - 0.25/r) {x, y} + 1), 1]]]], {x, -1, 1}, {y, -1, 1}, PlotPoints ->275, PlotRange -> {0, 1}, Mesh -> False, Frame -> False] Links "Monster of the Milky Way" - Nova video about black holes Black Holes - BBC Video with Sam Neill ART project - relativity simulations, see the movies “Tackling the Riddles of Gravity” and “Numerical Relativity” Werner Benger - Black Hole Raytracing, Einstein Ring Animation, and simulation of Earth as a black hole Real-time Relativistic Rendering - by Hansong Zhang Virtual Trips to Black Holes - simulations by Robert Nemiroff Black Hole Simulation - by Phil Armitage Black Hole Flight Simulator - screenshots and Java applet, by Andrew Hamilton |
1D Shock Tube - Fortran 90, Mathematica 4.2, POV-Ray 3.6.1, 3/15/07
(* runtime: 5 seconds *) gamma = 1.4; R[W_] := Module[{}, rho = W[[All, 1]]; u = W[[All, 2]]/rho; p = (gamma - 1)(W[[All, 3]] - rho u^2/2); F = u W + Transpose[{Table[0, {n}], p, u p}]; h = Table[(F[[Min[n, i + 1]]] + F[[i]])/2, {i,1, n}]; Q = Table[h[[Max[i, 2]]] - h[[Max[i, 2] - 1]], {i, 1, n}]; nu = Table[i = Max[2, Min[n - 1, i]]; Abs[(p[[i + 1]] - 2p[[i]] + p[[i - 1]])/(p[[i + 1]] + 2p[[i]] + p[[i - 1]])], {i, 1, n}]; S = Table[Max[nu[[Min[n, i + 1]]], nu[[i]]], {i, 1, n}]; alpha1 = 1/2; beta1 = 1/4;alpha2 = alpha1; beta2 = beta1; epsilon2 = Map[Min[alpha1, alpha2#] &, S];epsilon4 = Map[Max[0, beta1 - beta2#] &, epsilon2];dW = Table[W[[Min[n - 1, i] + 1]] - W[[Min[n - 1, i]]], {i, 1, n}];dW3 = Table[i = Max[2, Min[n - 2, i]]; -W[[i - 1]] + 3W[[i]] - 3W[[i + 1]] + W[[i + 2]], {i, 1, n}];d = (epsilon2 dW - epsilon4 dW3)(Abs[u] + a); Dflux = Table[d[[Max[2, i]]] - d[[Max[2, i] - 1]], {i, 1, n}]; (Q - Dflux)/dx]; n = 50; dx = 1.0/n; a = 1.0; dt = dx/(1.0 + a); W = Table[{If[i > n/2, 0.125, 1], 0, If[i > n/2, 0.1, 1]/(gamma - 1)}, {i, 1, n}]; Do[W -= dt R[W - dt R[W - dt R[W]/4]/3]/2;ListPlot[W[[All, 1]], PlotJoined -> True,PlotRange -> {0, 1}, AxesLabel -> {"i", "rho"}], {t, 0, 100dt, dt}]; |
Vibration Mode of a Spherical Membrane - Mathematica 4.2, MathGL3d, 9/28/04
This is a somewhat exaggerated example of how a spherical balloon might vibrate.(* runtime: 4 seconds *) Y := Re[SphericalHarmonicY[8, 4, theta, phi]]; r := (1 + 0.5Y); ParametricPlot3D[{r Sin[theta] Cos[phi], r Sin[theta] Sin[phi], r Cos[theta], {EdgeForm[], SurfaceColor[Hue[Y]]}}, {theta, 0, Pi}, {phi, 0, 2Pi}, PlotPoints -> 72, Boxed -> False, Axes -> None] Links Click here to read more about spherical harmonics. Spherical Harmonic Lighting - interesting technique for global illumination by Paul Baker |
Hydrogen Electron Orbital Probability Distribution Cross-Sections - Mathematica 4.2, 4/30/03
Here are some different shapes that a hydrogen atom can take. These plots are based on the same function as the vibrating balloon (shown above):
P = 4 π r2ψr2ψθ2, ψr = e-ρρlLn-1-l2l+1(ρ), ψθ = Ylm(θ,φ), ρ = 2r/n
|
(* runtime: ggg second *) P[n_, l_, m_, {x_, y_, z_}] := Module[{r = Sqrt[x^2 + y^2 + z^2]}, 4 Pi r^2(Exp[-r/n]r^l LaguerreL[n - 1 - l, 2l + 1, 2 r/n])^2 Abs[SphericalHarmonicY[l, m, ArcCos[z/r], ArcTan[x, y]]]^2]; DensityPlot[P[3, 1, 0, {x, 0, z}], {x, -20, 20}, {z, -20, 20}, Mesh -> False, PlotPoints -> 275, Frame -> False, Compiled -> False] |
Hydrogen 4f Electron Orbital Probability Distribution - Mathematica 4.2, 8/22/04
animated version: C++, POV-Ray 3.6.1, 11/27/05
Here is a 3D view of a hydrogren atom in the 4f state. The “gas-like” cloud is composed of a volume of 2D cross-sectional layers. The animations were made in POV-Ray using DF3 density files. The right animation shows what a hypothetical 12o orbital might look like.POV-Ray has a built-in internal function for the 3d orbital: // runtime: 4 seconds camera{location 16*z look_at 0} #declare P=function{internal(53)}; #declare P0=P(0,3,0,0); box{-8,8 pigment{rgbt t} hollow interior{media{emission 0.5 density{function{(P(x,y,z,0)-1.2)/(P0-1.2)} color_map{[0 rgb 0][1 rgb 1]}}}}} Links HydrogenLab - animations of hydrogen energy transitions |
Jet - C++, 6/6/07
This simulation is based on a simple method presented in this paper and sample C++ code by Jos Stam.(* runtime: * seconds *) |
Zoom Animation - Google Maps API, POV-Ray 3.6.1, 10/29/07
I made this using Google Maps.Links Google Sightseeing - discover the world via Google Maps Celestia - impressive free 3D space simulation software by Chris Laurel Celestia Addon Catalogs - Joseph Bolte, Bruckner’s Celestia Page JTrack 3D - fun satellite tracker Java applet Powers of 10 - Java applet Google Earth - fun program to explore satellite imagery around the globe Earth POV-Ray Code - by James Kuffner Google Maps API Developer's Guide |
Schwarz-Christoffel Flow - C++, 11/*/07
Here is a potential flow field of a vortex near a flat wall mapped to stair steps using a Schwarz-Christoffel transformation:(* runtime: 0.7 second *) Clear[z]; gamma = 1; z0 = 1 + I; z1 = I;z2 = 0; z3 = 1; zeta1 = -1; zeta2 = 0; zeta3 = 1; theta1 = 3Pi/2; theta2 = Pi/2; theta3 = 3Pi/2; z = z[zeta] /. DSolve[{z'[zeta] == K(zeta - zeta1)^(theta1/Pi - 1)(zeta - zeta2)^(theta2/Pi - 1)(zeta -zeta3)^(theta3/Pi - 1), z[zeta1] == z1}, z[zeta], zeta][[1]]; z = z /. K -> (K /. Solve[(z /. zeta -> zeta2) == z2, K][[1]]); zeta0 = zeta /. FindRoot[z == z0, {zeta, zeta2 + I}]; Z[w_] = z /. zeta -> (zeta /. Solve[w == (I gamma/(2Pi))(Log[zeta - zeta0] - Log[zeta - Conjugate[zeta0]]), zeta][[1]]); ToPoint[z_] := {Re[z], Im[z]}; grid = Table[ToPoint[Z[phi + I psi]], {phi, -1, 0, 0.05}, {psi, -1*^-6, -1, -0.02}]; Show[Graphics[{Map[Line, grid], Map[Line, Transpose[grid]]}, AspectRatio -> Automatic, PlotRange -> {{-1, 2}, {-1, 2}}]] Links My Educational Project - some basic potential flows using Matlab and Mathematica Schwarz-Christoffel Maps - nice plots by Matthias Weber used for minimal surfaces Schwarz-Christoffel Transformations - many examples by John Mathews |
Periodic Steady Vortices in a Stagnation Point Flow - C++, 11/17/07
I got the idea for this potential flow animation from this paper. The entrained fluid was integrated using the 4th order Runge-Kutta method. |
Fundamental Modes of Vibration for a Violin-Shaped Membrane - GeoStar 2.8, AutoCAD 2000, 4/21/04
Mode displacement plots of the 19th, 28th, and 45th modes, respectively, from left to right. Click here to see a holographic interferogram of a real violin.
|
Mode displacement plot of the 18th mode of a Sunset Moth-Shaped Membrane. Just for fun.
|
31st Fundamental Mode of Vibration of a Circular Membrane - Mathematica 4.2, MathGL3d, 11/10/03
This is an example of how a circular drumhead might vibrate. Click here to see a large animation. Click here to see a rotatable animated version.J3(k43r)cos(3θ) → frequency: f31 = f1k43/k10 = 6.74621 f1 (* runtime: 5 seconds *) Clear[r]; << NumericalMath`BesselZeros`; k = BesselJZeros[3, 4][[4]]; << MathGL3d`OpenGLViewer`; texture = Graphics[{RGBColor[0, 1, 1], Rectangle[{0, 0}, {10, 10}], RGBColor[0.5, 0, 1], Rectangle[{1, 1}, {9, 9}]}]; MVShow3D[ParametricPlot3D[{r Sin[theta], r Cos[theta], 0.25BesselJ[3, k r]Cos[3 theta]}, {r, 0, 1}, {theta, 0, 2Pi}, PlotPoints -> {56, 168}], MVNewScene -> True, MVTexture -> texture, MVTextureMapType -> MVMeshUVMapping, MVScaleTexture -> {14, 42}] Links Vibrating Drumheads - by Jim Swift Cymatics - wave patterns created by sound waves |
Magnifying Glass - POV-Ray 3.6.1, 3/20/07
This magnifying glass demonstrates refraction and dispersion. The magnifying glass was created by intersecting two paraboloids. |
Double Pendulum - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 1/9/06
This is a classical example of chaotic motion.(* runtime: 0.1 second *) Clear[theta1, theta2]; g = 9.81; soln = NDSolve[{theta1''[t](3 - Cos[2(theta1[t] - theta2[t])]) == -3g Sin[theta1[t]] - g Sin[theta1[t] - 2theta2[t]] - 2Sin[theta1[t] - theta2[t]](theta2'[t]^2 - theta1'[t]^2Cos[theta1[t] - theta2[t]]), theta2''[t](3 - Cos[2(theta1[t] - theta2[t])]) == 2Sin[theta1[t] - theta2[t]](2theta1'[t]^2 + 2g Cos[theta1[t]] + theta2'[t]^2Cos[theta1[t] - theta2[t]]), theta1[0] == 0.999Pi, theta2[0] == 0.999Pi, theta1'[0] == 0, theta2'[0] == 0}, {theta1[t], theta2[t]}, {t,0, 10}][[1]]; theta1[t_] = theta1[t] /. soln; theta2[t_] = theta2[t] /. soln; Do[p1 = -{Sin[theta1[t]], Cos[theta1[t]]}; p2 = p1 - {Sin[theta2[t]], Cos[theta2[t]]}; Show[Graphics[{Line[{{0, 0}, p1, p2}], Disk[p1, 0.1], Disk[p2, 0.1]}], Axes -> True, PlotRange -> 2{{-1, 1}, {-1, 1}}, AspectRatio -> Automatic], {t, 0,10, 0.1}]; Link: Double Pendulum - Java applet by Erik Neumann |
Spherical Pendulum - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 1/10/06
The balls in this frictionless bowl move like spherical pendulums. See also my magnetic pendulum fractals.(* runtime: 0.03 second *) Clear[theta, phi]; g = 9.81; soln = NDSolve[{theta''[t] == (phi'[t]^2Cos[theta[t]] - g)Sin[theta[t]], phi''[t]Tan[theta[t]] == -2phi'[t]theta'[t], theta[0] == Pi/2, phi[0] == 0, theta'[0] == 0, phi'[0] == 2.357}, {theta[t], phi[t]}, {t, 0, 8}][[1]]; theta[t_] = theta[t] /. soln; phi[t_] = phi[t] /. soln; ParametricPlot3D[{Sin[theta[t]]Cos[phi[t]], Sin[theta[t]]Sin[phi[t]], -Cos[theta[t]]}, {t, 0, 8}] Here’s some Mathematica code to animate it: (* runtime: 0.6 second *) << Graphics`Shapes`; Do[p = {Sin[theta[t]]Cos[phi[t]], Sin[theta[t]]Sin[phi[t]], -Cos[theta[t]]}; Show[Graphics3D[{EdgeForm[], Line[{{0, 0, 0}, p}], TranslateShape[Sphere[0.2, 10, 6], p]}, PlotRange -> 1.2{{-1, 1}, {-1, 1}, {-1, 1}}]], {t, 0, 8, 0.1}]; |
Cloth Simulation - Mathematica 4.2, MathGL3d, 7/11/05
This cloth is modelled as a net of small springs and masses. The following code still needs some improvement:(* runtime: 2 minutes *) Normalize[p_] := p/Sqrt[p.p]; r = 0.5; n = 13; dx = 2.0/(n - 1); dt = 0.05; g = {0, 0, -0.2}; cloth = Table[{{x, y, r}, {0, 0, 0}}, {y, -1, 1, dx}, {x, -1, 1, dx}]; Do[ Show[Graphics3D[{EdgeForm[], Map[Polygon[Flatten[#, 1][[{1,2, 4, 3}, 1]]] &, Partition[cloth, {2, 2}, 1], {2}]}, PlotRange -> {{-1,1}, {-1, 1}, {-1, 1}}]]; Do[cloth = Table[{p1, v1} = cloth[[i, j]]; If[(i == 1 || i == n) && (j == 1 || j == n), {p1, v1},v2 = v1 + g dt; Do[i2 = i + di; j2 = j + dj; If[! (i2 == i && j2 == j) && 0 < i2 <= n && 0 < j2 <= n, L0 = dx Sqrt[(j2 - j)^2 + (i2 - i)^2]; p2 = cloth[[i2, j2, 1]]; L = Sqrt[(p2 - p1).(p2 - p1)]; v2 += (L/L0 - 1)Normalize[p2 - p1]dt], {di, -2, 2}, {dj, -2, 2}]; v2 *= 0.9; {p1 + (v1 + v2) dt/2,v2}], {i, 1, n}, {j, 1, n}], {25}], {10}]; Links OpenGL program - by Paul Baker Cloth Animations - by Ron Fedkiw, Robert Bridson, et al. Cloth and Fur Energy Functions - Pixar presentation |
Fraunhofer Diffraction of Light Through a Circular Aperature - Mathematica 4.2, 8/24/04
When parallel light passes through a small hole and a shadow is cast on a wall, there will be small diffraction ridges along the edge of the shadow. These ridges become smaller when the frequency of the light is increased. See also my single slit diffraction of a green laser.(* runtime: 16 seconds *) Clear[i]; lambda = 500*^-9; k = 2Pi/lambda; a = 80*^-6; i[rho_] := (2Pi a BesselJ[1, k a rho]/(k rho))^2; DensityPlot[i[Sqrt[u^2 + v^2]], {u, -0.01, 0.01}, {v, -0.01, 0.01}, Mesh -> False, PlotPoints -> 275, Frame -> False] |
Vortices - Mathematica 4.2, C++, 10/31/06
Here is another code to solve Euler’s equations (inviscid flow) and plot the pathlines. This method was adapted from Stephen Montgomery-Smith’s Euler2D program. The vortex effects are found using the Biot-Savart law and the differential equations are solved using the Adams Bashforth method. Click here to download some Matlab code for this. Shown below is some Mathematica code:(* runtime: 19 seconds *) ToPoint[z_] := {Re[z], Im[z]};dt = 0.01; rcore = 0.1; dz0 = 0; Klist = {1, -1, 1, -1}; zlist = Join[{-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}, Flatten[Table[x + I y, {x, -1, 1, 0.1}, {y, -1, 1, 0.1}], 1]]; paths = Map[# Table[1, {10}] &, zlist]; Do[dz = Map[Sum[zj = zlist[[j]]; If[# != zj, r2 = Abs[# - zj]^2; I Klist[[j]](zj - #)/r2 (1 -Exp[-r2/rcore^2]), 0], {j, 1, Length[Klist]}] &, zlist]; zlist += dt (1.5dz - 0.5dz0); dz0 = dz; paths = Table[Prepend[Delete[paths[[i]], -1], zlist[[i]]], {i, 1, Length[zlist]}];Show[Graphics[Map[Line[Map[ToPoint, #]] &, paths], PlotRange -> {{-1, 1}, {-1, 1}}, AspectRatio -> Automatic]], {85}]; Link: math paper by Stephen Montgomery-Smith, seems fairly advanced to me |
Channel Flow - Fluent 5.5.14, Gambit, Mathematica 4.2, 5/14/05
Vibrating String with Damping - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 1/27/06
 This analytical solution for a vibrating string was solved using the method of separation of variables assuming small linear perturbations. The rope was drawn using Mike Williams’ rope macro. (* runtime: 2 seconds *) c = 4; L = 1; h = 0.1; b = Pi; k := (n - 0.5)Pi/L; omega := Sqrt[(c k)^2 - (0.5b)^2]; Do[Plot[Evaluate[Sum[(6h/k^2)Sin[k/3] Exp[-0.5b t](Cos[omega t] + (0.5b/omega) Sin[omega t])Sin[k x], {n, 1, 10}]], {x, 0,L}, PlotRange -> {{0, L}, {-h, h}}, AspectRatio -> Automatic], {t, 0, 1.5, 0.025}]; |
Wave Interference - Mathematica 4.2, 8/24/04
This interference pattern is the superposition of two wave sources, such as light passing through two narrow slits.(* runtime: 9 seconds *) Do[DensityPlot[Sin[2Pi(Sqrt[x^2 + (y - 3)^2] + t)] + Sin[2Pi(Sqrt[x^2 + (y + 3)^2] + t)], {x, 0, 16}, {y, -8, 8}, PlotPoints -> 275, Mesh -> False, Frame -> False], {t, 0, 1, 0.1}]; |
Physics Links
The Elegant Universe - interesting 3 hour online mini-series on controversial string theory by NOVA
Hydraulophone video - a hydraulophone works just like a flute except it uses water instead of air, the world's largest hydraulophone is at the Ontario Science Centre in Canada
Physics Lectures - free videos by Walter Lewin at MIT
Nanoscale Demonstrations - interesting physics movies
fun and simple homemade devices - by David Williamson
Physics Demonstrations - University of Maryland
HyperPhysics - Georgia State University
|