## Posts Tagged ‘Mathematica

10
Feb
09

### Scherk-Collins Surface

This surface can be formed by twisting and warping a singly-periodic Scherk’s minimal surface. This idea was originally attributed to Brent Collins. Technically, the surface is no longer considered exactly “minimal” after twisting but it still looks minimal (it is actually very difficult to find the exact shape for most minimal surfaces). Click here to download some POV-Ray code.

Here is some Mathematica code:
(* runtime: 0.3 second *) << GraphicsMaster; n = 5; r = 0.75n; Twist[{x_, y_, z_}, theta_] := {x Cos[theta] - y Sin[theta], x Sin[theta] + y Cos[theta], z}; Warp[{x_, y_, z_}, theta_] := {(x + r) Cos[theta], (x + r) Sin[theta], y}; f[z_] := Module[{t1 = Sqrt[2Cot[z]], t2 = Cot[z] + 1}, Warp[Twist[Re[{0.5xsign(Log[t1 - t2] - Log[t1 + t2])/Sqrt[2], ysign I(ArcTan[1 - t1] - ArcTan[1 + t1])/Sqrt[2], z}], 2Re[z]/n], 2Re[z]/n]]; DisplayTogether[Table[ParametricPlot3D[f[x + I y], {x, 0, n Pi}, {y, 0.001, 0.75}, PlotPoints -> {8n + 1, 5}, Compiled -> False], {xsign, -1, 1, 2}, {ysign, -1, 1, 2}]]

The following Mathematica code can be used to increase the number of edges (or “branches”). This code uses some complicated functions that were adapted from Matthias Weber’s Mathematica notebook:
(* runtime: 1.2 seconds *) << GraphicsShapes; k = 4; phi = Pi(0.6/k - 0.5)/(1 - k); f[z_] := Re[NIntegrate[Evaluate[{0.5 (w^(1 - k) - w^(k - 1)), 0.5 I (w^(1 - k) + w^(k - 1)), 1}/(w^(k + 1) + w^(1 - k) - 2w Cos[k phi])], {w, 0, z}]]; alpha = Pi/k; zbeta = Exp[I Pi(phi/alpha - 0.5)]; surface = ParametricPlot3D[Re[f[Exp[I alpha/2]((1 + I zbeta Exp[r + I theta])/(I Exp[r + I theta] -zbeta))^(alpha/Pi)]], {r, 0, 4}, {theta, 0, Pi}, PlotPoints -> 10, Compiled -> False, DisplayFunction -> Identity][[1]]; z0 = f[1][[3]]; surface = {surface, AffineShape[TranslateShape[surface, {0, 0, -2z0}], {1, 1, -1}]}; surface = {surface, AffineShape[surface, {1, -1, 1}]}; surface = Table[RotateShape[surface, 2Pi i/k, 0, 0], {i, 1, k}]; dz = Pi Csc[k phi]/k; Show[Graphics3D[Table[TranslateShape[surface, {0, 0, i dz}], {i, 0, 1}]]] 

21
Jan
09

### Punctured Helicoid

Here is a helicoidwith holes in it. The following Mathematica code uses some complicated functions that were adapted from Matthias Weber’sMathematica notebook:
(* runtime: 4 seconds *) << GraphicsShapes; tau0 = Exp[1.23409 I]; b0 = 0.629065; theta[z_] := EllipticTheta[1, Pi z, Exp[I Pi tau0]]; r1[z_] := theta[z + 0.5 (b0 - 2) (tau0 + 1)]/theta[z + 0.5 (b0 - 1) (tau0 + 1)]; r2[z_] := theta[z - 0.5 b0 (tau0 + 1)]/theta[z - 0.5 (b0 + 1) (tau0 + 1)]; omega3[z_] := r1[z] r2[z]/(0.386191 - 0.169839 I); G[z_] := (108.37 - 62.8417 I) Exp[I Pi (b0 - 2 z + 2 tau0 + b0 tau0)]r1[z]/r2[z]; f[z0_] := Re[NIntegrate[Evaluate[{-(G[z] omega3[z] - omega3[z]/G[z] )/2, I(G[z] omega3[z] + omega3[z]/G[z] )/2, omega3[z]}], {z, tau0/2, z0}]] + {0.434156, 0, -1}; a0 = -0.409956; r0 = 2.43051; g[z_] := (EllipticF[ArcSin[(a0 + r0 E^z)/(1 - a0 E^z)], 1/r0^2]/(2EllipticF[Pi/2, 1/r0^2]) + 0.5)(1 + tau0)/2; surface = ParametricPlot3D[f[g[x + I y]], {x, -2.5, 2.5 - 0.8881}, {y, 0.001,0.999Pi}, PlotPoints -> {15, 10}, Compiled -> False, DisplayFunction -> Identity][[1]]; surface = {surface, RotateShape[surface, 0, 0, Pi]}; Show[Graphics3D[{surface, TranslateShape[surface, {0, 0, 2}]}, ViewPoint -> {1, 6, 3}]];

