Jump to content

User:Sam Derbyshire/Mathematica

From Wikipedia, the free encyclopedia

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

[edit]
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]