function data = elgr(M)
    data  = zeros(729*3, 3);
    theta = zeros(16,1);
    a = zeros(3,1);
    b = zeros(3,1);
    c = zeros(3,1);
    
    theta(7) = M;
    
    for h=1:9
        theta(7+h)=theta(h+6)*(5-theta(h+6));
    end
    for h=-1:-1:-6
        theta(7+h) = (5-sqrt(25-4*theta(h+8)))/2;
    end
    
    a(3) = 1;
    t = 1;
    for n=7:12
        t = t * (1+4/(2-theta(n+4)));
    end
    b(3) = 1-(theta(11))*(1-t)/4;
    c(3) = b(3);
    for j=2:-1:1
        a(j) =((4-theta(2*j+8))*(a(j+1)+b(j+1))+2*c(j+1))/((2-theta(2*j+8))*(5-theta(2*j+8)));
        l(j) =((4-theta(2*j+8))*(c(j+1)+b(j+1))+2*a(j+1))/((2-theta(2*j+8))*(5-theta(2*j+8)));
        b(j) =((4-theta(2*j+7))*(a(j)+b(j+1))+2*l(j))/((2-theta(2*j+7))*(5-theta(2*j+7)));
        c(j) =((4-theta(2*j+7))*(a(j)+l(j))+2*b(j+1))/((2-theta(2*j+7))*(5-theta(2*j+7)));
    end
    
    X = zeros(7, 3^6, 3, 3);
    Y = zeros(7, 3^5, 6, 3);
    X(1,1,1,:) = [0,0,1];
    X(1,1,2,:) = [1,0,b(1)/a(1)];
    X(1,1,3,:) = [.5,sqrt(3)/2,c(1)/a(1)];
    
    for i = 2:7
        for j=1:3^(i-2)            
            Y(i,j,1,:) = [X(i-1,j,1,1), X(i-1,j,1,2), X(i-1,j,1,3)];
            Y(i,j,3,:) = [X(i-1,j,2,1), X(i-1,j,2,2), X(i-1,j,2,3)];
            Y(i,j,5,:) = [X(i-1,j,3,1), X(i-1,j,3,2), X(i-1,j,3,3)];            
            Y(i,j,2,:) = [(X(i-1,j,1,1) + X(i-1,j,2,1))/2, X(i-1,j,1,2), ...
                          ((4-theta(8-i))*(X(i-1,j,1,3)+X(i-1,j,2,3))+2*X(i-1,j,3,3))/((2-theta(8-i))*(5-theta(8-i)))];
            Y(i,j,4,:) = [(X(i-1,j,2,1)+X(i-1,j,3,1))/2,(X(i-1,j,2,2)+X(i-1,j,3,2))/2,...
                          ((4-theta(8-i))*(X(i-1,j,2,3)+X(i-1,j,3,3))+2*X(i-1,j,1,3))/((2-theta(8-i))*(5-theta(8-i)))];
            Y(i,j,6,:) = [(X(i-1,j,3,1)+X(i-1,j,1,1))/2,(X(i-1,j,3,2)+X(i-1,j,1,2))/2,...
                          ((4-theta(8-i))*(X(i-1,j,3,3)+X(i-1,j,1,3))+2*X(i-1,j,2,3))/((2-theta(8-i))*(5-theta(8-i)))];
    
           X(i,3*j-2,1,:) = [Y(i,j,1,1),Y(i,j,1,2),Y(i,j,1,3)];
           X(i,3*j-2,2,:) = [Y(i,j,2,1),Y(i,j,2,2),Y(i,j,2,3)];
           X(i,3*j-2,3,:) = [Y(i,j,6,1),Y(i,j,6,2),Y(i,j,6,3)];
           
           X(i,3*j-1,1,:) = [Y(i,j,2,1),Y(i,j,2,2),Y(i,j,2,3)];
           X(i,3*j-1,2,:) = [Y(i,j,3,1),Y(i,j,3,2),Y(i,j,3,3)];
           X(i,3*j-1,3,:) = [Y(i,j,4,1),Y(i,j,4,2),Y(i,j,4,3)];
           
           X(i,3*j,1,:) = [Y(i,j,6,1),Y(i,j,6,2),Y(i,j,6,3)];
           X(i,3*j,2,:) = [Y(i,j,4,1),Y(i,j,4,2),Y(i,j,4,3)];
           X(i,3*j,3,:) = [Y(i,j,5,1),Y(i,j,5,2),Y(i,j,5,3)];                           
        end
    end
    
    for i = 0:728
        data(3*i + 1, :) = X(7, i+1, 1, :);
        data(3*i + 2, :) = X(7, i+1, 2, :);
        data(3*i + 3, :) = X(7, i+1, 3, :);
    end
    plot3(data(:,1),data(:,2),data(:,3), '.', 'MarkerSize', 1);
    
function out = ThetaIndex(i)  
    out = 7+i;