29
Dec
08

### Polynomial Roots

Strange fractal patterns emerge when you plot the complex roots of high order polynomials. This picture shows all the roots for all possible combinations of 18th order polynomials with coefficients of ±1. You can easily find the roots using Mathematica’s Root function:

(* runtime: 34 seconds *)
n = 12; m = 275; image = Table[0.0, {m}, {m}];
Do[Do[z = N[Root[Sum[(2Mod[Floor[(t - 1)/2^i], 2] - 1) #^(n - i), {i, 0, n}], root]];
{j,i} = Round[m({Re[z], Im[z]}/1.5 + 1)/2];
If[0 < i <= m && 0 < j <= m, image[[i, j]]++], {root, 1, n}], {t, 1, 2^n}];
ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {0, 4}]

21
Oct
08

### Dodecaplex (120 Cell)

Polychorons are the 4D version of polyhedrons. One way to visualize a polychoron is to apply a 4D to 3D stereographic projection to it. A dodecaplex is a uniform 4D polychoron composed 120 dodecahedral cells. These cells can be divided into 12 rings (Hopf fibrations) of 10 cells each. This picture shows a stereographic projection of 6 rings of the dodecaplex. Each ring is shown in a different color, but only 5 rings are open to direct view because they are wrapped around the 6th ring. I first saw this concept on Matthias Weber’s book page. Click here to download some POV-Ray code.

25
Aug
08

### Cross Section of the Quintic Calabi-Yau Manifold

The left picture shows a 3D cross section of the quintic 6D Calabi-Yau Manifold proposed for String Theory. The right picture shows various degrees of complexity.

Here is some Mathematica code:
(* runtime: 54 seconds, change n for other degrees of complexity *) n = 5; CalabiYau[z_, k1_, k2_] := Module[{z1 = Exp[2Pi I k1/n]Cosh[z]^(2/n), z2 = Exp[2Pi I k2/n]Sinh[z]^(2/n)}, {Re[z1], Re[z2], Cos[alpha]Im[z1] + Sin[alpha]Im[z2]}]; Do[alpha = (0.25 + t)Pi; Show[Graphics3D[Table[ParametricPlot3D[CalabiYau[x + I y, k1, k2], {x, -1, 1}, {y, 0, Pi/2}, DisplayFunction -> Identity, Compiled ->False][[1]], {k1, 0, n - 1}, {k2, 0, n - 1}], PlotRange -> 1.5{{-1, 1}, {-1, 1}, {-1, 1}}, ViewPoint -> {1, 1, 0}]], {t, 0, 1, 0.1}];

12
Aug
08

### Clifford Torus

The Hopf map is a special transformation invented by Heinz Hopf that maps to each point on the ordinary 3D sphere from a unique circle of points on the 4D sphere. Taken together, these circles form a fiber bundle called a Hopf Fibration. If you apply a 4D to 3D stereographic projection to the Hopf Fibration, you get a beautiful 3D torus called a Clifford Torus composed of interlinked Villarceau circles.

By applying 4D rotations to the Hopf Fibration, you can transform the Clifford Torus into a Dupin cyclide or you can turn it inside-out.

Here is some Mathematica code:
(* runtime: 7 seconds *) HopfInverse[theta_, phi_, psi_] := {Cos[phi/2] Cos[psi], Cos[phi/2]Sin[psi], Cos[theta + psi]Sin[phi/2], Sin[theta + psi]Sin[phi/2]}; Ryw[theta_] := {{1, 0, 0, 0}, {0, Cos[theta], 0, Sin[theta]}, {0, 0, 1, 0}, {0, -Sin[theta], 0, Cos[theta]}}; StereographicProjection[{x_, y_, z_, w_}] := {x, y, z}/(1 - w); Table[Show[Graphics3D[Table[{Hue[(4 phi/Pi - 1)/3], Table[Line[Table[StereographicProjection[Ryw[alpha].HopfInverse[theta, phi, psi]], {psi, 0.0, 2Pi, Pi/18}]], {theta, 0.0, 2 Pi,Pi/9}]}, {phi, Pi/4, 3Pi/4, Pi/4}], PlotRange -> 3{{-1, 1}, {-1, 1}, {-1, 1}}]], {alpha, 0, Pi, Pi/18}];

15
Apr
08

### Moiré Pattern

A Moiré pattern is the interference of two similar overlapping patterns. Here is the Moiré pattern on a twisted IKEA wastepaper basket. The mesh on the wastepaper basket was ray-traced from 100,000 tiny cylinders.

Here is some Mathematica code to plot Moiré contours around radiating lines:
(* runtime: 1.7 seconds *) f[dx_] := Sin[200ArcTan[x - dx, y]]; DensityPlot[f[0.1] - f[-0.1], {x, -1, 1}, {y, -1, 1}, PlotRange -> {0, 1}, PlotPoints -> 275, Mesh -> False, Frame -> False]

Here is some Mathematica code to plot a Moiré pattern from rapidly varying contours of a function:
(* runtime: 0.8 second *) f[z_] := z^3; DensityPlot[Sin[20Pi Abs[f[x + I y]]], {x, -2.5, 2.5}, {y, -2.5, 2.5}, PlotPoints -> 275, Mesh -> False, Frame -> False]

27
Mar
08

### Weierstrass Function

The infamous Weierstrass function is an example of a function that is continuous but completely undifferentiable. $\displaystyle \sum_{k=1}^{ \infty}\frac{sin(\pi{k }^{a}x)}{\pi{k}^{a}}$

(* runtime: 0.7 second *) Plot3D[Sum[Sin[Pi k^a x]/(Pi k^a), {k, 1, 50}], {x, 0, 1}, {a, 2, 3}, PlotPoints -> 100];

08
Nov
07

### Joukowski Airfoil

These animations were created using a conformal mapping technique called the Joukowski Transformation. A Joukowski airfoil can be thought of as a modified Rankine oval. It assumes inviscid incompressible potential flow (irrotational). Potential flow can account for lift on the airfoil but it cannot account for drag because it does not account for the viscous boundary layer (D’Alembert’s paradox). In these animations, red represents regions of low pressure. The left animation shows what the surrounding fluid looks like when the Kutta condition is applied. Notice that the fluid separates smoothly at the trailing edge of the airfoil and a low pressure region is produced on the upper surface of the wing, resulting in lift. The lift is proportional to the circulation around the airfoil. The right animation shows what the surrounding fluid looks like when there is no circulation around the airfoil (stall). Notice the sharp singularity at the trailing edge of the airfoil.

Joukowski Airfoil – C++, 11/8/07

Here is an animation that shows how the streamlines change when you increase the circulation around the airfoil. (Please note: The background fluid motion in this animation is just for effect and is not accurate!) Here is some Mathematica code to plot the streamlines and pressure using Bernoulli’s equation:
(* runtime: 13 seconds *) U = rho = 1; chord = 4; thk = 0.5; alpha = Pi/9; y0 = 0.2; x0 = -thk/5.2; L = chord/4; a = Sqrt[y0^2 + L^2]; gamma = 4Pi a U Sin[alpha + ArcCos[L/a]]; w[z_, sign_] := Module[{zeta = (z + sign Sqrt[z^2 - 4 L^2])/2}, zeta = (zeta - x0 - I y0)Exp[-I alpha]/Sqrt[(1 - x0)^2 + y0^2]; U(zeta + a^2/zeta) + I gamma Log[zeta]/(2Pi)]; sign[z_] := Sign[Re[z]]If[Abs[Re[z]] < chord/2 && 0 < Im[z] < 2y0(1 - (2Re[z]/chord)^2), -1, 1];w[z_] := w[z, sign[z]]; V[z_] = D[w[z, sign], z] /. sign -> sign[z]; << GraphicsMaster; DisplayTogether[DensityPlot[-0.5rho Abs[V[(x + I y)Exp[I alpha]]]^2, {x, -3, 3}, {y, -3, 3}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (If[# == 1, Hue[0, 0, 0], Hue[(5# - 1)/6]] &)],ContourPlot[Im[w[(x + I y)Exp[I alpha]]], {x, -3, 3}, {y, -3, 3}, Contours -> Table[x^3 + 0.0208, {x, -2, 2, 0.1}], PlotPoints -> 100, ContourShading -> False], AspectRatio -> Automatic];

02
Nov
07

### Rayleigh Surface Waves

Rayleigh Surface Waves are associated with ocean waves and earthquakes.

Here is some Mathematica code:

(* 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}];

## 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

• 462,897 hits