User:Sam Derbyshire/Mathematica

Here is some sample Mathematica code I've used to generate some mathematical images for Wikipedia.

Back to my user page, my talk page.

3D plots

The Gamma function.

(*Gradient*)
yellowred[op_: 1] := 
  Blend[{{1, RGBColor[1, 0.15, 0, op]}, {0, 
      RGBColor[1, 0.8, 0, op]}}, #] &;
(*Antialiasing factor*)
aa := 2;
(* The function to plot! *)
f[x_, y_] := Abs[Gamma[x + I y]];
(* Bounds for the plot *)
{mx, Mx, my, My, mz, Mz} := {-5.5, 5, -5, 5, 0, 10};
(* Bounds for the gradient *)
{gz, Gz} := {0, 10};
(* Bounds for ticks/gridlines *)
{tx, Tx, ty, Ty, tz, Tz} := {-5, 5, -5, 5, 0, 10};
(* Number of big/small meshes *)
{nx, ny, nz} := {10, 10, 10};
{Nx, Ny, Nz} := {21, 20, 20};
(* Small/big line thicknesses *)
{w, W} := {0.001, 0.002};
(* Plot points and recursion *)
pp := 50;
recursion := 3;
(* Font options *)
fontsize := 6;
fontstyle := "jsMath-cmr10";
fontcolor := RGBColor[0.3, 0.6, 0.8];
font := Directive[(FontFamily -> fontstyle), fontsize, fontcolor];
(* Labels for axes *)
labels := {"x", "y", "z"};
(* Viewpoint for the 3D plot *)
viewpoint := {-5, -5, 3};
signs := Map[-Sign[#] &, viewpoint];
(* Aspect ratio *)
aspects := {1, 1, 0.8};
(* Lighting *)
lighting := {{"Ambient", RGBColor[0.9, 0.9, 0.9]}, {"Point", 
    RGBColor[0.15, 0.15, 0.15], {(mx + Mx)/2, (my + My)/2, 
     2 Mz + 1}}};

(* Gridlines *)
xfcs := Join[Table[i, {i, mx, tx, (Mx - mx)/Nx}], 
   Table[i, {i, tx, Tx, (Mx - mx)/Nx}], 
   Table[i, {i, Tx, Mx, (Mx - mx)/Nx}]];
xtcks := Table[i, {i, tx, Tx, (Tx - tx)/nx}]; 
yfcs := Join[Table[i, {i, my, ty, (My - my)/Ny}], 
   Table[i, {i, ty, Ty, (My - my)/Ny}], 
   Table[i, {i, Ty, My, (My - my)/Ny}]];
ytcks := Table[i, {i, ty, Ty, (Ty - ty)/ny}]; 
zfcs := Join[Table[i, {i, mz, tz, (Mz - mz)/Nz}], 
   Table[i, {i, tz, Tz, (Mz - mz)/Nz}], 
   Table[i, {i, Tz, Mz, (Mz - mz)/Nz}]];
ztcks := Table[i, {i, tz, Tz, (Tz - tz)/nz}]; 
st1[l_] := 
  Map[{#, Directive[Thickness[W], RGBColor[0.2, 0.2, 0.2, 0.3]]} &, l];
st2[l_] := 
  Map[{#, Directive[Thickness[w], RGBColor[0.4, 0.4, 0.4, 0.2]]} &, l];
ticky[l_, ts_: 0.01, fs_: fontsize] := 
  Map[{#, Style[#, (FontFamily -> fontstyle), fs], {ts, 0}} &, l];
xgrid := Join[st1[xtcks], st2[xfcs]]; ygrid := 
 Join[st1[ytcks], st2[yfcs]]; zgrid := Join[st1[ztcks], st2[zfcs]];

(* The 3D surface plot *)
surface[fn_, grad_] := Plot3D[fn[x, y], {x, mx, Mx}, {y, my, My},
  PlotPoints -> pp,
  PlotRange -> {{mx, Mx}, {my, My}, {mz, Mz}},
  MaxRecursion -> recursion,
  PlotRange -> {mz, Mz},
  MeshFunctions -> {#1 &, #2 &, #3 &},
  BoxRatios -> aspects,
  FaceGrids -> {{{signs[[1]], 0, 0}, {ygrid, zgrid}}, {{0, signs[[2]],
       0}, {xgrid, zgrid}}, {{0, 0, signs[[3]]}, {xgrid, ygrid}}},
  Mesh -> {xtcks, ytcks, ztcks},
  (* Note: 
  may need to remove some ticks manually to remove overlaps! *)
  (* Use this: Join[{{ytcks[[1]],"",{0.01,0}}},Drop[ticky[ytcks],1]] *)
  Ticks -> {Join[{{xtcks[[-1]], "", {0.01, 0}}}, 
     Drop[ticky[xtcks], -1]], ticky[ytcks], ticky[ztcks]},
  AxesLabel -> labels, 
  LabelStyle -> 
   Directive[(FontFamily -> fontstyle), fontsize + 2, fontcolor],
  Boxed -> False,(*BoxStyle->{Directive[Black,Opacity[0.9],Thickness[
  W]]},*)
  PlotRangePadding -> None, ClippingStyle -> {{Black, Opacity[0.5]}},
  MeshStyle -> {{Black, Opacity[0.3], Thickness[W]}, {Black, 
     Opacity[0.3], Thickness[W]}, {White, Opacity[0.2], Thickness[W]}},
  Lighting -> lighting,
  ColorFunction -> (grad[(#3 - gz)/(Gz - gz)] &) , 
  ColorFunctionScaling -> False,
  ViewPoint -> viewpoint
  ]

(* The 2D contours *)
contours[fn_, grad_] := ContourPlot[fn[x, y], {x, mx, Mx}, {y, my, My},
   		    Frame -> None,
   		    Axes -> None,
   		    PlotRange -> {{mx, Mx}, {my, My}, {mz, Mz}},
   		    PlotRangePadding -> None,
   	             BoundaryStyle -> None, ClippingStyle -> None, 
   ExclusionsStyle -> None,
   		    ContourShading -> None,
   		    PlotPoints -> pp,		         
   		    ContourStyle -> 
    Table[{grad[(i - 1)/(Nz - 1)], Thickness[W]}, {i, 0, Nz}],
   		    Contours -> Table[i, {i, mz, Mz, (Mz - mz)/Nz}]
   		    ];
(* Hack to project down...*)
proj[l_, fn_] := 
 GraphicsComplex[l[[1]] /. {x_?AtomQ, y_?AtomQ} :> {x, y, fn[x, y]}, 
  l[[2]]]
contours3d[fn_, grad_] := 
  Graphics3D[proj[Graphics[contours[fn, grad]][[1]], mz &]];

(* The plot itself. *)
plot := Show[surface[f, yellowred[]],
  		    contours3d[f, yellowred[]]
  		     ];

Image[ImageResize[
  Rasterize[plot, "Image", ImageResolution -> 150 aa, 
   Background -> None], Scaled[1/aa]], Magnification -> 1]