@tomleslie It won't work because the parser objects to having a loop (or any statement) split across multiple execution groups.

@acer 

This probably already occurred to you, but you can avoid the overhead of calls to an anonymous procedure v-> F.v by using map2. Change map( v -> F . v, CCXYZ_D50 ) to

map2(`.`, F, CCXYZ_D50 )

@Markiyan Hirnyk Thank you very much for the encouragement. That's a beautiful painting. The clouds are similar to the Gaussian Random Fields with large convolution radii. I find them also beautiful. Here's an example of a GRF(512, 128) (the 128 is the convolution radius). (Sorry for the excessive whitespace---these are 3D plots projected onto one of the coordinate planes; I don't know how to eliminate it.) It looks like this could be a model for how the shapes of continents, islands, and oceans develop.

In addition to the bracket problem, you need to change i. In Maple, sqrt(-1) is represented by default as I. You can change that to i if you want, but it requires an explicit command.

@Markiyan Hirnyk 

Okay, how about this?

restart:

Digits:= 15:

GRF:= proc(n,r)
uses AT= ArrayTools, IT= ImageTools;
local
     GRF:= IT:-Convolution(
          AT:-RandomArray(n, 'distribution'= 'normal'),
          Matrix(2*r+1$2, (i,j)-> `if`((i-r-1)^2 + (j-r-1)^2 <= r^2, 1, 0))        
     ),
     Pos:= map[evalhf](`>=`, GRF, 0),
     Neg:= map[evalhf](`<`, GRF, 0)
;
     (IT:-Preview@IT:-HSVtoRGB)(
          IT:-CombineLayers(.6*360*Pos, IT:-FitIntensity(Pos*~GRF), Pos) +~    
          IT:-CombineLayers(.1*360*Neg, IT:-FitIntensity(-Neg*~GRF), Neg)
     )
end proc:

GRF(200, 13);

The underlying greyscale image is the same as what I posted previously; only the color encoding is different. The two colors can be easily changed to any pair on the Hue (the H of HSV) scale, i.e., colors that correspond to specific wavelengths of light. I could put a gamma encoding on the Saturation (the S of HSV) layer to raise some of the whitish tints. (The normally distributed values close to zero are the colors close to white.) Make sure to view the image from a variety of oblique angles. It can make a big difference because of how the RGB layers are laid out in your monitor screen.

Here's a more-sophisticated color encoding that is fairly close to the Mathematica-generated plot that you posted:

GRF:= proc(n,r)
uses AT= ArrayTools, IT= ImageTools;
local
     GRF:= IT:-Convolution(
          AT:-RandomArray(n, 'distribution'= 'normal'),
          Matrix(2*r+1$2, (i,j)-> `if`((i-r-1)^2 + (j-r-1)^2 <= r^2, 1, 0))        
     ),
     Pos:= map[evalhf](`>=`, GRF, 0),
     Neg:= map[evalhf](`<`, GRF, 0)
;
     (IT:-Preview@IT:-HSVtoRGB)(
          IT:-CombineLayers(
               .6*360*Pos*~IT:-FitIntensity(Pos*~GRF, .85..1.2),
               IT:-FitIntensity(IT:-FitIntensity(Pos*~GRF)^~.75, .08..1),
               .93*Pos
          ) +~    
          IT:-CombineLayers(
               .1*360*Neg*~IT:-FitIntensity(Neg*~GRF, .4..1.5),
               IT:-FitIntensity(IT:-FitIntensity(-Neg*~GRF)^~.75, .08..1),
              .93*Neg
          )
     )
end proc:

GRF(200, 12);

This uses two ranges of Hues (positive (blue): 0.51..0.72, negative (orange): 0.04..0.15; on a 0..1 scale) rather than two specific Hues. A gamma (or exponent) of .75 is applied to the Saturation layer and then the minimal value is set to .08 (on a 0..1 scale). The Value layer is constant .93 (on a 0..1 scale).

Again, this is just a translation (plus color encoding) of the first of the many Matlab programs given in that paper, and this is not the algorithm for the plot that you posted, which, I believe, was done with FFT rather than convolution. I may translate some of the others.

@tomleslie I really like that as a quick and dirty solution, but it's difficult to generalize it to an arbitrary number of matrices.

@Carl Love 

So, how did I get those expressions for the partial derivatives of F, the key to the previous computation? Working with multivariate chain rules and the inverse function theorem to get derivatives by hand is difficult. Here are two ways to do the symbolics with Maple:

F:= (y,z)-> RootOf(MF(x,y,z)-11, x):
D[1](F)(y,z);

D[2](F)(y,z);

The following way is perhaps easier to understand:

restart:
alias(F(y,z)= RootOf(MF(x,y,z)=11, x)):
diff(F(y,z), y);

and likewise for diff(F(y,z), z).

@Haley_VSRC 

I just substantially updated my previous Reply. (I was doing it while you were writing yours.) Please read that. I asked a few more questions. And I gave formulas for the partial derivatives of f in terms of f and the partial derivatives of MF, which I think is specifically what you're asking for. Let me know if that answers your Question. If so, I'll upgrade my Reply to an Answer.

So, I take your P.S. to mean that the isosurface is topologically disk-like for given intervals of y and z and that there's a unique x for each (y,z) in those intervals. In other words, x = f(y,z) is a function (in the strict sense of the word) defined over a rectangle (y0 <= y <= y1) x (z0 <= z <= z1) in the (y,z) plane. Is that correct? That'd probably be the easiest situation for finding the surface area by numeric computation.

I changed your title to reflect what may be possible with Maple. If you don't approve, suggest another, and I'll change it again.

It's an interesting problem. In order to wrap by brain around it, it'd help me to know the topology of the surface. Is it a simple closed surface, i.e., is it topologically homeomorphic to a sphere? I don't know if that ultimately will make a difference in the computation; it may just help me think about it. Being homeomorphic to a disk is the other case that's relatively easy to think about.

Let me verify what you have currently: You currently have a procedure that will give you a numeric approximation of MF(x,y,z) for any x, y, and z (within certain ranges, of course). Is that correct? So we could use numeric differentiation (command fdiff) to approximate all partial derivatives of MF.

You suggested expressing x as a function of y and z. Is there a reason that you chose x, such as that there's a unique x for every (y,z) pair (which would necessarily mean that the topology is not sphere-like)? Given that you have a procedure MF as described in the previous paragraph and given that there's a unique x for each (y,z) pair, then we can write the "x return" function for the isosurface as

F:= proc(y,z) local x; fsolve(MF(x,y,z)=11, x) end proc

using a minor modification to MF to allow it to return unevaluated when passed any symbolic input. If the x isn't unique, then it's trickier. The partial derivatives of F (which you called Fy and Fz) can easily be expressed in terms of F and the partial derivatives of MF. They are 

D[k](F)(y,z) = -(D[k+1]/D[1])(MF)(F(y,z), y, z))  (for k = 1,2)

