how would you do it?

I have expected that the code generating the polygons were capable to optimize and avoid overlappings.

how would you do it?

I have expected that the code generating the polygons were capable to optimize and avoid overlappings.

I find that ?solve,details is precisely one of the help pages (of common use) that refers more heavily to environment/global variables: _SolutionsMayBeLost, _EnvConditionalSolutions, _EnvTryHard, _EnvExplicit, _MaxSols,  _SolutionsMayBeLost, _EnvAllSolutions, ...

Observing the combination of points and polygons:

with(plots):
with(plottools):
data2:=...
p1:=pointplot3d(data2, color = red):
p2:=polygon(data2):
display([p1,p2],axes=boxed,orientation=[74,60.]);

I think that there are two problems here:
 

1. Some points are not joined by segments, making the borders ragged.

2. Some polygons are not correctly painted (probably related to a longstanding bug), paradoxically making the borders a bit less ragged.

Observing the combination of points and polygons:

with(plots):
with(plottools):
data2:=...
p1:=pointplot3d(data2, color = red):
p2:=polygon(data2):
display([p1,p2],axes=boxed,orientation=[74,60.]);

I think that there are two problems here:
 

1. Some points are not joined by segments, making the borders ragged.

2. Some polygons are not correctly painted (probably related to a longstanding bug), paradoxically making the borders a bit less ragged.

whether Maplesoft has, for internal use, such a reference document listing all the available enviromental variables, with their description and all other details.

for '_EnvMaxCosetsToddCoxeter ' finds it.

I bet that there is no published document that has the complete list of environment variables. It seems an issue fit to add to the Maple Wiki.

In Mathematica 5:

N[Hypergeometric1F1[1+I,2*I,-10^6]]

6.739100053151782*10^(-9) + 2.0757700378175314*10^(-7)*I

is very fast.

In Mathematica 5:

N[Hypergeometric1F1[1+I,2*I,-10^6]]

6.739100053151782*10^(-9) + 2.0757700378175314*10^(-7)*I

is very fast.

ie

KummerM(a,b,z); subs(a=1,b=3/2,z= -10^6,%); evalf(%);
 
it is also slow.

ie

KummerM(a,b,z); subs(a=1,b=3/2,z= -10^6,%); evalf(%);
 
it is also slow.

It is a=1 and b=3/2:

convert(hypergeom([1],[3/2],-z),erfc);
                          1/2                1/2
                        Pi    (-1 + erfc((-z)   ))
                   -1/2 --------------------------
                                  1/2
                              (-z)    exp(z)

It is a=1 and b=3/2:

convert(hypergeom([1],[3/2],-z),erfc);
                          1/2                1/2
                        Pi    (-1 + erfc((-z)   ))
                   -1/2 --------------------------
                                  1/2
                              (-z)    exp(z)

I do not know whether there is an algorithm to find all the roots for eta, taking into account all possible cases with positive parameters. 

As a first thing, I would change to new variable/parameters without squares to avoid duplication of cases for each of them:

eq:= alpha^2*sqrt(1-(alpha*eta/kappa[B])^2)
*((2-eta^2)^2-4*sqrt(1-(eta/kappa)^2)*sqrt(1-eta^2))
+mu[B]*sqrt(1-(eta/kappa)^2)*((2-(alpha*eta)^2)^2
-4*sqrt(1-(alpha*eta/kappa[B])^2)*sqrt(1-(alpha*eta)^2))/mu:
eq1:=subs(eta^2=y,alpha^2=A,kappa[B]=sqrt(B),kappa=sqrt(K),eq);

           /    A y\1/2         2              1/2        1/2
  eq1 := A |1 - ---|    ((2 - y)  - 4 (1 - y/K)    (1 - y)   ) +
           \     B /

                       1/2 /         2     /    A y\1/2          1/2\
        mu[B] (1 - y/K)    |(2 - A y)  - 4 |1 - ---|    (1 - A y)   |
                           \               \     B /                /

        /mu

The positivity constraint will have to be applied to these variable/parameters. Then I think that finding roots expressing the parameters in terms of y is computationally easier than the other way round. So I would start with the solution sets from:

solve(eq1)

and choose values of y>0 that make the parameters positive. Eg. This one seems fine:

{A = A, B = A*y, K = y, mu = mu, y = y, mu[B] = mu[B]}

I do not know whether there is an algorithm to find all the roots for eta, taking into account all possible cases with positive parameters. 

As a first thing, I would change to new variable/parameters without squares to avoid duplication of cases for each of them:

eq:= alpha^2*sqrt(1-(alpha*eta/kappa[B])^2)
*((2-eta^2)^2-4*sqrt(1-(eta/kappa)^2)*sqrt(1-eta^2))
+mu[B]*sqrt(1-(eta/kappa)^2)*((2-(alpha*eta)^2)^2
-4*sqrt(1-(alpha*eta/kappa[B])^2)*sqrt(1-(alpha*eta)^2))/mu:
eq1:=subs(eta^2=y,alpha^2=A,kappa[B]=sqrt(B),kappa=sqrt(K),eq);

           /    A y\1/2         2              1/2        1/2
  eq1 := A |1 - ---|    ((2 - y)  - 4 (1 - y/K)    (1 - y)   ) +
           \     B /

                       1/2 /         2     /    A y\1/2          1/2\
        mu[B] (1 - y/K)    |(2 - A y)  - 4 |1 - ---|    (1 - A y)   |
                           \               \     B /                /

        /mu

The positivity constraint will have to be applied to these variable/parameters. Then I think that finding roots expressing the parameters in terms of y is computationally easier than the other way round. So I would start with the solution sets from:

solve(eq1)

and choose values of y>0 that make the parameters positive. Eg. This one seems fine:

{A = A, B = A*y, K = y, mu = mu, y = y, mu[B] = mu[B]}