Check this:

plot3d options : lightmodel, ambientlight, color

Perhaps lightmodel=light1 does what you're looking for. Play around with this and let us know which options you found best.

Notice how line breaks (SHIFT+RETURN) make it easier to read your code.

And creating a separate "execution group" for the plot allows you to play with your options quicker, in that you don't have to recompute everything each time.

I have inserted extra options that are available by right-clicking, in the way suggested by jakubi above, and which you can also manually force. See if you can set the color to Z (Hue) (I couldn't off the top of my head just now).

> restart:

with(LinearAlgebra):

with(plots):

> N := 10:

n1 := 1.5:

n2 := 3.4:

n3 := 3.5:

d1 := 0.16e-6:

d2 := d1:

d3 := d1:

c := 0.3e9:

theta2 := arcsin(n1*sin(theta1)/n2):

theta3 := arcsin(n2*sin(theta2)/n3):

`#mover(mi("n1"),mi("~"))` := n1*cos(theta1):

`#mover(mi("n2"),mi("~"))` := n2*cos(theta2):

`#mover(mi("n3"),mi("~"))` := n3*cos(theta3):

Lambda := d1+d2+d3:

`#mover(mi("n"),mo("&uminus0;"))` := (n1*d1+n2*d2+n3*d3)/Lambda:

`?b` := c/(`#mover(mi("n"),mo("&uminus0;"))`*(2*Lambda)):

`evalf/constant/?B` := proc () options operator, arrow:

`?b` end proc:

phi1 := 2*Pi*n1*d1*nu*cos(theta1)/c:

phi2 := 2*Pi*n1*d1*nu*cos(theta2)/c:

phi3 := 2*Pi*n1*d1*nu*cos(theta3)/c:

M111 := (`#mover(mi("n2"),mi("~"))`+`#mover(mi("n1"),mi("~"))`)*exp(-I*phi1):

M112 := (`#mover(mi("n2"),mi("~"))`-`#mover(mi("n1"),mi("~"))`)*exp(I*phi1):

M121 := (`#mover(mi("n2"),mi("~"))`-`#mover(mi("n1"),mi("~"))`)*exp(-I*phi1):

M122 := (`#mover(mi("n2"),mi("~"))`+`#mover(mi("n1"),mi("~"))`)*exp(I*phi1):

M211 := (`#mover(mi("n3"),mi("~"))`+`#mover(mi("n2"),mi("~"))`)*exp(-I*phi2):

M212 := (`#mover(mi("n3"),mi("~"))`-`#mover(mi("n2"),mi("~"))`)*exp(I*phi2):

M221 := (`#mover(mi("n3"),mi("~"))`-`#mover(mi("n2"),mi("~"))`)*exp(-I*phi2):

M222 := (`#mover(mi("n3"),mi("~"))`+`#mover(mi("n2"),mi("~"))`)*exp(I*phi2):

M311 := (`#mover(mi("n1"),mi("~"))`+n3)*exp(-I*phi3):

M312 := (`#mover(mi("n1"),mi("~"))`-n3)*exp(I*phi3):

M321 := (`#mover(mi("n1"),mi("~"))`-n3)*exp(-I*phi3):

M322 := (`#mover(mi("n1"),mi("~"))`+n3)*exp(I*phi3):

M1 := (Matrix(2, 2, {(1, 1) = M111, (1, 2) = M112, (2, 1) = M121, (2, 2) = M122}))/(2*`#mover(mi("n2"),mi("~"))`):

M2 := (Matrix(2, 2, {(1, 1) = M211, (1, 2) = M212, (2, 1) = M221, (2, 2) = M222}))/(2*`#mover(mi("n3"),mi("~"))`):

M3 := (Matrix(2, 2, {(1, 1) = M311, (1, 2) = M312, (2, 1) = M321, (2, 2) = M322}))/(2*`#mover(mi("n1"),mi("~"))`):

Mcell := `.`(`.`(M3, M2), M1):

t := 1/Mcell[2, 2]:

r := t*Mcell[1, 2]:

h := abs(Re(Mcell[2, 2])):

R := abs(r)^2:

`ϕ` := arccos(Re(Mcell[2, 2])):

`ϕi` := arccosh(abs(Re(Mcell[2, 2]))):

`?N1` := sin(N*`ϕ`)/sin(`ϕ`):

`?N2` := sinh(N*`ϕi`)/sinh(`ϕi`):

RN1 := `?N1`^2*R/(1-R+`?N1`^2*R):

RN2 := `?N2`^2*R/(1-R+`?N2`^2*R):

f := piecewise(h <= 1, RN1, RN2):

> implicitplot3d(z-f, nu = 0 .. 6*`?b`, theta1 = 0 .. (1/2)*Pi, z = 0 .. 1,

tickmarks = [default, [seq((1/18)*i*Pi = (10*i)^`𝔬`, i = 1 .. 9)], default],

resolution = 10000,

axes = boxed,

style = patch,

orientation = [30,45],

lightmodel = light1,

labels = ["?", "?", "R"]);