Let me know if you need help understanding Maple's D notation that I just used.

What do you mean by "special" derivatives of MF?

@Markiyan Hirnyk 

My blues are graded from medium blue to black. The plot that you showed used medium blue to white. Mathematically, it's the same thing, but it's harder to see the gradation blue to black. I'll try to make it blue to white using the HSV scale. It's a slightly more difficult computation.

@Markiyan Hirnyk 

Okay, granted, the color is important. Anyway, that's just what I did. The blues are the positives and the yellow-greens are the negatives. I think that Maple's color scaling may be wrong. Tomorrow I'll try using HSV instead of RGB.

Note that I only implemented the first of the many, many Matlab programs given in the paper, and this isn't the algorithm used for the plot that you showed.

@Markiyan Hirnyk 

I used the Matlab code from page 4 of this paper. I didn't use or even look at the Mathematica code. The mathematical object is represented by a greyscale image. Adding the color is only "art for art's sake" as you like to say, and I only did it because the image that you posted was split into orange and blue. I think that you are not seeing the variation in the blue fields in my plot. Try looking at the underlying greyscale images. Use this code:

restart:

Digits:= 15:
GRF:= proc(n::posint, r::posint)
uses AT= ArrayTools, IT= ImageTools;
     (IT:-View@IT:-FitIntensity@IT:-Convolution)(
          AT:-RandomArray(n, 'distribution'= 'normal'),
          Matrix(2*r+1$2, (i,j)-> `if`((i-r-1)^2 + (j-r-1)^2 <= r^2, 1, 0))        
     )
end proc:

GRF(200, 12);

Compare with the left greyscale image on page 2 of the paper.

All of the numeric information from the greyscale images is included in the color images; it just may be difficult to see variations in some color ranges.

@Preben Alsholm 

Okay, the spam-deleters strike is on.

@acer Thanks for the correction!

@rlopez 

I appreciate where you're coming from with that desk-drawer analogy, but I'm afraid that it might give some readers the wrong idea about sets for some uses unrelated to this post. Sets are stored in a strict, reliable, and predictable order that is determined by Maple---the user has no control over it. Long ago, that order was session dependent (being partially determined by memory addresses, I believe). In modern Maple, the order is consistent across sessions and even across different computers.