Elgr2:=proc(m)    theta[0]:=m:
>     for h from 1 to 9 do
>     theta[h]:=theta[h-1]*(5-theta[h-1]):
>     od:
>     for i from -1 to -5 by -1 do
>     theta[i]:=(5-sqrt(25-4*theta[i+1]))/2:
>     od:
>     a[2]:=1: b[2]:=1-(theta[4])*(1-product( (1+4/(2-theta[n+4])), n=0..5))/4: c[2]:=b[2]:
>     for j from 1 to 0 by -1 do
>      a[j]:=((4-theta[2*j+1])*(a[j+1]+b[j+1])+2*c[j+1])/((2-theta[2*j+1])*(5-theta[2*j+1])):
> l[j]:=((4-theta[2*j+1])*(c[j+1]+b[j+1])+2*a[j+1])/((2-theta[2*j+1])*(5-theta[2*j+1])):
> b[j]:=((4-theta[2*j])*(a[j]+b[j+1])+2*l[j])/((2-theta[2*j])*(5-theta[2*j])):
> c[j]:=((4-theta[2*j])*(a[j]+l[j])+2*b[j+1])/((2-theta[2*j])*(5-theta[2*j])):
>     unassign('a[j+1]'): unassign('l[j]'): unassign('b[j+1]'): unassign('c[j+1]'):
>     od:
>     X:=array(1..6): Y:=array(2..6):
>     X[1]:=array(1..1): X[1][1]:=array(1..3): X[1][1][1]:=array(1..3,[0.,0.,1.]): X[1][1][2]:=array(1..3,[1.,0.,b[0]/a[0]]*1.): X[1][1][3]:=array(1..3,[1/2.,sqrt(3)/2.,c[0]/a[0]*1.]):
>     for k from 2 to 6 do
>     X[k]:=array(1..(3^(k-1))): Y[k]:=array(1..(3^(k-2))):
>     for l from 1 to (3^(k-2)) do
>     X[k][3*l-2]:=array(1..3): X[k][3*l-1]:=array(1..3): X[k][3*l]:=array(1..3):        Y[k][l]:=array(1..6):
>     Y[k][l][1]:=array(1..3,[X[k-1][l][1][1],X[k-1][l][1][2],X[k-1][l][1][3]]): Y[k][l][3]:=array(1..3,[X[k-1][l][2][1],X[k-1][l][2][2],X[k-1][l][2][3]]): Y[k][l][5]:=array(1..3,[X[k-1][l][3][1],X[k-1][l][3][2],X[k-1][l][3][3]]): Y[k][l][2]:=array(1..3,[(X[k-1][l][1][1]+X[k-1][l][2][1])/2,X[k-1][l][1][2],((4.-theta[1-k])*(X[k-1][l][1][3]+X[k-1][l][2][3])+2*X[k-1][l][3][3])/((2.-theta[1-k])*(5.-theta[1-k]))]): Y[k][l][4]:=array(1..3,[(X[k-1][l][2][1]+X[k-1][l][3][1])/2,(X[k-1][l][2][2]+X[k-1][l][3][2])/2,((4.-theta[1-k])*(X[k-1][l][2][3]+X[k-1][l][3][3])+2*X[k-1][l][1][3])/((2.-theta[1-k])*(5.-theta[1-k]))]): Y[k][l][6]:=array(1..3,[(X[k-1][l][3][1]+X[k-1][l][1][1])/2,(X[k-1][l][3][2]+X[k-1][l][1][2])/2,((4-theta[1-k])*(X[k-1][l][3][3]+X[k-1][l][1][3])+2*X[k-1][l][2][3])/((2.-theta[1-k])*(5.-theta[1-k]))]):
>     X[k][3*l-2]:=array(1..3): X[k][3*l-2][1]:=array(1..3,[Y[k][l][1][1],Y[k][l][1][2],Y[k][l][1][3]]): X[k][3*l-2][2]:=array(1..3,[Y[k][l][2][1],Y[k][l][2][2],Y[k][l][2][3]]): X[k][3*l-2][3]:=array(1..3,[Y[k][l][6][1],Y[k][l][6][2],Y[k][l][6][3]]):     X[k][3*l-1]:=array(1..3): X[k][3*l-1][1]:=array(1..3,[Y[k][l][2][1],Y[k][l][2][2],Y[k][l][2][3]]): X[k][3*l-1][2]:=array(1..3,[Y[k][l][3][1],Y[k][l][3][2],Y[k][l][3][3]]): X[k][3*l-1][3]:=array(1..3,[Y[k][l][4][1],Y[k][l][4][2],Y[k][l][4][3]]): X[k][3*l]:=array(1..3): X[k][3*l][1]:=array(1..3,[Y[k][l][6][1],Y[k][l][6][2],Y[k][l][6][3]]): X[k][3*l][2]:=array(1..3,[Y[k][l][4][1],Y[k][l][4][2],Y[k][l][4][3]]):
> X[k][3*l][3]:=array(1..3,[Y[k][l][5][1],Y[k][l][5][2],Y[k][l][5][3]]):
>     unassign('Y[k]'): unassign('X[k-2]'):
>     od:
>     od:
>     vals:={}:
>     for p from 1 to 243 do
>     vals:=vals union {array(1..3,[X[6][p][1][1],X[6][p][1][2],log(X[6][p][1][3])])}:
>     vals:=vals union {array(1..3,[X[6][p][2][1],X[6][p][2][2],log(X[6][p][2][3])])}:
>     vals:=vals union {array(1..3,[X[6][p][3][1],X[6][p][3][2],log(X[6][p][3][3])])}:
>     od:
>     with(plots):
>     pointplot3d(vals,axes=BOXED,title="Lambda0="||m);
> end: