% M. Sukop and E. Perfect, University of Kentucky, 2000
% Program computes, displays, and stores randomized Menger Sponges. 
% The homogeneous or heterogeneous (binomial/truncated binomial) algorithms 
% can be specified. 
% b is fractal scaling parameter (3 for standard Menger Sponge).
% p is the probability of a pore at each iteration level (7/27 for standard sponge).
% Output to Junk.dat is ASCII file of x,y, and z coordinates and 1's (solids) or 0's (pores)
% in the following order:
% row 1, column 1, layer 1; 
% row 1, column 1, layer 2; ...

clear
rand('state',sum(100*clock))

b=input('b value ? ')
p=input('p value ? ')
n=round(b^3-p*b^3)

D=log(n)/log(b)

maxit=input('Maximum number of iterations ? ')
typemode=input('Select fractal type: 0 = homogeneous; 1 = heterogeneous; 2 = non-random ')

if typemode == 1

	hettypemode=input('Select heterogeneous mode: 0 = regular binomial; 1 = truncated binomial ')

elseif typemode == 2
      
   a=input('Specify generator;    e.g., [1 1 1, 1 0 1, 1 1 1, 1 0 1, 0 0 0, 1 0 1, 1 1 1, 1 0 1, 1 1 1] for standard Menger sponge')
      
end


it=1
if typemode==0 %Homogeneous Algorithm

a=randperm(b^3);
for i = 1:b^3;
	if a(i)<=n 
	 	c(i)=1;%solid
	else 
		c(i)=0;%pore
	end
end
elseif typemode==1 %Heterogeneous Algorithm
   
   if hettypemode==0 %Regular binomial
      
   	a=rand(b^3,1);
		for i = 1:b^3;
			if a(i)<=1-p 
	 			c(i)=1;
			else 
      		c(i)=0;
      	end
      end
      
   else %Truncated binomial
      a=2;
      while min(a)>1-p
      	a=rand(b^3,1);
      end
		for i = 1:b^3;
			if a(i)<=1-p 
	 			c(i)=1;
			else 
      		c(i)=0;
      	end
   	end
   end
   
else %non-random
   c=a;
end %end algorithm selection


for row=1:b;
	for col=1:b;
      		for lay=1:b;
         			d(row,col,lay)=c(row*b^2-(b^2-b*col)-(b-lay));
      		end;
	end;
end;
oldcol=1

while it<maxit
	for oldrow=1:b^it;
		disp(100*oldrow*oldcol/b^(2*it))
		disp('% done current iteration')
		for oldcol=1:b^it;
      	for oldlay=1:b^it;
				if d(oldrow,oldcol,oldlay)==0 
					for newrow=oldrow*b-(b-1):oldrow*b;
            		for newcol=oldcol*b-(b-1):oldcol*b;
               		for newlay=oldlay*b-(b-1):oldlay*b;
              				e(newrow,newcol,newlay)=d(oldrow,oldcol,oldlay);
               		end   
						end
       			end  
            else	
               if typemode==0 %Homogeneous Algorithm
						a=randperm(b^3);
						for i = 1:b^3;
							if a(i)<=n 
		 						c(i)=1;
							else 
								c(i)=0;
							end
                  end
     			elseif typemode==1 %Heterogeneous Algorithm
   
				   if hettypemode==0 %Regular binomial
      
   					a=rand(b^3,1);
						for i = 1:b^3;
							if a(i)<=1-p 
					 			c(i)=1;
							else 
		      				c(i)=0;
	   	   			end
	  		    		end
      
   				else %Truncated binomial
   			   	a=2;
		      		while min(a)>1-p
			      		a=rand(b^3,1);
      				end
						for i = 1:b^3;
							if a(i)<=1-p 
	 							c(i)=1;
							else 
      						c(i)=0;
      					end
   					end
            	end
         	else %non-random
					c=a;
				end %end algorithm selection

						newrow=oldrow*b-(b-1);
         			newcol=oldcol*b-(b-1);
         			newlay=oldlay*b-(b-1);
						for row=1:b;
            			for col=1:b;
            				for lay=1:b;
                 				e(newrow+row-1,newcol+col-1,newlay+lay-1)=c(row*b^2-(b^2-b*col)-(b-lay));
  				     			end   
							end
						end
         		end
      		end   
			end
		end
		it=it+1
		d=e;
	end
	close
	figure(1)
hold on
fid = fopen('junk.dat','w');
col=0;
for row=1:b^maxit;
  	disp(100*row*col/b^(2*it))
	disp('% done generating plot')
 	for col=1:b^maxit;
   	for lay=1:b^maxit;
   		if d(row,col,lay)==1;
   			xmin=(row-1)/b^maxit;
   			xmax=row/b^maxit;
   			ymin=(col-1)/b^maxit;
   			ymax=col/b^maxit;
   			zmin=lay/b^maxit;
   			zmax=(lay-1)/b^maxit; 
fprintf(fid,'%6.5f %6.5f %6.5f \n',xmin+(xmax-xmin)/2,ymin+(ymax-ymin)/2,zmin+(zmax-zmin)/2); vert=[xmin,ymin,zmin;xmax,ymin,zmin; xmax,ymax,zmin; xmin,ymax,zmin; xmin,ymin,zmax; xmax,ymin,zmax; xmax,ymax,zmax; xmin,ymax,zmax];
fac=[1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8; 1 2 3 4; 5 6 7 8];
rgbcolors=[.5 .5 .5 ; 0 0 0; .5 .5 .5 ; 0 0 0; 1 1 1; 1 1 1];
%rgbcolors=[1 0 0 ; 1 0 0 ; 1 0 0  ; 1 0 0; 1 0 0 ; 1 0 0];

patch('faces',fac,'vertices',vert,'FaceVertexCData',rgbcolors,'FaceColor','flat','EdgeColor',[0 0 0],'BackFaceLighting', 'reverselit');
view(3)
			end;   
   	end;
	end;
end;
fclose(fid)
light('Position',[-1,0,0]);
light('Position',[0,-1,0]);
light('Position',[0,0,1]);
 
hold off

