title: Is this the most effective and efficient programming style for prototyping (a ray tracer) in Mathematica?

text: I am trying out Mathematica as a prototyping tool. As an initial exercise I have put together a brute force ray tracer, with a view to using Mathematica's built in probability distribution functions to test out some fancy shaders. (My starting point was the F# ray tracer available http://code.msdn.microsoft.com/ParExtSamples.) As this is a first attempt at using Mathematica, other than as a fancy calculator, I would welcome some **guidance or critique on whether this coding style is going to be effective**. (You will detect from the code that I am accustomed to an object-oriented paradigm.) My own observations, in no particular order, are:   1. It's an enjoyable way to work, as a lot of the busy type definitions and additional syntax required by other languages is not needed here; the code is quite compact and the intention fairly readable (IMHO).   2. I have built the ray tracer from scratch, without using the **_built-in graphics primitives_** as that would defeat the learning purpose of the exercise and I would probably not be able to roll my own shaders. As a simple example, addition is not defined for RGBColor, for example.   3. Probably as a consequence, the tracer is very slow, even if I use all the 4 kernels that are available to me. I dare say that the existing code could be speeded up considerably by removing all the type / pattern matching, but the code would be much harder to follow, I think, and certainly much harder to debug. My approach seems to mean that I cannot take advantage of Compile optimization, so far as I can see. If I wanted a fast ray tracer I would write it in C++, but I wonder whether there are any **easy optimizations that I have missed** (and that don't involve mangling what I have to the extent that the prototyping benefits of Mathematica -- easy refactoring -- would be lost). However, I am surprised that even with 4 kernels, the cpu of this 4-core (8-core, if you count hyperthreading) never maxes out. For example, is it worth changing some of the `:=` to `=`? Things would be even slower if I aliased the results by averaging 4 adjacent traces, for example. **Real time ray tracing** seems out of reach.   4. Ray tracing results in a lot of "corner cases" that are dealt with naturally with **IEEE maths**. Unfortunately Mathematica does produce +/-Infinity for +/-1/0, for example, so some of the code should really be extended to treat those cases properly.   5. It would be great if Mathematica had more built-in **_vector algebra_** so that I could write the equations defining the objects and the rays involved and getting Mathematica to calculate the ray intersection points that are at the heart of the ray tracer. As things stand, `Reduce` and `Solve` have not helped me to find better intersection algorithms, producing either nothing at all or something large an unintelligible, depending on how I posed the problem. Anyway, here is what you get after 110s from `raytrace[400, 400, basescene, 6]`: ![Here is what you get after 110s from raytrace\[400, 400, basescene, 6\]](http://i.stack.imgur.com/AZE08.png) ... and here is the code:                   (* Colour helpers *)         black = {0., 0., 0.};         darkgrey = {.25, .25, .25};         grey = {.5, .5, .5};         white = {1., 1., 1.};         background = black;         defaultcolor = black;              brightness[{r_, g_, b_}] = Mean[{r, g, b}];               scale[k_, c: {r_, g_, b_}] = k * c;              zero = {0.,0.,0.};              (* Mainly for reference; pattern matching normally used instead *)              ray /: start[ray[s_, d_]] = s;         ray /: dir[ray[s_, d_]] = d;              camera /: pos[camera[p_, l_]] = p;         camera /: lookat[camera[p_, l_]] = l;         camera /: forward[camera[p_, l_]] = Normalize[l - p];         camera /: down[camera[p_, l_]] = {0., -1., 0.};         camera /: right[c : camera[p_, l_]] := 1.5 * Normalize[Cross[forward[c], down[c]]];         camera /: up[c : camera[p_, l_]] := 1.5 * Normalize[Cross[forward[c], right[c]]];                   light /: pos[light[p_, c_]] = p;         light /: color[light[p_, c_]] = c;              scene /: things[scene[t_, l_, c_]] = t;         scene /: lights[scene[t_, l_, c_]] = l;         scene /: camera[scene[t_, l_, c_]] = c;              surface /: diffuse[surface[d_, s_, re_, ro_]] = d;         surface /: specular[surface[d_, s_, re_, ro_]] = s;         surface /: reflect[surface[d_, s_, re_, ro_]] = re;         surface /: roughness[surface[d_, s_, re_, ro_]] = ro;              intersection /: thing[intersection[t_, r_, d_]] = t;         intersection /: ray[intersection[t_, r_, d_]] = r;         intersection /: dist[intersection[t_, r_, d_]] = d;         miss = intersection[nothing, ray[zero, zero], Infinity];              sceneobject /: surface[sceneobject[s_, i_, n_]] = s;         sceneobject /: intersect[sceneobject[s_, i_, n_]] = i;         sceneobject /: normal[sceneobject[s_, i_, n_]] = n;              sphere /: center[sphere[c_, r_, s_]] = c;         sphere /: radius[sphere[c_, r_, s_]] = r;         sphere /: surface[sphere[c_, r_, s_]] = s;         normal[sphere[center_, _, _], pos_] = Normalize[pos - center];              plane /: normal[plane[n_, o_, s_]] = n;         plane /: offset[plane[n_, o_, s_]] = o;         plane /: surface[plane[n_, o_, s_]] = s;         normal[plane[n_, _, _], _] = n;              (* Axis-aligned bounding box *)         (* TODO: not yet used; integrate into tracer *)         box /: lowerb[box[l_, u_]] := l;         box /: upperb[box[l_, u_]] := u;              extendby[box[l_, u_], pt_] :=             box[MapThread[Min, {l, pt}], MapThread[Max, {u, pt}]];         size[box[l_, u_]] :=             u - l;         majoraxis[b : box[l_, u_]] :=             Ordering[size[b], -1];                    (* TODO: This does not work for cases where dir has 0 compnent as Mathematic returns ComplexInfinity,            not +/-Infinity for +/-1/0 *)         intersectboxQ[b : box[l_, u_], r : ray[start_, dir_]] :=             Module[ {tl = (l - start) / dir, tu = (u - start) / dir, tmin, tmax},                      (* Swap u[i] and l[i] if dir[i] < 0 to avoid                     erroneous result because 0 == -0 *)                 tmin = Max[MapThread[Min, {tu, tl}]];                 tmax = Min[MapThread[Max, {tu, tl}]];                 Not[tmax < 0 && tmin > tmax];             (* Use Not to cover some Infinity comparisons *)                  (*  Which[                 tmax < 0, False, (* Intersection at t = tmax, but it's behind us *)                 tmin > tmax, False, (* No intersection *)                 True, True (* Interesection at t = tmin *)                 ]*)             ];                        intersect[s : sphere[center_, radius_, _], r : ray[start_, dir_], i : intersection[_, _, currentdist_]] :=             Module[ {eo = center-start, v, dist, disc},                 v = eo.dir;                 dist = If[ v < 0.,                            0.,                            disc = radius * radius - (eo.eo - v * v);                            If[ disc < 0.,                                0.,                                v - Sqrt[disc]                            ]                        ];                 If[ dist == 0. || dist > currentdist,                     i,                     intersection[s, r, dist]                 ]             ];                   intersect[p : plane[norm_, offset_, _], r : ray[start_, dir_], i : intersection[_, _, currentdist_]] :=             Module[ {denom =  norm . dir, candidatedist},                 If[ denom >= 0.,                     i,                     candidatedist = (norm . start + offset) / (-denom);                     If[ candidatedist > currentdist,                         i,                         intersection[p, r, candidatedist]                     ]                 ]             ];                   testray[ray_, scene_] :=             dist[Fold[intersect[#2, ray, #1]&, miss, things[scene]]];                   traceray[ray_, scene_, depth_, maxdepth_] :=             shade[Fold[intersect[#2, ray, #1]&, miss, things[scene]], scene, depth, maxdepth];                   shade[miss, _, _, _] :=             background;              shade[intersection[thing_, ray[start_, dir_], dist_], scene_, depth_, maxdepth_] :=             Module[ {pos = dist * dir + start, n, reflectdir, naturalcolor, reflectedcolor},                 n = normal[thing, pos];                 reflectdir = dir - 2. * n . dir * n;                 naturalcolor = defaultcolor + getnaturalcolor[thing, pos, n, reflectdir, scene];                 reflectedcolor = If[ depth >= maxdepth,                                      grey,                                      getreflectioncolor[thing, pos + (0.001*reflectdir), n,                                           reflectdir, scene, depth, maxdepth]                                  ];                 naturalcolor + reflectedcolor             ];              getreflectioncolor[thing_, pos_, n_, rd_, scene_, depth_, maxdepth_] :=             reflect[surface[thing]][pos] *             traceray[ray[pos, rd], scene, depth + 1, maxdepth];              getnaturalcolor[thing_, pos_, n_, rd_, scene_] :=             Module[ {addlight, normraydir = Normalize[rd], howrough = roughness[surface[thing]]},                 SetAttributes[addlight, Listable];                 addlight[light[p_, c_]] :=                     Module[ {ldis = p - pos, livec, neatisect, isinshadow, illum, lcolor, spec, scolor},                         livec = Normalize[ldis];                         neatisect = testray[ray[pos, livec], scene];                         isinshadow = neatisect <= Norm[ldis];                         If[ isinshadow,                             defaultcolor,                             illum = livec . n;                             lcolor = If[ illum > 0.,                                          illum * c,                                          defaultcolor                                      ];                             spec = livec . normraydir;                             scolor = If[ spec > 0.,                                          (spec ^ howrough) * c,                                          defaultcolor                                      ];                             diffuse[surface[thing]][pos] * lcolor + specular[surface[thing]][pos] * scolor                         ]                     ];                 defaultcolor + Total[addlight[lights[scene]]]             ];                         raytrace[screenwidth_ : 64, screenheight_ : 64, scene_ : basescene, maxdepth_ : 1] :=              Module[ {getpoint},                  getpoint[x_, y_, camera_] :=                      With[ {recenterx = (x - screenwidth / 2.) / (2. * screenwidth),                       recentery = -(y - screenheight / 2.) / (2. * screenheight)},                          Normalize[forward[camera] + recenterx * right[camera] + recentery * up[camera]]                      ];                  Image[ParallelArray[traceray[ray[pos[camera[scene]], getpoint[#2-1, #1-1, camera[scene]]],                       scene, 0, maxdepth]&, {screenheight, screenwidth}]] // AbsoluteTiming              ];                    (* Harness *)              (* surface diffuse, specular reflect, roughness *)         uniformsurface[diffuse_, specular_, reflect_, roughness_] = surface[diffuse&, specular&, reflect&, roughness];              shiny = uniformsurface[white, grey, .7, 250.];         matteshiny = uniformsurface[white, darkgrey, .7, 250.];              checkerboard =             surface[             If[ OddQ[Floor[#[[3]]] + Floor[#[[1]]]],                 white,                 black             ]&,              white &,              If[ OddQ[Floor[#[[3]]] + Floor[#[[1]]]],                 .1,                 .7             ]&,              150.];                         basescene = scene[{                sphere[{0., 1., -.25}, 1., shiny],                    sphere[{-.5, 1.3, 1.5}, 0.5, matteshiny],              plane[{0., 1., 0.}, 0., checkerboard]},               {light[{-2., 2.5, 0.}, {.5,.45,.41}], light[{2.,4.5,2.},{.99,.95,.8}]},               camera[{2.75, 2.0, 3.75}, {-.6, .5, 0.}]];      * * * To illustrate where I am going with this, here is an example of how one could generate a mesh for a more complex object (a polysphere):                   polyspherepoints[rad_Real, divs_Integer] :=             With[ {u = -(Pi/2.), v = -Pi,                du = Pi/divs, dv = (2.*Pi)/divs},                 rad *                   Flatten[Table[{Cos[du*i + u]*Cos[dv*j + v], Sin[du*i + u],                      Cos[du*i + u]*Sin[dv*j + v]}, {j, 0, divs}, {i, 0,                  divs}], 1]             ];              (* Put the polygon vertices in the right order *)         ordervertices[{{a_, b_}, {c_, d_}}] :=             {a, b, d, c};              orderverticestotriangeles[{{a_, b_}, {c_, d_}}] :=         {{a, b, d}, {a, c, d}}              (* Generate a list of (polyspherepoint) vertice numbers,             partition them cyclically, and then into quads, and associate them            with Polygons *)              polyspheremeshtriangles[rad_Real, divs_Integer] :=             Normal @ GraphicsComplex[polyspherepoints[rad, divs],                Map[Polygon,                 Map[orderverticestotriangeles,                  Partition[Partition[Range[(divs+1)^2], divs+1], {2, 2}, 1, 1], {2}], 1]];               polyspheremeshtriangles[rad_Real, divs_Integer] :=             Normal @ GraphicsComplex[polyspherepoints[rad, divs],                Map[Polygon,                 Map[orderverticestotriangeles,                  Partition[Partition[Range[(divs+1)^2], divs+1], {2, 2}, 1, 1], {2}], 1]];      (It would have been satisfying to use some of the geometric transform functions built into Mathematica to generate the vertices, but life was too short.) And here is what `Graphics3D @ polyspheremeshtriangles[1., 8]` generates: ![PolySphere, of trianges](http://i.stack.imgur.com/RU6dD.png)


title: How to implement Donald Knuth's up-arrow notation

text: How does one enter very large numbers in up-arrow notation? Up-arrow notation was created by Donald Knuth to write very large numbers in it iterated exponentiation form, for example 6↑↑3 = 6^6^6.


title: Prime factor counting function

text: Is there any way I can speed up this prime factor counting function? (I am looking for all numbers in a range with 3 prime factors (counted with multiplicity).)               Omega3[n_Integer] := \[Not] FreeQ[PrimeOmega[n], _?(# == 3 &)]     Omega3Count[n_] := Count[Range@n, _?Omega3]


title: Mathematica Library

text: I created a Mathematica packagr and then I used the Needs function to add the library to my .nb file, but when use a function in my .nb file and the function is built in my library it comes up with the functions output multiplied by "mylibraryname`Private`". My Package:               BeginPackage["Afak10`"]      Unprotect @@ Names["Afak10`*"];      ClearAll @@ Names["Afak10`*"];      f::usage = "f[x]"      Begin["Private"]      f[x_] := Module[{}, x^2 + s];      End[]      Protect @@ Names["Afak10`*"];      EndPackage[]      Output:               Afak10`Private`s + x^2      What if I define f[x] as:               f[x_] := Module[{asdfs},asdfs];      The Output is:               Afak10`Private`asdfs$4980


title: Does Mathematica know if an integral is convergent or not?

text: Consider the following integrals that I asked _Mathematica_ to do,               Assuming[m >= 1 && q >= 1 && k >= 1 && m ∈ Integers && k ∈ Integers,        Integrate[Exp[-2 π[m q Sqrt[λ + 1/4] + k Sqrt[λ]]], {λ, 0, ∞}]]          Integrate[Tanh[π Sqrt[λ]] Log[1 - Exp[-2 π q Sqrt[λ + 1/4]]], {λ, 0, ∞}]         for various different positive values of $m,k,q$ or in general In each of these cases, _Mathematica_ just rewrote back the integral doing nothing. Does this mean that the integral converges, but _Mathematica_ doesn't know how to do it? * * * It would be an added bonus if someone can point out if and how these integrals can be done!


title: ContourPlot - unequal contour spacing

text: how can i specify the contour spacing or how the color range is distributed amongst the function values? Let's say i plot some function depending on x and y. Function values differ between [0, 1]. I want that mathematica makes ...let's say 10 contours for the function value range [0,0.9] and another 10 contours for the rest of the range [0.9,1]. Further more the last Contour of the [0,0.9] range is supposed to be thicker, in order to show that from here on the contour spacing is finer. thx for your help thats alot so far for your help, my actual question has been answered already. but a follow-up question has arisen for me since the implementation of your codes does not work with my actual function and i don't know why. I have purged that function and the functions it calls from most variables and replaced them by constants for simplicity reasons. I hope it's ok so.               eislexpl[epsilon_, n2_, x_] :=  28 (epsilon)^2 - 1 +      ((-2 Log[E^(1/2) x])/x^2 + (2.7(0.3 - epsilon))/x +      (10 epsilon^2 n2^(3/2))/x^(3/2));          eislminx2[epsilon_?NumericQ, n2_?NumericQ] :=      (minx = FindArgMin[{eislexpl[epsilon, n2, x], 0 < x}, x];     eislexpl[epsilon, n2, minx]);          ueislexpl[epsilon_, h_, n1_, n2_] :=      n1 (40  epsilon^2 - 1) + (1 -1.27 ) (1 + (n1 - 1) HeavisideTheta[1 - n1] -      (-1 + Exp[-(n1 - 1)]) HeavisideTheta[n1 - 1]) + n2  eislminx2[epsilon, n2] +      (h - n1 - n2) (28  epsilon^2 - 1);      ueislexpl[epsilon_, h_, n1_, n2_] is the function i want to plot with ContourPlot. If i use :               ContourPlot[      ueislexpl[0.1, 5, n1, n2], {n1, 0, 5}, {n2, 0, 5 - n1},       Contours -> (Function[{min, max},       Module[{r1, r2},       Join[r1 = Range[min, .1 max, (.1 max - min)/10],            r2 = Range[.1 max + .01, max, (max - .1 max - .01)/10]]]]),       ContourStyle -> (Join @@ {Table[Thick, {Length@r1}],           Table[Thin, {Length@r2}]})]      It should plot a contourplot with 10 divisions for the value range [min,0.1 max] and 10 more divisions for [0.1 max, max]. But my result is 3 divisions. What do i have to change?


title: Differential Equations: Solving a second order ODE with DSolve

text: I have the following equation               h''[η] + (2 a'/a) h'[η] + (k)^2 h[η] == 0     2 a'/a = 1.3551      with boundary conditions               h[0] == 1, h'[0] == 0      where `k` is some constant factor. I have to plot `h[η]` @ `η = c` versus `k η` @ `η = c`. So far this is what I have done:               sol = DSolve[{h''[η] + (2 a'/a) h'[η] + (k)^2 h[η] == 0, h[0] == 1, h'[0] == 0},                 h[η], η] /. η -> 1.48* 0.0195          Plot[h[0.02886`] /. sol, {k, 0, 1000}]      As a note: I was expected to plot vs `k η` but since _Mathematica_ is refusing to do so I plotted vs `k` only. The result looks like this > ![Mathematica graphics](http://i.stack.imgur.com/dOBx7.png) However the solution should be damped with 3 oscillations only as shown below:![Expected plotsshown by dotted lines](http://i.stack.imgur.com/I58MS.png)


title: Remove CellMargins option from many cells

text: I have a very old notebook in which many cells had the `CellMargins` option set manually. What's the easiest way to remove all those manual option settings?


title: How can I prevent graphics scaling with CDF notebook size?

text: Recently we saw how to produce xkcd-style graphs. I evaluated the same code and deployed the result to CDF. Everything is working fine except that, in the CDF notebook, the graph is scaled according to the size of the window. To start with, we have: ![Original size](http://i.stack.imgur.com/BMoxW.png) But, on resizing the window, the graphic also gets smaller: ![Resized window](http://i.stack.imgur.com/a0Bqp.png) I don't want the graph to be scaled with the window. How can I prevent this?


title: Perimeter and area are positive integers

text: In Geometry 3D, How can I find the vertices with integer coordinates of a triangle whose perimeter and area are positive integers with _Mathematica_? Suppose its vertices $(x,y,z)$ has coordinates belong to interval $[-50,50]$ If the triangle with three sides 9, 10, 17, I choose a vertex is $(1,2,3)$ and tried               a = {x, y, z};     b = {x1, y1, z1};     c = {1, 2, 3};     {a, b, c} /.       Solve[{SquaredEuclideanDistance[a, b] == 9^2,         SquaredEuclideanDistance[a, c] == 10^2,         SquaredEuclideanDistance[c, b] == 17^2, -50 <= x <= 50, -50 <= y <=          50, -50 <= x1 <= 50, -50 <= y1 <= 50}, {x, y, z, x1, y1, z1},        Integers]      I tried with perimeter. Time is too long.               a = {x1, y1, z1};     b = {x2, y2, z2};     p = Norm[a] + Norm[b] + Norm[a - b];     Solve[{p == k,        12 <= k <= 50, -5 <= x <= 5, -5 <= y <= 5, -5 <= z <= 5,        5 <= x1 <= 5, -5 <= y1 <= 5, -5 <= z1 <= 5}, {x, y, z, x1, y1,        z1, k}, Integers]      I used Maple and found some triangles. For example $A(-6,1,2)$, $B(-9,1,2)$, $C(-9,1,6)$ or $A(1, 2, 3)$, $B(13, 21, 51)$, $C(49, 18, 15)$; The answer of Maple http://www.mapleprimes.com/questions/200319-Perimeter-And- Area-Of-A-Triangle-Are I love this answer. * * *


title: What do the options of SmoothKernelDistribution do?

text: The function `SmoothKernelDistribution` has three options that are not described in too much detail in the Mathematica's help window. **InterpolationPoints** : What is interpolated by Mathematica in function `SmoothKernelDistribution`? **MaxMixtureKernels** : As far as my limited knowledge of kernel estimates goes, there will be as many kernels as grid points on my discretized domain. If my data lives on $x$, and I am attempting to approximate probability density function $f(x)$ then I may choose to do so at N evenly-spaced discrete points along $x$, e.g. $x_i, i=1,2,...,N$. **MaxRecursion** : What recursion does this option pertain to?


title: Comfortable Edge Labeling of Undirected Graph

text: I'm working on a program that finds special [called "interval"] edge-colorings for graphs. Output is as follows:               Coloring = {1 \[UndirectedEdge] 10 -> 1, 1 \[UndirectedEdge] 2 -> 2,                   9 \[UndirectedEdge] 10 -> 2, 1 \[UndirectedEdge] 9 -> 3,                   2 \[UndirectedEdge] 10 -> 3, 1 \[UndirectedEdge] 3 -> 4}      Now I want to visualize the coloring. I use `EdgeLabels->Coloring`, but it writes the labels in the middle, which is not comfortable, especially in case of symmetric graphs like Complete Graphs. I want to display 2 labels on both ends of each edge. Right now I can display it only on one end of each edge using the following code:               1 \[UndirectedEdge] 10 -> Placed[1, 1/6]      and I get the following result: ![CompleteGraph\[10\]](http://i.stack.imgur.com/yjsSz.png) Documentation suggests to use `Placed[{1,1},{"Start","End"}]` to display 2 labels on two ends. That particular example works, but when I change `"Start"` or `"End"` to `1/6` or `5/6`, no labels are displayed at all.


title: Coloring each face of a Cuboid with a different color

text: Is there a way to color each face of a `Cuboid` with a different color? I looked in the help but couldn't find anything. I found this Coloring the faces of a 3D object different colors (in Mathematica) which proposes building the cuboid with `Polygon` primitives; any other way?


title: Subplots with connector lines

text: I am looking for advice from people who have more experience in this area on what is the best (simplest, least effort) way to create a graphic like the following: ![Mathematica graphics](http://i.stack.imgur.com/xyJXQ.png) This is a rough mockup made in a drawing program. There is a central graph, surrounded by smaller ones, each of which is showing some information about a point in the main graph. Those points are connected to the subgraphs with lines. Requirements:   * Each plot must be able to have their own axes/frame   * Proper alignments of the connector lines (red dashed lines on the mockup)---I have the coordinates of one end in the coordinate system of the central plot, while the other end must point at the smaller plots.   * Consistent font sizes and line widths (i.e. everything must be 8 pt when printed)   * Vector graphics (I'd like to avoid rasterizing to bitmaps) Possible approaches:   * `GraphicsGrid` with `Epilog` (`GraphicsGrid` seems to be based on `Inset`.)   * Lots of `Inset`s in a graphic (the main issue is aligning the coordinate system of the central plot with that of the whole graphic)   * Learn to use LevelScheme (I didn't use it for anything serious yet, but when I tried it last time it seemed to have issued with alignment). Whenever I start doing something like this, and the details must be accurate, lots of small issues tend to come up. I'd like to know which approach is likely to prove the least troublesome. * * * **Answers summary** The main difficulty was the correct positioning of connector lines. The usual way of including subplots is by using `Inset` (which is also used by `GraphicsGrid`). One endpoint of the lines is in the main graphics coordinate system, while the other is in the central subplot coordinate system. Converting between the two is very difficult and depends on the scaling of graphics. **Heike's solution** uses `FullGraphics` to expand the axes/frames of subplots. Then all subplots can be directly included in the main graphic and scaled to size. There will be a single coordinate system to deal with. **Chris Degnen's solution** uses image processing to align the main graphic coordinate system with the inset coordinate system. It places a red dot at the desired endpoints, rasterizes the graphic, measures the position of the dot, and then uses this information to compose a vector graphic with the connector lines going between these positions. The result is a vector graphic that looks correct only at a certain scale, but can be exported to PDF. The other solutions recommend adding the connector lines manually.


title: Finding equation of tangent line at a point

text: my function is defined as f[x_] = x^1.1 - 2.5 x^.5 I need to find the equation of the tangent line to the curve using the point slope form which would be something like f'[a][x+a]+f[a] a=1 and x=1 in the instance how would i type this?


title: How to show the coordinates of a point and the equation of a line dynamically?

text: I am constructing a diagram that consists of a sine curve and its tangent line at a given argument. However, I want to show the coordinates of the point of tangency and the equation of the tangent line dynamically. My work so far:               f[x_] := Sin@x;     l[x_, a_] := f[a] + f'[a] (x - a);     Manipulate[        Show[         Plot[Sin@x, {x, 0, 2 \[Pi]},           PlotRange -> {-2, 2}, AxesStyle -> Arrowheads[0.03],            AxesLabel -> {"x", "Sin[x]"}, AspectRatio -> Automatic,            PlotLabel -> f[a] <> "+" <> f'[a] ("x-" <> a)],          Plot[l[x, a], {x, a - 2, a + 2}, PlotStyle -> Pink],         Graphics[{PointSize[Medium], Red, Point[{a, Sin@a}]}]],       {a, 0, 2 \[Pi]}]      Unfortunately, what I've done doesn't work. So my question is: how do I get the information I need to make it work. ![enter image description here](http://i.stack.imgur.com/O3TFk.png)


title: Computing eigenvectors and eigenvalues

text: I have a (non-sparse) $9 \times 9$ matrix and I wish to obtain its eigenvalues and eigenvectors. Of course, the eigenvalues can be quite a pain as we will probably not be able to find the zeros of its characteristic polynomial. Actually, what I really want is find the eigenvector belonging to the largest eigenvalue. Would the following code give me this: `Eigenvectors[matrix, 1]`? How does Mathematica do this? Is there some kind of algorithm of which's existence I am unaware that computes the eigenvector belonging to the largest eigenvalue? Can the largest eigenvalue be computed numerically (in modulus or even throw away all the complex ones), and then the eigenvectors pseudo-analytically using `Eigenvectors`? Added: The matrix is symbolical.


title: Decoupling system of differential equations

text: Here I have one task and it is preparation for small exam. I solved it by hand for first case 1), but I need to check it in $Mathematica$ and to try to implement it for both cases 1) and 2) automatically. After decoupling I am looking for solutions of the decoupled system and how to get them and come back to the source variables because I need solutions of them? System of two differential equations of 4th order and I need to transform them in 1) equivalent system of 4 equations of second order with appropriate substitutions and 2) 8 equations of first order. How to do that in Mathematica automatically, in one case 4 equations second order maximum and in another case 8 equations of first order maximum? Example is with known constants a and b. I know that I must use some substitutions to reduce the system. Firstly can be helpful homogenous solution when g1[x] and g2[x] are equal to zero. first eq                 g1[x]==a1*X1[x] + a2*X2[x] + a3*Derivative[2][X1][x] -        a4*Derivative[2][X2][x] + a5*Derivative[4][X1][x]      second eq                 g2[x]==b2*X1[x] + b1*X2[x] - b4*Derivative[2][X1][x] +        b3*Derivative[2][X2][x] + b5*Derivative[4][X2][x]      1) 4 equations with adopted substitutions of second order differential equations. 2) 8 equations of first order I did a) by hand in the form $(\text{x1}=\text{X1}(x)) (\text{x3}=\text{X2}(x)) \left(\text{x2}=\text{X1}''(x)\right) \left(\text{x4}=\text{X2}''(x)\right) $ $ X(x)=\left( \begin{array}{c} \text{x1}(x) \\ \text{x2}(x) \\ \text{x3}(x) \\ \text{x4}(x) \\ \end{array} \right) $ $, A=\left( \begin{array}{cccc} 0 & 1 & 0 & 0 \\ -\frac{\text{a3}}{\text{a1}} & -\frac{\text{a2}}{\text{a1}} & -\frac{\text{a5}}{\text{a1}} & -\frac{\text{a4}}{\text{a1}} \\ 0 & 0 & 0 & 1 \\ -\frac{\text{b5}}{\text{b1}} & -\frac{\text{b4}}{\text{b1}} & -\frac{\text{b3}}{\text{b1}} & -\frac{\text{b2}}{\text{b1}} \\ \end{array} \right)$ $, B=\left( \begin{array}{c} 0 \\ \frac{\text{g1}(x)}{\text{a1}} \\ 0 \\ \frac{\text{g2}(x)}{\text{b1}} \\ \end{array} \right)$ $$ X''(x)=A X(x)+B$$


title: Easier program for period of Fibonacci sequence modulo p

text: For a little project I need to calculate the period of a Fibonacci sequence modulo p, for which p is a prime number. For example, the Fibonacci sequence modulo 19 would be: $$0, 1, 1, 2, 3, 5, 8, 13, 2, 15, 17, 13, 11, 5, 16, 2, 18, 1, 0, 1, 1, 2...$$ As you can see, the sequence repeats after 18 numbers. The formula for a Fibonacci sequence is given by: $f(n)=f(n-1)+f(n-2)$. In my attempt at making a program, I came up with the following:               newrow[a_, b_, c_] := Module[{f, g},       f[n_] := f[n] = f[n - 1] + f[n - 2];       f[0] = a;       f[1] = b;       g[n_] := Mod[f[n], Prime[c]];       L = Table[g[n], {n, 0, Prime[c]^2} ];       number1 = Flatten[Position[L, g[0]]];       number2 = Flatten[Position[L, g[1]]];       number3 = Select[number1, (MemberQ[number2, # + 1] ) &];       period = number3[[2]] - 1;       Return[         period       ]     ]      However, if I want to calculate it for quite some prime numbers at once like this:               Table[newrow[0, 1, i], {i, 1, 50}]      the computer can't handle it. Might there be an easier algorithm for this? I also tried this, but this seems to be even harder for the computer:               anderrijtje[a_, b_, c_] := Module[{f, g},       f[n_] := f[n] = f[n - 1] + f[n - 2];       f[0] = a;       f[1] = b;       g[n_] := Mod[f[n], c];       Return[         For[           k = 1, Not[{g[0], g[1]} == {g[k], g[k + 1]}], k++, {}         ];       k       ]     ]      Thanks in advance for any help.


title: Exported PDF figure looks different in Apple's Preview than in Adobe PDF Reader

text: I'm trying to export figures into pdf, and the results looks good in Preview, but very bad in Adobe pdf reader. I'm using Mathematica version 8 on Mac 10.8.5. The Preview version is 6.0.1 (765.6) and the Adobe PDF Reader version is 11.0.04.63. For a simple three point polygon, it looks good in both pdf viewer               Graphics3D[{Red, Opacity[0.15], EdgeForm[Lighter[Black, 0.8]],        Polygon[{{1, 0, 0}, {1, 1, 1}, {0, 0, 1}}]}, Boxed -> False]     Export["~/Downloads/12345.pdf", %, ImageResolution -> 100]      ![enter image description here](http://i.stack.imgur.com/cfopC.png) but for polygon with more points, it looks very bad in Adobe reader               Graphics3D[{Red, Opacity[0.15], EdgeForm[Lighter[Black, 0.8]],        Polygon@Join[{{-1, 0, 0}},          Table[{-1, t, Cos[t]^2}, {t, 0, \[Pi]/2, 0.01}]]},       BoxRatios -> {1, 1, 0.5}, Boxed -> False, ViewPoint -> Right]     Export["~/Downloads/1234.pdf", %, ImageResolution -> 100]      ![enter image description here](http://i.stack.imgur.com/EqX3C.png) Why they looks different in the two pdf viewer and is there a way to fix it so that pdf figure works in both pdf viewers?


title: DSolve output - understanding how to extract result

text: I'm a little confused with how to use the output from `DSolve`. The following works               s = DSolve[{y'[x] == x, y[0] == 0}, y[x], x];     Plot[y[x] /. s, {x, 1, 10}]      but I don't understand how it works. I am familiar with `/. a -> b` which will replace a with b, but I don't get what `y[x] /. s` is doing?


title: Suggested Settings for Mathematica over X11/SSH Tunnel

text: I am happy to discover that I can run a Mathematica session tunneled over SSH. However, as expected for X11 applications over the network, Mathematica's GUI does not have the responsiveness that I wish it had. I have enabled SSH Compression to the highest level, and I am using ethernet connections where possible, though my cable modem is probably the bottleneck. For various reasons, I am unable to install NX, x2go, xpra, etc. Also, it turns out that remote MathKernel/local UI will not be an option. It seems that some stock X applications (such as evince and gnome-terminal) run acceptably over my current setup. My wish is to trim Mathematica GUI as much as possible, so as to minimize re-draws and the like which would slow things down. I have been through the appearance preference pane, and turned off various "code intelligence" features which seem to cause things to redraw, and lowered graphics quality to their lowest. Also, I have disabled niceties such as the Suggestion Bar. But it is possible that I am missing something that others have found. Which settings/modifications have you found to reduce the latency of Mathematica GUI sessions?


title: Plot a 3D graph while changing the range of the parameter

text: I wish to plot the graph of a surface while being able to change a parameter. My current code is               ContourPlot3D[(Cos[Pi/4]+x Cos[0]Sin[Pi/4]+y Sin[0]Sin[Pi/4])^2 ==       (1+x^2+y^2)Cos[a]^2, {x, -10, 10}, {y, -10, 10}, {z, -10, 10},       {a = 0 .. Pi}, AxesLabel -> {x, y, z}]      I don't want to do it with `Manipulate[]`, since I don't a the graph of the surface of a given "a" between 0 and Pi, but the graph of the surface with "a" from 0 to Pi.


title: How to evaluate Modules in a notebook from another notebook?

text: I would like to create a notebook populated with several modules and be able to evaluate these modules by running a separate notebook. I have attempted several uses of the `NotebookEvaluate` command with no luck. Can someone please explain to me how to do this? More specifically, I would like to have a notebook named "functions.nb" which is populated with several functions, and evaluate this notebook in a notebook named "output.nb" which will call my written functions in "functions.nb", evaluate the functions, and provide the output. As an example, if functions.nb contains:               Add[x0_,y0_]:=     Module[{x=x0,y=y0},     z=x+y     ]      And output.nb contains:               Add[2,5]      The output would be:               z=7


title: RSolve doesn't solve this recurrence equation

text: RSolve[{1 == b ((c k[t - 1]^a - k[t])/(c k[t]^a - k[t + 1])) c a k[t]^(a - 1),              k[1] == b, k[3] == a}, k[t], t]      I tried to read the documentation in Help and searched through Google. But I still can't make Mathematica to solve my recurrence function. I am new to Mathematica. Is there any place where I made mistakes? And is there any way that I can solve this function? I'm using Mathematica 8.0.4 for Mac.


title: Integration of Two Manipulatable Cylinders

text: ## Introduction I created the following code to simulate the many possible interactions between a cylindrical protein crystal and a x-ray beam during serial femtosecond crystallography.               Manipulate[       Graphics3D[         {          { EdgeForm[None], Directive[Yellow, Opacity[0.4], Specularity[White, 20]],             GeometricTransformation[Cylinder[{{-h + x, y, z}, {h + x, y, z}}, 1*10^-6],             RotationTransform[α Pi/2, {1, 0, 0}].RotationTransform[β Pi/2, {0, 1, 0}]]},          { EdgeForm[None], Directive[Blue, Opacity[0.3], Specularity[White, 20]],             Cylinder[{{-1*10^-5, 0, 0}, {1*10^-5, 0, 0}}, r]} },              Axes -> True, AxesLabel -> {"X", "Y", "Z"} ],            {{α, 0}, -1, 1}, {{β, 1}, -1, 1},        {{r, 1*10^-7}, 1*10^-7, 6*10^-6} ,        {{h, 5*10^-6}, 1*10^-6, 10*10^-6},       {{x, 0}, -1*10^-5, 1*10^-5},       {{y, 0}, -1*10^-5, 1*10^-5},        {{z, 0}, -1*10^-5, 1*10^-5} ]      ## Problem I would like to be able to display the volume of the intersection of the two objects. As far as I have looked, this doesn't seem to be possible using only `Cylinder` and `Graphics3D`. If there is a way to find this volume using only `Graphics` objects that someone knows of, please enlighten me. Also, the use of `RegionPlot3D` and `NIntegrate` as demonstrated by this example looks promising.               RegionPlot3D[x^2 + y^2 < 1, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]     NIntegrate[Boole[x^2 + y^2 < 1 ], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]      But my problem with using these functions is I don't know how to apply the `Manipulate` function to these raw inequalities in order to get the geometric transforms, variable height, translation, and change in beam radius like I have in my `Graphics3D` model. Really all I care about it that it is manipulatable like the `Graphics3D` model I have shown and that the volume of the intersection is displayed as the model is manipulated. I greatly appreciate any help on this. ## Solution               cyl[{x_, y_, z_}, r_, h_] := y^2 + z^2 < r^2 && -h < x < h     xform[x_, y_, z_, a_, b_] := RotationTransform[a Pi/2, {1, 0, 0}].RotationTransform[b Pi/2, {0, 1, 0}].TranslationTransform[{x, y, z}]          Manipulate[      Graphics3D[       {        {EdgeForm[None], Directive[Yellow, Opacity[0.4], Specularity[White, 20]],        GeometricTransformation[ Cylinder[{{-h, 0, 0}, {h, 0, 0}}, 1*10^-6], xform[x, y, z, a, b]]},        {EdgeForm[None], Directive[Blue, Opacity[0.3], Specularity[White, 20]], Cylinder[{{-1*10^-5, 0, 0}, {1*10^-5, 0, 0}}, r]}      },     PlotLabel -> Chop[10^21 NIntegrate[       Boole[cyl[{x0, y0, z0}, r, 1*10^-5] &&          cyl[InverseFunction[xform[x, y, z, a, b]][{x0, y0, z0}],           1*10^-6, h]], {x0, -2*10^-5, 2*10^-5}, {y0, -2*10^-5,         2*10^-5}, {z0, -2*10^-5, 2*10^-5},       Method -> {"MultidimensionalRule", "Generators" -> 9},        MaxRecursion -> ControlActive[2, 4], AccuracyGoal -> 2,        PrecisionGoal -> 2]],      Axes -> True,      AxesLabel -> {"X", "Y", "Z"}],      {{a, 0}, -1, 1},      {{b, 1}, -1, 1},      {{r, 1*10^-7}, 1*10^-7, 6*10^-6} ,      {{h, 5*10^-6}, 1*10^-6, 10*10^-6},      {{x, 0}, -1*10^-5, 1*10^-5},      {{y, 0}, -1*10^-5, 1*10^-5},      {{z, 0}, -1*10^-5, 1*10^-5}     ]


title: Symbolically solving a system of nonlinear equations

text: I am trying to solve the following system of equations and did the following               Solve[x (1 - (y/(1 + y))) == ((1 - alpha) A y^alpha -            alpha A y^(alpha - 1))/sigma1 &&        1 - x == (0.0725 - alpha A y^(alpha - 1))/sigma2, {x, y}]      I want to solve for x and y symbolically, and then try to do the following. First, I want to find the numerical solution for the above system, given the values for `sigma1, sigma2, alpha, A`. Second, I want to assess whether the derivatives $\frac{\partial y}{\partial \sigma_i}$ and $\frac{\partial z}{\partial \sigma_i}$ are positive or negative by assigning values to the relative parameters. I am using _Mathematica_ 10, and once i tried to run the code above, i got this message: > Solve::inex: Solve was unable to solve the system with inexact coefficients > or the system obtained by direct rationalization of inexact numbers present > in the system. Since many of the methods used by Solve require exact input, > providing Solve with an exact version of the system may help. >> Any suggestions ?


title: NonlinearStateSpaceModel: Trying to parse the documentation

text: I am trying to follow the documentation on `NonlinearStateSpaceModel` here, and the first basic example throws me for a loop. It says it's defining the nonlinear system $x_1'(t)=u(t)+x_1(t)x_2(t)$ $x_2'(t)=u(t)x_2(t)$ with output $y(t)=x_1(t)$ The example call is written               nsys = NonlinearStateSpaceModel[{{u-x1 x2, u x2+1}, {x1}}, {x1, x2}, u]      This appears to be the arity-3 overload for `NonlinearStateSpaceModel`, which, if right, must fit the template $\textrm{NonlinearStateSpaceModel}[\{f,g\},x,u]$ representing the model $x'(t)=f(x(t), u(t)), y(t)=g(x(t),u(t))$ That means that $x$ is `{x1, x2}` $f$ is `{u - x1*x2, u*x2 + 1}` $g$ is `{x1}` which I can only read as $x_1'(t)=u(t)-x_1(t)x_2(t)$ $x_2'(t)=1+u(t)x_2(t)$ with output $y(t)=x_1(t)$ which doesn't look anything like the system they say they're defining. I've blinked my eyes several times and stared at it for too long. It's either a horrendous cluster of typos or I just don't know how to read these models.


title: Adjusting Density Plot Color Function for Diffraction Simulation

text: I'm trying to simulate Fraunhofer diffraction by a single slit using `DensityPlot` I run the code:               DensityPlot[(Sinc[β])^2, {β, -6 π, 6 π}, {y, -6 π, 6 π}, PlotPoints -> 200,                 ColorFunction -> GrayLevel, AspectRatio -> 9/16, Frame :> False]      and get: ![Image](http://i.imgur.com/iGkhRtX.png?1) Which is great but I'd like the coloring to be based on a red scale so that the "lightest" areas are bright red which corresponds to a 632.8 nm He-Ne Laser. The following image shows the type of red scale I'd like to mimic: ![Red](http://tsgphysics.mit.edu/pics/Q%20Diffraction/Q2-Single-Slit- Diffraction.jpg) Thanks


title: Replacing variable name in TeXForm

text: I have some macros in a tex-file I'm writing and I would like Mathematica to use the same notation in its output from the TeXForm command. For example, if I write the equation: a = b I would like the output when I use TeXForm on this expression to be: a = \bias That is, I want Mathematica to output the variable name b as \bias with the TeXForm command. How can this be achieved?


title: How to stop the Kernel from running wild?

text: How to abort or stop this **(accidently large defined)** SparseArray production?               s = N[SparseArray[Table[{2^i, 4} -> i, {i, 30}]]]        Alt+. or Alt+, seem to give up. Even your Windows task manager is struggling.


title: Suspecting bug in Parallel Framework

text: I've reduced my actual problem to:               ParallelEvaluate[       Module[{slot=1},         Slot[slot]       ]     ]      gives error messages that I think shouldn't be there. It seems I can work around them with               ParallelEvaluate[       Module[{slot=1},         Slot@@{slot}       ]     ]      and               ParallelEvaluate[       Module[{slot=1,s=Slot},         s[slot]       ]     ]      but I still think this is a bug in the parallel framework. Am I missing something obvious? It seems the numeric check happens with the local Module placeholder values passed to the parallelism framework in a syntactic way, not in a semantic way. But that's wrong, I should be able to pass `Slot[anything]` v8.0.4 and v9 alike, Windows and Linux.


title: Identify columns of a huge matrix just polynomially many rows chosen at random

text: I am interested in the following problem in combinatorics $\Cap$ Probability. Let $\lambda \in \mathbb{N}$ be a parameter. Consider a matrix of $2^\lambda$ rows and $2^\lambda$ columns. Each column differs from every other column in at least superpolynomial (in $\lambda$) rows. Does there exist some $t$ ($\in \mathbb{N}$) that is polynomial in $\lambda$ such that the following holds? Condition: For every column $C$, rows $R_1,...,R_t$, let $Same_{C,R_1,...,R_t}$ denote the set of columns that have the same entries in Rows $R_1,...,R_t$ as the column $C$. Now choose $t$ random rows $R_1,...,R_t$. Can I be assured that, with all but negligible probability over the choice of the rows, there exists $m$ ($\in \mathbb{N}$) that is polynomial in $\lambda$ such that for every column $C$ its corresponding set $Same_{C,R_1,...,R_t}$ contains at most $m$ columns? (Intuitively, I want to identify a column in this matrix -- upto a precision of polynomially many columns -- by choosing at most polynomially many rows at random.). Thanks!


title: How can one compute a 3D Periodic Delaunay Triangulation using DelaunayMesh?

text: Given a non-periodic point set, one can easily tetrahedralize it using the new (in **V10** ) `DelaunayMesh` function. e.g.:               SeedRandom[0]     pts3d = RandomReal[4, {200, 3}];     del = DelaunayMesh[pts3d];     style = MapThread[Style[#1, Directive[#2]] &, {{1, 0, 2}, {{Thin, Purple},                              {PointSize[0.02], Red}, {Opacity[0.5], Green}}}];     HighlightMesh[del, style]      ![Mathematica graphics](http://i.stack.imgur.com/1Lgww.png) Now consider the following periodic point set in 3D:               cavs = {{0., 1.199, 2.53}, {0., 1.265, 2.53}, {0.067, 1.199, 2.53}, {3.263,            1.199, 2.53}, {0.067, 1.265, 2.53}, {3.263, 1.265, 2.53}, {0.133,            1.199, 2.53}, {3.196, 1.199, 2.53}, {0.133, 1.265, 2.53}, {3.196,            1.265, 2.53}, {0.2, 1.199, 2.53}, {3.13, 1.199, 2.53}, {3.196,            1.199, 2.464}, {0.133, 1.265, 2.597}, {0.133, 1.332, 2.53}, {0.2,            1.265, 2.53}, {3.13, 1.265, 2.53}, {0.2, 1.199, 2.597}, {0.266,            1.199, 2.53}, {3.13, 1.199, 2.464}, {0.2, 1.265, 2.597}, {0.2,            1.332, 2.53}, {0.266, 1.265, 2.53}, {3.063, 1.265, 2.53}, {3.13,            1.265, 2.464}, {0.266, 1.132, 2.53}, {3.063, 1.199, 2.464}, {0.2,            1.265, 2.664}, {0.2, 1.332, 2.597}, {0.2, 1.398, 2.53}, {0.266,            1.332, 2.53}, {0.333, 1.265, 2.53}, {3.063, 1.265, 2.464}, {3.063,            1.332, 2.53}, {3.063, 1.199, 2.397}, {0.2, 1.398, 2.597}, {0.266,            1.332, 2.597}, {0.266, 1.398, 2.53}, {0.266, 1.332, 2.464}, {0.333,            1.332, 2.53}, {0.333, 1.265, 2.464}, {0.4, 1.265, 2.53}, {2.997,            1.265, 2.464}, {3.063, 1.265, 2.397}, {3.063, 1.332, 2.464}, {2.997,            1.332, 2.53}, {2.997, 1.199, 2.397}, {3.063, 1.132, 2.397}, {3.063,            1.199, 2.331}, {0.2, 1.465, 2.597}, {0.266, 1.398, 2.597}, {0.266,            1.332, 2.664}, {0.333, 1.332, 2.597}, {0.266, 1.465, 2.53}, {0.333,            1.398, 2.53}, {0.333, 1.332, 2.464}, {0.4, 1.332, 2.53}, {0.333,            1.199, 2.464}, {0.4, 1.265, 2.464}, {0.466, 1.265, 2.53}, {2.997,            1.265, 2.397}, {2.997, 1.332, 2.464}, {3.063, 1.265, 2.331}, {2.997,            1.332, 2.597}, {2.997, 1.398, 2.53}, {2.997, 1.199, 2.331}, {0.266,            1.465, 2.597}, {0.266, 1.398, 2.664}, {0.333, 1.398, 2.597}, {0.4,            1.398, 2.53}, {0.4, 1.332, 2.464}, {0.466, 1.332, 2.53}, {0.333,            1.132, 2.464}, {0.4, 1.199, 2.464}, {0.466, 1.265, 2.464}, {2.93,            1.265, 2.397}, {2.997, 1.332, 2.397}, {2.93, 1.332, 2.464}, {2.997,            1.398, 2.464}, {2.997, 1.398, 2.597}, {2.997, 1.132, 2.331}, {0.266,            1.465, 2.664}, {0.333, 1.465, 2.597}, {0.333, 1.398,            2.664}, {0.466, 1.332, 2.464}, {0.533, 1.332, 2.53}, {0.333, 1.065,            2.464}, {0.466, 1.265, 2.397}, {2.93, 1.332, 2.397}, {2.997, 1.332,            2.331}, {2.863, 1.332, 2.464}, {2.93, 1.132, 2.331}, {0.266, 1.465,            2.73}, {0.533, 1.332, 2.464}, {2.863, 1.332, 2.397}, {2.997, 1.398,            2.331}, {2.863, 1.132, 2.331}, {0.266, 1.465, 2.797}, {0.333, 1.465,            2.73}, {0.533, 1.332, 2.397}, {0.599, 1.332, 2.464}, {2.93, 1.398,            2.331}, {2.863, 1.199, 2.331}, {0.266, 1.465, 2.863}, {0.266, 1.532,            2.797}, {0.599, 1.332, 2.397}, {0.666, 1.332, 2.464}, {2.797,            1.199, 2.331}, {0.266, 1.532, 2.863}, {0.333, 1.532, 2.797}, {0.666,            1.332, 2.397}, {0.666, 1.332, 2.53}, {0.732, 1.332, 2.464}, {2.73,            1.199, 2.331}, {2.797, 1.265, 2.331}, {0.266, 1.598, 2.863}, {0.333,            1.532, 2.863}, {0.666, 1.398, 2.397}, {0.732, 1.332,            2.397}, {0.732, 1.332, 2.53}, {0.732, 1.398, 2.464}, {0.799, 1.332,            2.464}, {2.664, 1.199, 2.331}, {2.73, 1.132, 2.331}, {2.73, 1.265,            2.331}, {0.333, 1.598, 2.863}, {0.732, 1.398, 2.397}, {0.799, 1.332,            2.397}, {0.732, 1.332, 2.597}, {0.799, 1.398, 2.464}, {2.597,            1.199, 2.331}, {2.664, 1.199, 2.264}, {2.664, 1.199, 2.397}, {2.664,            1.265, 2.331}, {2.73, 1.132, 2.264}, {2.73, 1.265, 2.397}, {0.799,            1.398, 2.397}, {0.799, 1.398, 2.53}, {0.799, 1.465, 2.464}, {0.866,            1.398, 2.464}, {2.597, 1.199, 2.397}, {2.664, 1.199, 2.197}, {2.664,            1.265, 2.264}, {2.664, 1.265, 2.397}, {0.799, 1.398,            2.331}, {0.799, 1.465, 2.397}, {0.866, 1.398, 2.397}, {2.597, 1.132,            2.397}, {2.597, 1.265, 2.397}, {2.664, 1.332, 2.264}, {0.799,            1.465, 2.331}, {0.866, 1.398, 2.331}}      Here's what it looks like:               Graphics3D[{Red, PointSize[0.02], Point[cavs]}, BoxRatios -> {1, 1, 1}]      ![Mathematica graphics](http://i.stack.imgur.com/qPNx4.png) We tetrahedralize it naively:               cavdel = DelaunayMesh[cavs];      Visualize using the same `style` from before:               Show[HighlightMesh[cavdel, style], BoxRatios -> {1, 1, 1}]      ![Mathematica graphics](http://i.stack.imgur.com/tzMP5.png) Obviously, this is wrong as it's assuming the end of the box to be the end region and will lead to a larger volume than the true volume. One can re-align the points in a box that gives them the minimum distance. When this is done the true arrangement of the points looks like this: ![Mathematica graphics](http://i.stack.imgur.com/Z0YnL.png) One can then easily use `DelaunayMesh` to tetrahedralize and obtain the following: ![Mathematica graphics](http://i.stack.imgur.com/ZzeH2.png) Which is the true Delaunay triangulation. The Delaunay tetrahedralization of the original periodic point set should look something like this: ![Mathematica graphics](http://i.stack.imgur.com/bvHf1.png) Clearly, this is not quite right, since there are regions below and above the points that should fill up to the box edges and continue on the other side (wrap-around effect). Since we can't give `DelaunayMesh` a `Distance` function AFAIK, my question is, given a set of periodic points in 3D how can one tetrahedralize it using `DelaunayMesh`? **Note** : The length of the box is `3.2629` in each of `x`, `y` and `z` direction. Origin is `(0, 0, 0)` and Minimum image periodic boundary conditions were applied in all directions.


title: Finding maximum or minimum of implicit functions

text: is there any built in function that can be used to find maximum or minimum of implicit functions? For example, if we have the equation $$x^2 + y^2 = (2 x^2 + 2 y^2 - x)^2,$$ then we can visualize the set of all $(x,y)$ making the equation true using `ContourPlot`.               ContourPlot[      x^2 + y^2 == (2 x^2 + 2 y^2 - x)^2, {x, -1, 2}, {y, -1, 1},       AspectRatio -> Automatic]      ![enter image description here](http://i.stack.imgur.com/tOhA4.png) Clearly, $y$ is not a function of $x$ but, in the neighborhood of most points on the graph, a function is _implied_ , i.e. $y$ is implicitly a function of $x$. Is there any built in way to find the maximum and/or minimum value of this function (like what we have for the explicit functions)?


title: ContourPlot Legend with colorbar for arbitrary contour values

text: I have a ContourPlot in which I have specified the contour-values myself as Contours->{0.01, 0.05, 0.1, 0.2, 0.3}. I would like to get the corresponding contour-values to be shown in the color-bar in the Legend. How do I do this? To get the max and min values, I am using the code in the following post which defines the reportColorRange[...] function: ShowLegend values . Part of my question is answered by another post:   How can I label a ListDensityPlot with a color bar? , but probably this only works for a linear spacing of contours? An example code is below.               Needs["PlotLegends`"]     {plTest, colors, range} = reportColorRange[       ContourPlot[xv^2 + yv^2, {xv, 0, 1}, {yv, 0, 1},          Contours -> {0.01, 0.05, 0.1, 0.2, 0.3}, ColorFunction -> "Rainbow"]];     LMin = range[[1]];  LMax = range[[2]];     ShowLegend[plTest, {ColorData["Rainbow"][1 - #1] &, 6,        ToString[Round[LMax, 0.01]], ToString[Round[LMin, 0.01]],        LegendPosition -> {0.6, 0}, BaseStyle -> {FontSize -> 14}}]       I would like the 5 contour-values to appear in the Legend next to the "ticks" which are blank right now. I have left the max and min values shown right now, but I don't necessarily need them. It's sufficient (and probably better) to show only the contour values. Thanks. PS: This is also related to another post of mine: Mathematica 8: ContourPlot finding out max and min function values Sorry for the multiple posts, but I am still converging on what appearance I want for the plot and legend.


title: How to write this without For loop

text: Suppose I have a few lists of numbers and want to exponentiate element-wise, then sum up everything into a polynomial. For example, if I have               a = {1,2,3}     b = {2,4,6}     c = {3,2,1}      I would like to get               poly = xy^2z^3 + x^2y^4z^6 + x^3y^2z      Unfortunately, I am not too familiar with the functional programming style Mathematica is intended to be used for, and I can only write something like               For[i = 1, i < Length[a] + 1, i++,      term = x^Part[a, i]*y^Part[b, i]*z^Part[c, i];     poly = poly + term;]      As you may expect, this is taking quite a long time to evaluate. Would anybody have any thoughts as to how I would rewrite this to make it more efficient?


title: Performance on Binomial Deviate is slow

text: I hope I won't get spanked for reposting. This is related to a previous question discussed here: Performance on Multinomial Deviate is slow. but this version is simpler because it is a binomial rather than multinomial which is my justification for the new question. I've kluged together a real solution for this problem that someone might find useful. I'm hoping someone will improve on the readability and/or performance. The performance is acceptable but still a bit slow compared to C. I think there might be another factor of 2 or 3 speedup to be had. My question is: does anyone have suggestions that will improve on my solution? The problem is that I need a bunch of deviates where the parameters of the distribution are different for each deviate. Here is the simple way to do it which is unacceptably slow.               Nboxes = 8000;     p = RandomReal[{0, .01}, Nboxes];     rt = RandomInteger[{0, 10}, Nboxes];     rp = RandomInteger[{0, 10}, Nboxes];     probs = 1 - p rp;     out = MapThread[         RandomVariate[BinomialDistribution[#1, #2]] &, {rt, probs}]; // Timing      The timing result was over 4.7 seconds on my laptop. I create one below that is about 25 times faster. Now to generate a binomial deviate with parameters n and p you generate n uniform deviates on (0,1) and count how many are < p. The total number of uniform deviates needed then is Total[rt]. Those can be generated in a small fraction of a second. The following code implements a binomial deviate generator and runs in a total of under 0.2 seconds. The deviates are stored in "counts". If anyone can improve on this code for performance, readability, etc, I'm all ears. First I generate all the uniform deviates I need. Then I "box" those up according to what rt is and count the number of deviates in each partition that are less than that boxes p value.               unitdeviates = RandomReal[1, Total[rt]]; // Timing     nt = ConstantArray[0, Nboxes + 1];     Table[nt[[i]] = nt[[i - 1]] + rt[[i - 1]], {i, 2, Length[rt] + 1}]; // Timing          counts = ConstantArray[0, Nboxes];     Do[counts[[nboxes]] = Count[unitdeviates[[nt[[nboxes]] + 1 ;;         nt[[nboxes + 1]]]], _?(# < probs[[nboxes]] &)], {nboxes, 1,Nboxes}] // Timing      Here I show that the mean and variance of the two sets of deviates are consistent               N[Mean[out]]     N[Mean[counts]]     N[Variance[out]]     N[Variance[counts]]


title: Strategies to solve an oscillatory integrand only known numerically

text: I have an integrand that looks like this: ![enter image description here](http://i.stack.imgur.com/1h7vg.jpg) the details of computation are complicated but I only know the integrand numerically (I use `NDSolve` to solve second order ODE). The integrand is not simply the solution of my ODE either; calling the two solutions of my ODE osc1[s], osc2[s] then schematically the integrand I have looks something like exp(-is)[g(s)osc1[s]osc2*[C-s]+f(s)osc2[s]osc2*[C-s]]. The exp bit is only very slowly oscillating over my integration range, it is really osc1,osc2 that give wild oscillation, as a certain parameter they depend on gets larger. **More explicitely**               rstar[r_] := r + 2 M Log[r/(2 M) - 1];     M=1;     rinf=10000;     rH = 200001/100000;     r0 = 10;     wp=40;     ac=wp-8;     \[Lambda][l_] = l (l + 1);          eq[\[Omega]_,l_] := \[CapitalPhi]''[r] + (2 (r - M))/(     r (r - 2 M)) \[CapitalPhi]'[     r] + ((\[Omega]^2 r^2)/(r - 2 M)^2 - \[Lambda][l]/(     r (r - 2 M))) \[CapitalPhi][r] == 0;     init=-0.0000894423075560122420468703835499 +      0.0000447222944185058822813688948339 I;     dinit=-4.464175354293244250869336196691640386266791`30.*^-6 -      8.950483248390306670770345406047835993931665`30.*^-6 I;          osc1 := \[CapitalPhi] /.      Block[{$MaxExtraPrecision = 100},      NDSolve[{eq[1/10, 1], \[CapitalPhi][rinf] ==      init, \[CapitalPhi]'[rinf] == dinit}, \[CapitalPhi], {r, r0,      rinf}, WorkingPrecision -> wp, AccuracyGoal -> ac,      MaxSteps -> \[Infinity]]][[1]];      osc2 is obtained simiarly. Note these are for non problematic params and it will run quick quickly, and not be too badly behaved. The problem I have is that I only know the integrand to maybe 6-12 digits of precision (dp), depending on the parameters. This is computing the NDSolve with a WorkingPrecision of 50-60, AccuracyGoal->42-52 and it takes around 2 hrs. I want to integrate this with `NIntegrate`, but when my parameters are large (and the oscillation is very high) I usually only know the integrand around the 6 dp end of scale, and `NIntegrate` wants a greater `WorkingPrecision` than this otherwise it complains (since oscillation is also getting very large). I can force it to do the integral by making the `WorkingPrecision` higher, but I think this is cheating if I don't believe my integrand any higher than 6 dp? The only ideas I've had so far are to try different rules. Are there any rules people would recommend for doing such oscillatory integrands? So far I've tried `"LevinRule"`, `"ClenshawCurtisRule"`, `"GaussKronrodRule"` but none seem to compute it any quicker than just the default. They all agree up to a reasonable number of dp, so no idea if I should just stick to the default, or if there is something better one could do with such an integrand. Speed is not a concern just accuracy. **UPDATE** Let's say I managed to split my integral into a few different integrals. First give the definitions:               vbar[tau_?     NumericQ] := (4 M) ((tau/tauh)^(1/3) + 1) Exp[-(tau/tauh)^(1/3) +       1/2 (tau/tauh)^(2/3) - 1/3 (tau/tauh)];     ubar[tau_?     NumericQ] := -(4 M) ((tau/tauh)^(1/3) - 1) Exp[(tau/tauh)^(1/3) +       1/2 (tau/tauh)^(2/3) + 1/3 (tau/tauh)];     rtau[tau_?NumericQ] := (2 M) (tau/tauh)^(2/3);      in addition to those made above, then I think I can give my integral as a sum of integrands that look like this               Exp[-I s] (ubar[tau_f - s])^(-i 4/10)Exp[+i 1/10 rstar[tau_f - s]]osc1[rtau[tauf-s]]*      here tau_f constant. The first part is an amplitude, the osc1 satisfies the linear ODE given above. I think this has Levin potential if I can work out how to input the LevinRules given the above second order ODE? (Here and in the above I fix my parameters the ODE depends on to (1/10,1) to simplify giving the ICs but I don't that detracts from the main problem). Would need to work out what the Kernal is from the ODE above.


title: Clearing distributed definitions from remote kernels

text: The documentation shows the following example to "Remove a distributed definition..." ![Example](http://i.stack.imgur.com/MQWkv.png) Oh, course this doesn't work without redefining the function prior to re distributing the definition. The following does work:               ClearAll[f];     f[x_] := Labeled[Framed[x^2], $KernelID]     DistributeDefinitions[f];     ParallelMap[f, {1, 2, 3, 4}]      and will give the same output as the above screen clip. OK, but... This seems a very cumbersome way to do this and just doesn't seem right at all. It would seem clearer to me if one could call back the distribution of a function definition, something like `UnDistributeDefinitions[]`, to withdrawal the function from remote kernels rather than needing to clear the function in the main kernel, redefine it, then redistribute it. I came across this issue as I looked to move a DistributeDefinitions[] into a function that will run a bunch of other functions from a standalone cell in a notebook. To test whether it would work, I need to clear all the functions I wanted to distribute or restart the kernel. Just seems a messy way to do things. So, does a way exist to clear definitions in remote kernels without clearing them in the main kernel?


title: Selecting subsets of a list containing a specific element

text: I'm trying to generate all the subsets of a set containing a specific element. For example, for the set {1, 2, 3, 4}, with the specific element I want inside each subset being 4, the required answer would be               {{4}, {1, 4}, {2, 4}, {3, 4}, {1, 2, 4}, {1, 3, 4}, {2, 3, 4}, {1, 2, 3, 4}}      My current code is to do `Select[Subsets[Range[n]], MemberQ[#, n] &]`, where for the above example we would have `n=4`. However, the code seems inefficient and my computer runs out memory for `n>25`. Is there a more efficient way to do so, or is this already the best? NB: I will only need to work with one of the subsets at any point of time, but the `Select` code above doesn't seem to allow me to pick an arbitrary $i^{th}$ subset.


title: Why does compiling a function with ConstantArray give an error when used in parallel?

text: Here is a sample compiled function:               f := Compile[{x}, ConstantArray[0., 10]; x^2, RuntimeAttributes -> {Listable},             Parallelization->True,  CompilationTarget -> "C"]      Using with a single argument works fine:               f[10.0] = 100.0;      But using it in parallel gives an error message:               f[{10.0,20.0}]      > `CompiledFunction::pext:` Instruction 1 in «original function» calls > ordinary _Mathematica_ code that can be evaluated on only one thread at a > time. >> The problem is in the `ConstantArray[0.,10]`, which is not even being used. Replacing it with `Table[0.,{10}]` fixes the problem. Any ideas?


title: Mathematica Error Message using CylindricalPlot3D

text: I have a full version of _Mathematica Home Edition 9_ installed on a high end desktop. I'm trying to learn how to use it and keep coming up with an error message that I cannot figure out. I want to plot a `3D` graph in cylindrical coordinates. I thoroughly know calculus `3`, so the math isn't the problem. The problem is figuring out what I'm doing wrong in _Mathematica_. I searched many websites and found that I need to load a graphics package called `ParametricPlot3D` in order to get the graphics routines necessary. I copy and pasted the exact commands, and I still get the error message. Below is an exact copy & paste of what I typed and what error I got from _Mathematica_. Can anyone please guide me to correct this?               Needs["Graphics`ParametricPlot3D`"]          CylindricalPlot3D[ 9 - r^2, {r, 0, 3}, {theta, 0, 2*Pi}]      > >     Get::noopen: Cannot open Graphics`ParametricPlot3D. >   >     Needs::nocont: "Context \!\(\"Graphics`ParametricPlot3D`\"\) was not > created when >     Needs was evaluated. " >   I have also tried a different way to load the graphics routines, and here's the message I'm getting:               << Graphics``ParametricPlot3D`          CylindricalPlot3D[ 9 - r^2, {r, 0, 3}, {theta, 0, 2*Pi}]      > >     Get::noopen: "Cannot open \!\(\"Graphics`ParametricPlot3D`\"\). " >   I realize that the `character was missing from the ends of the`Graphics` statement, but this text editor won't let me put one at the end and keep that same character at the beginning. Trust me, it is there. I suspect it might be some sort of setting I have to enable, or maybe a directory path that is incorrect, or even I need to download graphics routines. If someone can point me in the right direction, I'd be truly grateful.


title: ContourPlot of a complicated function

text: I have the function $pot_1(x,y,z)$ which is given in analytical form. So, in order to plot the contours at the primary planes $(x,y)$, $(x,z)$ and $(y,z)$, I use the following `Mathematica` code:               Clear["Global`"];     Off[General::spell];     pot1 = -Mn/Sqrt[x^2 + y^2 + z^2 + cn^2];     V = pot1;     Vxy = V /. {z -> 0};     Vxz = V /. {y -> 0};     Vyz = V /. {x -> 0};     Mn=0.08;cn=0.25;     S1 = ContourPlot[Vxy, {x, -10, 10}, {y, -10, 10},      FrameLabel -> {"x", "y"}, RotateLabel -> False, Contours -> 20,      ContourStyle -> Black, ContourShading -> False, PlotRange -> All]     S2 = ContourPlot[Vxz, {x, -10, 10}, {z, -10, 10},      FrameLabel -> {"x", "z"}, RotateLabel -> False, Contours -> 20,      ContourStyle -> Black, ContourShading -> False, PlotRange -> All]     S3 = ContourPlot[Vyz, {y, -10, 10}, {z, -10, 10},      FrameLabel -> {"y", "z"}, RotateLabel -> False, Contours -> 20,      ContourStyle -> Black, ContourShading -> False, PlotRange -> All]      OK, so far so good. Now, I have one more function $pot_2(x,y,z)$ which is much more complicated than $pot_1$ and we can obtain its value only on specific points $(x_0,y_0,z_0)$. This function is given by               pot2[G_, Mt_, a_, b_, c_, x_, y_, z_] := Module[{},     m[u_] := Sqrt[x^2/(a^2 + u) + y^2/(b^2 + u) + z^2/(c^2 + u)];     λ = If[m[0] > 1, u /. FindRoot[m[u]^2 == 1, {u, 1}], 0];     Δ[u_] := Sqrt[(a^2 + u) (b^2 + u) (c^2 + u)];     w000 = NIntegrate[1/Δ[u], {u, λ, \[Infinity]}];     w100 = NIntegrate[1/Δ[u] 1/(a^2 + u), {u, λ, \[Infinity]}];     w001 = NIntegrate[1/Δ[u] 1/(c^2 + u), {u, λ, \[Infinity]}];     w010 = 2/Δ[λ] - w100 - w001;     w110 = (w010 - w100)/(a^2 - b^2); w011 = (w001 - w010)/(b^2 - c^2);      w101 = (w100 - w001)/(c^2 - a^2);     w200 = 1/3 (2/(Δ[λ] (a^2 + λ)) - w110 - w101);     w020 = 1/3 (2/(Δ[λ] (b^2 + λ)) - w011 - w110);     w002 = 1/3 (2/(Δ[λ] (c^2 + λ)) - w101 - w011);     w111 = (w110 - w011)/(c^2 - a^2);     w120 = (w020 - w110)/(a^2 - b^2); w012 = (w002 - w011)/(b^2 - c^2);      w201 = (w200 - w101)/(c^2 - a^2);     w210 = (w110 - w200)/(a^2 - b^2); w021 = (w011 - w020)/(b^2 - c^2);      w102 = (w101 - w002)/(c^2 - a^2);     w300 = 1/5 (2/(Δ[λ] (a^2 + λ)^2) - w210 - w201);     w030 = 1/5 (2/(Δ[λ] (b^2 + λ)^2) - w021 - w120);     w003 = 1/5 (2/(Δ[λ] (c^2 + λ)^2) - w102 - w012);     cc = 15/16 G Mt;     pot = -(cc/6) (w000 - 6 x^2 y^2 z^2 w111 +       x^2 (x^2 (3 w200 - x^2 w300) +           3 (y^2 (2 w110 - y^2 w120 - x^2 w210) - w100)) +        y^2 (y^2 (3 w020 - y^2 w030) +           3 (z^2 (2 w011 - z^2 w012 - y^2 w021) - w010)) +        z^2 (z^2 (3 w002 - z^2 w003) +           3 (x^2 (2 w101 - x^2 w201 - z^2 w102) - w001)));     Print["pot2(x,y,z) = ", pot]     ]      When $G=1$, $Mt=0.1$, $a=7$, $b=1.5$, $c=0.6$ and $-10 \leq x,y,z \leq 10$ I want to plot the contours on the three primary planes of the function $V(x,y,z) = pot1 + pot2$. Any suggestions?


title: Convert a string into a list of numbers

text: I got the following problem I don't know how to solve or if it's at all possible to solve(I'm quite new to _Mathematica_ and don't have an IT background, so please explain slowly if possible:)) I'm trying to convert a simple string which I imported from a homepage into a list of numbers. The string looks for example the following:               {"1,2,3,5,10,12,13,17,26,30,32,41,42,43,113,115,121,125"}      _Mathematica_ sees this as one string, this is important, not as a list of strings. (I used the head function to check that) My question: is it possible to convert this one string into a list of numbers so that each number in the string above is recognized as an individual element and can be calculated with? Another important note: these strings don't all contain the same amount of numbers, a couple of them are shorter than other ones. Is there any way I can convert these different strings into lists of numbers?


title: Color coding based on 4D data within a building

text: I am trying to create something similar to a Heatmap within a building. Suppose I have a room with dimensions, {30,40,50} as length, width, and height respectively. I have a set of data that is of the form {x,y,z,temp} where x,y,z represent the coordinates of the room along the axis of the room, and temp is the temperature at a particular point. My idea is to create a color coding scale like this. If the temp is between 25-30 degrees, then the color at the cordinates within the building is blue, if the temp is between 30-35 degrees, then the color is green, if the temp is between 35-40, the color at that point is red. If there is no temp data available, that particular point should be just a transparent block. I need to arrange cuboids in rows that are predetermined like say, 30 cubes along the length of the building, 40 cubes along the width of the room, and 50 cubes along the height of the building. More correctly, I dont the shapes to be points or cuboids of equal sides but rather readjusting based on the dimensions of the building that I can rearrange along the building. I am trying to achieve something like this, but to readjust itself along the dimensions of the building. I am able to generate the point codes, cuboid data, hues. But I want more precise color coding like the way I have mentioned in the question. Edit:                w = 30;      l = 40;      h = 50;      m = 40;      backwall = {{0, l, 0}, {w, l, 0}, {w, l, h}, {0, l, h}};      side1 = {{0, 0, 0}, {0, 0, h}, {0, l, h}, {0, l, 0}};      side2 = {{w, 0, 0}, {w, 0, h}, {w, l, h}, {w, l, 0}};      floor = {{0, 0, 0}, {w, 0, 0}, {w, l, 0}, {0, l, 0}};      top = {{0, 0, h}, {w, 0, h}, {w, l, h}, {0, l, h}};      front = {{0, 0, 0}, {w, 0, 0}, {w, 0, h}, {0, 0, h}};      leftRoof = {{0, 0, h}, {w/2, 0, (h + m)}, {w/2, l, (h + m)}, {0, l, h}};      rightRoof = {{w, 0, h}, {w/2, 0, (h + m)}, {w/2, l, (h + m)}, {w, l,h}};      roofBack = {{w, l, h}, {w/2, l, (h + m)}, {0, l, h}};      roofFront = {{w, 0, h}, {w/2, 0, (h + m)}, {0, 0, h}};      data = RandomReal[{0, 60}, {10000, 3}];      binning1 = {0, 50, 2};      binning2 = {0, 60, 2};      binning3 = {0, 40, 2};      binned = BinCounts[data, binning1, binning2, binning3];      dims = Dimensions@binned;      coordswithdataAll = Table[{x, y, z, x + y + z}, {x, 1, dims[[1]]}, {y, 1, dims[[2]]}, {z, 1, dims[[3]]}];      coordswithdata = coordswithdataAll;      myTemperatureMap[Null] := {Opacity[0]};      myTemperatureMap[t_] := {Opacity[.2],ColorData["TemperatureMap"][Rescale[t, {25, 40}]]}      temps = Flatten[coordswithdata, 2][[All, 4]];      coords = Flatten[coordswithdata, 2][[All, 1 ;; 3]];      cubes = Graphics3D[{Flatten@MapThread[{#1, Cuboid[#2]} &, {myTemperatureMap /@ temps, coords}]}];      fig = Graphics3D[{Opacity[0.3],Style[Polygon[{backwall, side1, side2, floor, front, leftRoof, rightRoof, roofBack, roofFront}],Lighting -> {{"Ambient", Green}}]}, AxesLabel -> {x, y, z}, Boxed -> False, RotationAction -> "Clip"];      Show[fig, cubes]       What I want to do is to is to arrange the cuboids in a very nice along the length,width, breadth of the warehouse to provide for a visualization for the heat maps within the warehouse. Any pointers, thanks and regards,Anu


title: Generating continuous products or sum

text: This is kind of a continuation of this question. I want to automate these 2 steps on _Mathematica_ once I give an integer $N$.   * generate a list of variables $\{a_i\}_{i=1}^{i=N}$   * then for some function like $f(x) = e^{x}$ or $f(x) = \sin^2(x)$ compute an expression like $\prod\limits_{i \neq j} f(a_i - a_j)$ or $\sum\limits_{i < j} f(a_i - a_j)$ or $\prod\limits_ {i,j} f(a_i - a_j)$ ( - and then hopefully integrate such expressions over all $a_i$...) I am having to type these expressions by hand for every $N$ and that is very hard once $N$ gets large (one would typically have $\binom{N}{2}$ terms to type by hand!) I would like to know how this can be automated - since I would typically need to to use large $N$. (...also, help with that previous question would be great!)


title: Using previous two solutions as a start point in FindRoot

text: I need to solve a sequence of equations and I want the starting point of each new FindRoot to be the solution of the earlier FindRoot, increased or decreased by a a bit in the following way: Some data               d={..,.,...}      Step1: Solve for counter c =100,                {x100,y100}={x,y}/.FindRoot[{f1(x,y,c,d[[c]])==0,f2(x,y,c,d[[c]])==0},{x,a100},{y,b100}]      Step2: Solve for counter c =99                {x99,y99}={x,y}/.FindRoot[{f1(x,y,c,d[[c]])==0,f2(x,y,c,d[[c]])==0},{x,a99},{y,b99}]      The starting points $(a100,b100)$ and $(a99,b99)$ are something I find by trial and error. Step3: For counter c=98,               a98 = x99 + x99-x100;     b98 = y99 + y99-y100;      I know that my roots are monotonic in $c$, therefore I want my new start point to be the previous solution perturbed by how much more or less the previous solution was relative to the solution before that. Finally when I have all my solutions, I would like to store them as a list so I could do more math with them. I already know how to use Table and feed in the previous solution as the new start point.                d = {0.5, 0.4, 0.5, 0.4, 0.5};     sol = {0.9, 0.4};     Table[sol={x,y}/.FindRoot[{x^2+y^2==a,x-y==d[[a]]},{x,sol[[1]]},{y,sol[[2]]}], {a,1,5}]      But I don't know enough about Mathematica Loop commands to figure out how to use the previous solution + the difference between the previous solution and the one before that, as the start point. Also, I don't know when the counter should stop. I'd like for the counter to stop when the solutions first hit a particular region. For example: Stop when $-4 <x < y < 0.1$ or if counter $c = 0$. So basically, I'd like my loop stopping condition to be some function of the counter and/or all the roots calculated so far. Suppose more complicated stopping conditions are hard to figure out, I'd appreciate it if someone to tell me how run a loop with simple stopping conditions like $x<0.1$ unless ofcourse $c=0$. _**EDIT(New stuff):Below is an implementation of Kguler's solution method_** Ok, following Kgluer's comments, I have tried to mimic what he has suggested and it is sort of working but I am not entirely sure of what's going on and I wanna be absolutely certain that I am error free before I implement it on my actual problem. Ok here goes:               f = Table[{x + y == a*a, x - y == 0}, {a, 1, 5}]          ContourPlot[Evaluate[f], {x, 0, 25}, {y, 0, 25}]      ![enter image description here](http://i.stack.imgur.com/SZhhW.jpg) I want to give my FindRoot a start point of $(0,0)$ and a start distance of $(0.5,0.5)$ to look for the first solution when $a=1$. Then when it computes the solution as $(0.5,0.5)$, I want it to look for the next solution, for $a=2$, near the first solution $(0.5,0.5)$+the difference between the first solution and the previous start point$(0.5,0.5)-(0,0)$. So it should look for solutions near $(1,1)$ when $a=2$ and so on and so forth. The parameter $d$ below is superfluous for this problem but I just threw it in there because in my actual problem, $d$ like data shows up non-trivially.               d = {0, 0, 0, 0, 0};     sol = {0, 0};     add = {0.5, 0.5}     counter = 5;          data = NestWhileList[{#[[1]] + 1,{x, y} = {x, y} /.       FindRoot[{x + y == #[[1]]*#[[1]], x - y == d[[#[[1]]]]}, {x, #[[2]][[1]] + #[[3]]  [[1]]}, {y, \#[[2]][[2]] + #[[3]][[2]]}],     {x - #[[2]][[1]], y - #[[2]][[2]]}} &,{1, sol, add},Not[2 < #[[2]][[2]] < #[[2]][[1]] < 3] &, 2, counter]      This code gives the following output.               {{1, {0, 0}, {0.5, 0.5}}, {2, {0.5, 0.5}, {0.5, 0.5}}, {3, {2., 2.}, {1.5, 1.5}}, {4, {4.5, 4.5}, {2.5, 2.5}}, {5, {8., 8.}, {3.5, 3.5}}, {6, {12.5, 12.5}, {4.5, 4.5}}}      So it looks like it is doing what I want it to do.. but I still don't fully understand many thing. Like what would change if I altered the 4th input of NestWhileList from 2 to 3? Or how NestWhileList test function knows to look at the counter. I am sure there are more things I don't understand so any help with writing good code and explanations will be greatly appreciated. Thanks.


title: Is there a convenient/automatic way of uploading a notebook into SE?

text: This post discusses how to streamline the other way - extracting code from a fragmented and longish SE answer but I imagine a two-way workflow like: Fashion an answer/question in an SE-templated notebook with a submit button that converts both text/input cells into their SE equivalents. The notebook could also be effectively uploaded to SE by having the conversion include a link back to a public repository to which the original notebook is also sent.


title: How can I plot the complex graph of $x^x$ in Mathematica?

text: For example, this Wolfram Alpha query shows this graph: ![enter image description here](http://i.stack.imgur.com/zFEXx.png) But it does not show the code for plotting it in _Mathematica_. `Plot[x^x, {x, -1, 1}]` only plots the real values. How can I do this in _Mathematica_?


title: Numerical solution of two coupled transcendental equations with three variables

text: I have following set of coupled transcendental equations. $\tan^{-1} \left(\frac{\sigma+\beta}{K'-1}\right)+\tan^{-1} \left(\frac{\sigma-\beta}{K'-1}\right)=\frac{2 \sigma}{K'}\cos(\beta\tau)$   and   $\log\left(\frac{(K'-1)^2+(\sigma-\beta)^2}{(K'-1)^2+(\sigma+\beta)^2}\right)=\frac{4 \sigma}{K'}\sin(\beta\tau)$ Here K', $\sigma$ and $\beta$ are variables and $\tau=0.07$ is a constant. I want to obtain $\sigma$ vs K' curve (as shown in the attached figure) by eliminating $\beta$.   I tried to eliminate $\beta$ by using _Eliminate_ function of _Mathematica_ , but it gives:   _Eliminate::ifun: Inverse functions are being used by Eliminate, so some solutions may not be found; use Reduce for complete solution information._   When a tried using _Reduce_ , _Mathematica_ keeps on running and does not give anything.   I also tried using FindRoot and NSolve to find values of K' and $\sigma$ that satisfy these equations for a range of $\beta$ going from $-2\pi/\tau$ to $2\pi/\tau$, but I keep on get either $\sigma\approx 0$ or very high.                 Do[Print[FindRoot[{ArcTan[(sigma + beta)/(k - 1)] +        ArcTan[(sigma - beta)/(k - 1)] == ((2 sigma) Cos[0.07  beta])/k,      Log[((k - 1)^2 + (sigma - beta)^2)/((k - 1)^2 + (sigma +           beta)^2)] == ((4 sigma) Sin[0.07  beta])/k}, {{sigma,       2.5}, {k, 2.0}}]], {beta, -2 Pi/0.07, 2 Pi/0.07, 0.01}]      I am relatively new to Mathematica, so please suggest me ways to handle such equations.![enter image description here](http://i.stack.imgur.com/5Cq5j.jpg)


title: Front End Find feature broken in Version 10 with OSX Mavericks?

text: I installed Mathematica 10 and on my Macbook-Pro with OSX Mavricks. I was most frustrated to find that the front-end Find feature doesn't work. Do others have this problem? Is it possible something got messed up during download & installation?


title: Image Processing: Detecting and counting small objects

text: I would like to analyze close-up images of human hair. My goal is to detect and count the number of hair strands and determine their width, but I am having a hard time binarizing the image. Seems like there is too much noise. Any suggestions ? Image is below: ![enter image description here](http://i.stack.imgur.com/Oq8j2.jpg) I used the following code and got a result that does not detect all the hairs:               arm2 = ColorNegate[DeleteSmallComponents[Binarize[arm, 0.42]]]      ![enter image description here](http://i.stack.imgur.com/NYMyW.jpg)


title: Split absolute value in chambers

text: I usually expand expression containing absolute values, with the PiecewiseExpand[] function.               Delta[m1_, m2_, m3_] := -Abs[m2 - m3] + 1/2*(Abs[m1 - m2] + Abs[m1 - m3]) + 3/2 (Abs[m2] + Abs[m3])          PiecewiseExpand[Delta[m1, m2, m3], {m1 ∈ Integers, m2 ∈ Integers, m3 ∈ Integers, m2 >= m3}] // FullSimplify      This usually splits the absolute value in the various chambers in which it assumes a certain functional form which does not contain absolute values anymore, and also the conditions which define the cells are given by the PiecewiseExpand[] function. However, the last chamber is never given, but it corresponds to the bool vale "True". How can I recover it?


title: How can AbsoluteTimes be split by date faster?

text: I have data where the date is stored as an `AbsoluteTime`. Here's a dummy dataset:               sampleData=Transpose[{#,Range[Length[#]]}]&[Range[AbsoluteTime[{2012,10,01}],         AbsoluteTime[{2013,09,30}],60]];      I want to split it by date, so I use:               splitData=SplitBy[sampleData,DateList[First[#]][[;;3]]&];      This takes ~22 s on my machine. Often I have bigger data sets than this. Is there a faster way to get `splitData`?


title: How to speed up multi-curve plotting?

text: The following example is a very simplified version of the problem I'm working with, though I hope it still captures the essentials. First, the definitions:               sigmoid[s_, a_][x_] := (a (-s + x))/Sqrt[1 + a^2 (s - x)^2];          plots[params_] := Module[{fns},        fns = sigmoid @@ # & /@ params;        Plot[#[x] & /@ fns, {x, -20, 20}, ImageSize -> Small]];      The `plots` function produces a plot containing as many sigmoid curves as there are parameter sets in the list `params`. E.g., with 5 parameter sets (differing only in the first parameter):               plots[{{-10, 1}, {-5, 1}, {0, 1}, {5, 1}, {10, 1}}]      `plots` produces ![5-curve plot](http://i.stack.imgur.com/kTwgA.png) More typically, however, I need to run `plots` with an argument consisting of ~100 sets of parameters, and I'm finding it to take considerably longer than I had expected. (Maybe I've gotten spoiled!) For example, running               plots[{#, 1} & /@ Array[# &, 100, {-15, 15}]]      takes ~40s on my desktop. Does this sound right? Currently I am in the process of producing several figures each consisting of a grid containing ~100 plots of comparable complexity (i.e. ~100 curves/plot), so I need to find ways to significantly speed this up. Any suggestions would be appreciated. NB: The details of the `sigmoid` function should not be an important aspect of the proposed optimization, since I chose this function for this example only because it is significantly simpler than the function that I am actually working with. (For one thing, the latter is parametrized by three rather than two numbers. Nevertheless, it produces s-shaped curves looking very similar those produced by `sigmoid`.)


title: Why can NSolve only solve this system of equations after trivial substitution?

text: The first system immediately yields {}, but the second, equivalent, system yields the solution `{{k -> 4.96048, tl1 -> 1.64842, tl2 -> 2.07466}}` :               NSolve[{       (5 - k)(2940/100) == tl1^0.3,        (5 - k)(3150/100) == tl2^0.3,        4200 tl1^0.7 + 4500 tl2^0.7 == 1000000},        {k, tl1, tl2}]          NSolve[{       (5-k)(2940/100) == tl1^0.3,       (5-k)(3150/100) == tl2^0.3,        (5-k) 4200(2940/100)^(7/3) + (5-k) 4500 (3150/100)^(7/3) == 1000000},        {k, tl1, tl2}]      The second system is created by noting that `tl1^0.7 = (tl1^0.3)^(7/3)` and `tl2^0.7 = (tl2^0.3)^(7/3)` and substituting the results into the third equation. Is there a way to get _Mathematica_ to do such substitutions on its own?


title: Understanding a slowdown associated with wrapping a set of nested Do loops with a ParallelDo expression

text: I have a set of set of nested Do loops that require data from three different arrays (though the position from each array is fixed per cycle of the outer Do loop). Let's call these arrays: ListA, ListB, ListC. My attempts at parallelizing the outer Do loop (i.e. using ParallelDo) seem to work, but there is a ~100x slowdown on a 12 core processor. I was hoping someone here could help me understand why. Here's some nonsense pseudocode that should capture the structure of the bit of code I evaluated:               Do[        Do[           Do[              ListA[[i]] = Append[ListA[[i]], ListB[[i, a]] + ListB[[i, b]]];             , {b, (a+1), M}             ];          , {a, 1, M}          ];               ListA[[i]] = Sort[ListA[[i]]];                  Do[             If[Abs[ListA[[i, j]] - ListA[[i, j + 1]]] < ListC[[i]],                  ListC[[i]] = Abs[ListB[[i, j]] - ListB[[i, j + 1]]];               ];           , {j, 1, Length[PairwiseSums[[i]]] - 1}          ];               , {i, 1, Length[AllSubsets]}     ];      Previously defining ListA, ListB, and ListC, and writing:               SetSharedVariable[ListA]     SetSharedVariable[ListB]     SetSharedVariable[ListC]          ParallelDo[        Do[           Do[              ListA[[i]] = Append[ListA[[i]], ListB[[i, a]] + ListB[[i, b]]];             , {b, (a+1), M}             ];          , {a, 1, M}          ];               ListA[[i]] = Sort[ListA[[i]]];                  Do[             If[Abs[ListA[[i, j]] - ListA[[i, j + 1]]] < ListC[[i]],                  ListC[[i]] = Abs[ListB[[i, j]] - ListB[[i, j + 1]]];               ];           , {j, 1, Length[PairwiseSums[[i]]] - 1}          ];               , {i, 1, Length[AllSubsets]}     ];      Yields the ~100x slowdown on a 12 core processor. Why?


title: how to efficiently jump between sections in a notebook?

text: I have functions defined in some section of the notebook and some cells where I use these functions maaaaany pages below in the same notebook. If I want to go and change the definition of the function I need to climb up in the notebook, find the definition and change it, then climb down to the place where I was working and using these functions. This is very inconvenient and looks like "navigation by sight". I was wondering how this can be made more efficient. For instance if it was a C++ or TeX source I would use the table of content of the file that my editor builds for me so that I can jump to the section where I define function clicking the TOC and then jump back to the original location either using "back" or again clicking the TOC. Can anything like this be done in mathematica?


title: Better method for creating a tuple out of two lists

text: I just combined two lists using this method               molefrac = {1., 0.787402, 0.478927, 0.282592, 0.}     calibrationIndexOfRefraction = {1.440, 1.421, 1.395, 1.380, 1.351}     Partition[Riffle[molefrac, calibrationIndexOfRefraction], 2]      > >     {{1., 1.44}, {0.787402, 1.421}, {0.478927, 1.395}, {0.282592, 1.38}, > {0., 1.351}} >   The purpose of this was to prepare the lists to be plotted using `ListPlot`. Is there some way to use `ListPlot` directly (or something more functional) to do this without going through all the extra steps?


title: Efficiently finding the maximum value of a column in a matrix

text: Imagine I have a matrix of the following form: TestArray = $((a_1,b_1,...),(a_2,b_2,...),(a_3,b_3,...),...)$, where a specific example might be: TestArray = $((10^6,1,...),(10^7,7,...),(10^3,4,...))$ In the above example, we have $Max[a_i] = 10^7$, and $Max[b_i] = 7$. In a single line of Mathematica, how can I find the maximum $a_i$, $b_i$, etc.? Moreover, what is the most computationally efficient manner in which to do this?


title: Undirected Network and its Monte Carlo Simulation

text: So I have this array of points, in the form {x, {y,z}}, where x is the node number assigned when I used the locator pane to make a network, and y and z are its unitary coordinates. The graph is made up of two distinct node types, let's call them A and B. I want to be able to distinguish these two by assigned different shapes to their nodes and also want to connect the network how it originally was and using Primm's Algorithm (Minimum Spanning Tree). I was wondering how I go about it. I should also probably mention that from this network I am going to run some Monte Carlo Simulations on it. The code I have so far is:               graphicslist =        Table[Rectangle[mpg[[i, 2]] , mpg[[i, 2]] + {.005, .005}], {i, 1, Length@mpg}];       Graphics[graphicslist]       where mpg is the array of my nodes and their coordinates. This returns: ![enter image description here](http://i.stack.imgur.com/juI6d.png) Small sample of the mpg:               mpg == ({{1, {0.0996, 0.4195}}, {2, {0.1046, 0.392}}, {3, {0.1144,       0.468}}, {4, {0.1164, 0.4415}}, {5, {0.1432,       0.4455}}, {6, {0.1442, 0.376}}, {7, {0.147, 0.5085}}, {8, {0.16,       0.4965}}, {9, {0.1628, 0.5305}}, {10, {0.1708, 0.484}}.......      Here the integer is the node number and the following decimals are the unitary coordinates of its corresponding vertex. In the real Maryland Power Grid, there are two distinct vertices; loads and generators. I plotted, using the locator pane, all of its vertices and wish to differentiate which one are loads and which are generators rather than just have squares for both of them.


title: How can I extract when a select condition is met in a do loop the value of the iterator?

text: I use the following code to find out the number of consecutive prime numbers using a formula $n^2+n+i$ found out by Euler (starting from n=0):               Nbs = {};     Do[Nbs = Union[Nbs,      Select[Range[5000], (PrimeQ[#^2 + # + i] == False &), 1]], {i, 1,      5000}];     Nbs       How can I also get in the output list the value of the `Do` iterator where $i$ is corresponding to each number of consecutive primes? I would like to get something like that:               ({1,2},...,{40,41}}


title: symbolic integration of a definite integral, how to speed up

text: I am wondering what are the general tips & hints for speeding up symbolic integration. As far as I understand, it solely depends on the form a function is represented before it is processed by `Integrate`. So, `Simplify`, `ExpandAll`, `PowerExpand` _etc._ should be used. Are there any other strategies? My monster function (it looks relatively simple, though) is:               f[r_]:=(A^(2) E^((2 r α)/(3)) α^(3) (A^(8/3) \     (-0.002396-0.001198 r α)+E^((8 r α)/(3)) α^(4) \     (0.00001755-(8.776*10^(-6)) r α)+A^(4/3) E^((4 r \     α)/(3)) α^(2) (-0.0002734+0.0004785 r α)))/(r \     (1.000 A^(4/3)+0.02853 E^((4 r α)/(3)) α^(2))^(4))      or (TeX) $$ \frac{A^2 \alpha ^3 e^{\frac{2 \alpha r}{3}} \left(\alpha ^4 \ e^{\frac{8 \alpha r}{3}} \left(0.00001755-\alpha r \left(8.776 \times 10^{-6}\right)\right)+\alpha ^2 A^{4/3} (0.0004785 \alpha \ r-0.0002734) e^{\frac{4 \alpha r}{3}}+A^{8/3} (-0.001198 \alpha \ r-0.002396)\right)}{r \left(0.02853 \alpha ^2 e^{\frac{4 \alpha \ r}{3}}+1.000 A^{4/3}\right)^4} $$ where $\alpha > 0$ and $A > 0$ are some parameters. What I am trying to do: `Integrate[f[r]*r^2*4*Pi,{r,O,Infinity}]` Sure, I could have used `NIntegrate` for some specific $\alpha$ and $A$ (and belive me, I did it:), but now I am interested in the general result. Seems like the situation is hopeless, because my laptop can not handle it.


title: Plotting discrete data but not using discreteplot function

text: For a given function:               Plot[Sqrt[Abs[x]], {x, -Pi, Pi}]      I have the code to draw the function (with its Abs remove), partial sums and cesaro means as:               f[x_] := Sqrt[x]     s[k_, x_] := \frac{2\sqrt{\pi}}{3}+(-Sqrt[2] FresnelS[Sqrt[2] Sqrt[n]] + 2 Sqrt[n] Sin[n \[Pi]])/(n^(3/2) Sqrt[\[Pi]]) Cos[n x], {n, 1, k}]      partialsums[x_] = Table[s[n, x], {n, {4}}];      c[n_, x_] := (1/n) Sum[s[m, x], {m, 0, n - 1}]     Plot[Evaluate[{f[x], partialsums[x], c[4, x]}], {x, -Pi, Pi},     PlotLegends -> {"f(x)=x", "Fourier, 4 terms", "Cesaro, 4 terms"},      PlotStyle -> {{Blue}, {Dashed, Thickness[0.006]}, {Red, Thickness[0.006]}}]      This code fails on my computer and hence I resolve to manual computation. Updates: It turn out that I can easily solve this issue by removing the k with any number rather than letting it to be indefinite. Although I am not certain the graph is right for k=4 as both graphs(Partial and Cesaro) coincides with each other.


title: Rookie mistake in defining a function with Modules?

text: I can't figure out what's wrong with this piece of code. I'm trying to define a function that counts the occurrences of `x`, `y` and `l` or `lc` in expressions like `x[x[y[l]]`, where `x` will add 1 to an "eigenvalue", `y` will subtract 1 and likewise for `l` and `lc`. Then for example I would get `H[x[x[y[l]]] = 2`. Now I'm stuck at the first phase with the function defined below:               H[z_] := Module[{list1 = Join[Level[z, Infinity], {z}], lam = 0},       Which[list1[[1]] == l, lam = lam + 1,              list1[[1]] == lc, lam = lam - 1       ];       lam     ]      This should give `H[l] = 1` and `H[lc] = -1`, but instead I get just the initial value `lam = 0` for the "lc" case... When I discard the last line and use instead               H[z_] := Module[{list1 = Join[Level[z, Infinity], {z}], lam = 0},       Which[list1[[1]] == l, lam = lam + 1,              list1[[1]] == lc, lam = lam - 1       ]     ]      I get               Which[lc == l, lam$538474 = lam$538474 + 1,            list1$538474[[1]] == lc, lam$538474 = lam$538474 - 1     ]      Clearly something's wrong with the `list1$538474[[1]]` above. What am I missing?


title: Plot, extract data to a file

text: I need to save data that _Mathematica_ uses inside the `Plot` command. The format that I would like to have is:               x1 y1     x2 y2     .....      i.e. basically, two columns, _x_ and _y_ (sorted, if possible) so that I can use the data for other programs. I already played with `Export` and got the data enclosed in `{}` and tons of other information I don't need: color, axis, _etc._ What is the nicest way to do what I am trying to do?


title: How to extract data from a plot?

text: Here is the equation I want to plot,               t = -((E^-η h1^2 γ)/(16 (1 + γ)^7)) -      (E^-η h1^3 γ)/(32 (1 + γ)^7) - (     E^-η h1^2 Nr γ)/(6 (1 + γ)^7) - (     E^-η h1^3 Nr γ)/(8 (1 + γ)^7) - (     E^-η h1^2 Nr^2 γ)/(9 (1 + γ)^7) - (     E^-η h1^3 Nr^2 γ)/(6 (1 + γ)^7) - (     2 E^-η h1^3 Nr^3 γ)/(27 (1 + γ)^7) - (     4 E^-η h1^2 Pr γ)/(9 (1 + γ)^7) - (     85 E^-η h1^3 Pr γ)/(216 (1 + γ)^7) - (     35 E^-η h h1 M^2 Pr γ)/(72 (1 + γ)^7);               Clear[η, M, γ, h, h1, Nr, Pr];     θ[η_, M_, γ_, h_, h1_, Nr_, Pr_] = t;     Plot[{  θ[η, 1.5, 0.6, -0.5, -0.5, -0.5, 0.1],      θ[η, 1.5, 0.6, -0.5, -0.5, -0.5, 0.1],      θ[η, 1.5, 0.6, -0.5, -0.5, -0.5, 0.3],      θ[η, 1.5, 0.6, -0.5, -0.5, -0.5, 0.5],      θ[η, 1.5, 0.6, -0.5, -0.5, -0.5, 0.7],      θ[η, 1.5, 0.6, -0.5, -0.5, -0.5, 0.9],      θ[η, 1.5, 0.6, -0.5, -0.5, -0.5, 1.1]},      {η, 0, 8}, PlotRange -> All, AxesLabel -> {"\!\(\*StyleBox[\"η\",\nFontSize-   >22]\)", " \!\(\*StyleBox[\"θ\",\nFontSize->20]\)\!\(\*StyleBox[\"(\",\nFontSize->22]\)\!\(\*StyleBox[\"η\",\nFontSize->22]\)\!\(\*StyleBox[\")\",\nFontSize->22]\)"}, PlotStyle -> Directive[Black],AxesStyle -> {Directive[Black, 12, Bold], Directive[Black, 12, Bold]}]      How to save data from this plot as a .txt file? Thanks


title: Manipulation of Automatic settings in a given case

text: Often I use the `Automatic` option to specify plot styles, plot ranges, image sizes, etc. I'd like to be able to apply some modifications to these automatically defined values. For example, if a plot with the `PlotRange -> Automatic` would draw the function on the interval $[0,1]$, I'd like to be able to use something like `PlotRange -> #*2 & /@ Automatic` to extend the plot range to $[0,2]$. Is this possible in general, or maybe every case would need a specific solution? _I'm not sure what tags are most appropriate for this topic, please add them if you think they're more suitable than the one I've chosen._


title: Coloring a honeycomb structure

text: I am trying to generate a honeycomb structure. With the help of this MMA.SE question & this website, I got               hexgon = Polygon[       N@Transpose[Through[{Re, Im}[Exp[2 Pi I Range[6]/6]]]]];               points[x_, y_] :=        Flatten[Table[{{3 m, Sqrt[3] n}, {3 m + 3/2,            Sqrt[3] n + Sqrt[3]/2}}, {m, 0, x}, {n, 0, y}], 2];          Graphics[Translate[       Rotate[{RGBColor[RandomReal[], RandomReal[], RandomReal[]],          hexgon}, \[Pi]/6], ({#[[2]], #[[1]]} &) /@ points[1, 4]]]      ![Honeycomb](http://i.stack.imgur.com/Z9WkI.png) How do I fill every hexagon with different color?