User:Zeracles/Background/Code/Structures

See also this quick download

%	structures.m
%	by Anthony Smith (Zeracles); astroant@hotmail.com, asmith@physics.usyd.edu.au, http://www.physics.usyd.edu.au/~asmith/, http://wiki.uqm.stack.nl/User:Zeracles
%	see http://starcontrol.classicgaming.gamespy.com/forum/index.php?topic=1660.0

plotpoints=1;
savepoints=1;
volumetype=1; % 1 for a cube, 2 for a sphere, 3 for a cylinder

%	type
%	1	single point
%	2	gaussian
%	3	hard-edged sphere
%	4	sheet / disc
%	5	filament

%	format of input file by column is:	type, x, y, z, scale, density, outerradius (types 1,2,3)
%						type, x, y, z, (half-) thickness, density, radius, azimuth, elevation (type 4)
%							azimuth relative to positive x-axis, elevation relative to negative z-axis, in radians
%						type, x1, y1, z1, x2, y2, z2, number of points (type 5)

load('structureproperties')

sizestructureproperties=size(structureproperties);
numberofstructures=sizestructureproperties(1);

if volumetype==1
	xmin=0;
	xmax=1;
	ymin=0;
	ymax=1;
	zmin=0;
	zmax=1;
	
	xextent=xmax-xmin;
	yextent=ymax-ymin;
	zextent=zmax-zmin;
	
	volume=xextent*yextent*zextent;
end
if volumetype==2
	centrex=0.5;
	centrey=0.5;
	centrez=0.5;
	
	radius=0.5;

	volume=(4/3)*pi*radius^3;
end
if volumetype==3
	centrex=0.5;
	centrey=0.5;
	
	radius=0.5;
	
	zmin=0;
	zmax=1;
	
	zextent=zmax-zmin;
	
	volume=pi*radius^2*zextent;
end

mapx=[];
mapy=[];
mapz=[];

for i=linspace(1,numberofstructures,numberofstructures)
	structureinfo=structureproperties(i,:);
	
	structuretype=structureinfo(1);
	structurex=structureinfo(2);
	structurey=structureinfo(3);
	structurez=structureinfo(4);
	structurescale=structureinfo(5);
	structuredensity=structureinfo(6);
	structureouterradius=structureinfo(7);
	if structuretype==1
		structurexs=structurex;
		structureys=structurey;
		structurezs=structurez;
	end
	if structuretype==2
		initialrandomx=rand(structuredensity,1);
		initialrandomy=rand(structuredensity,1);
		initialrandomz=rand(structuredensity,1);
		
		clusterboxxmin=structurex-structureouterradius;
		clusterboxxmax=structurex+structureouterradius;
		clusterboxymin=structurey-structureouterradius;
		clusterboxymax=structurey+structureouterradius;
		clusterboxzmin=structurez-structureouterradius;
		clusterboxzmax=structurez+structureouterradius;
		
		clusterboxmask=(clusterboxxmin<=initialrandomx).*(initialrandomx<=clusterboxxmax).*(clusterboxymin<=initialrandomy).*(initialrandomy<=clusterboxymax).*(clusterboxzmin<=initialrandomz).*(initialrandomz<=clusterboxzmax);
		
		clusterboxrandomx=initialrandomx(find(clusterboxmask>0));
		clusterboxrandomy=initialrandomy(find(clusterboxmask>0));
		clusterboxrandomz=initialrandomz(find(clusterboxmask>0));
		
		distance=((clusterboxrandomx-structurex).^2+(clusterboxrandomy-structurey).^2+(clusterboxrandomz-structurez).^2).^(1/2);
		
		clusterspheremask=distance<=structureouterradius;
		
		clustersphererandomx=clusterboxrandomx(find(clusterspheremask>0));
		clustersphererandomy=clusterboxrandomy(find(clusterspheremask>0));
		clustersphererandomz=clusterboxrandomz(find(clusterspheremask>0));
		
		distance=nonzeros(clusterspheremask.*distance);
		
		admissionprobability=rand(sum(clusterspheremask),1);
		probabilitythreshold=1-exp(-(distance.^2)./(2*structurescale^2));
		admissionmask=admissionprobability>probabilitythreshold;
		
		structurexs=clustersphererandomx(find(admissionmask>0));
		structureys=clustersphererandomy(find(admissionmask>0));
		structurezs=clustersphererandomz(find(admissionmask>0));
	end
	if structuretype==3
		initialrandomx=rand(structuredensity,1);
		initialrandomy=rand(structuredensity,1);
		initialrandomz=rand(structuredensity,1);
		
		clusterboxxmin=clusterx-structureouterradius;
		clusterboxxmax=clusterx+structureouterradius;
		clusterboxymin=clustery-structureouterradius;
		clusterboxymax=clustery+structureouterradius;
		clusterboxzmin=clusterz-structureouterradius;
		clusterboxzmax=clusterz+structureouterradius;
		
		clusterboxmask=(clusterboxxmin<=initialrandomx).*(initialrandomx<=clusterboxxmax).*(clusterboxymin<=initialrandomy).*(initialrandomy<=clusterboxymax).*(clusterboxzmin<=initialrandomz).*(initialrandomz<=clusterboxzmax);
		
		clusterboxrandomx=initialrandomx(find(clusterboxmask>0));
		clusterboxrandomy=initialrandomy(find(clusterboxmask>0));
		clusterboxrandomz=initialrandomz(find(clusterboxmask>0));
		
		distance=((clusterboxrandomx-structurex).^2+(clusterboxrandomy-structurey).^2+(clusterboxrandomz-structurez).^2).^(1/2);
		
		clusterspheremask=distance<=structureouterradius;
		
		structurexs=clusterboxrandomx(find(clusterspheremask>0));
		structureys=clusterboxrandomy(find(clusterspheremask>0));
		structurezs=clusterboxrandomz(find(clusterspheremask>0));
	end
	if structuretype==4
		initialrandomx=rand(structuredensity,1);
		initialrandomy=rand(structuredensity,1);
		initialrandomz=rand(structuredensity,1);
		
		clusterboxxmin=structurex-structureouterradius;
		clusterboxxmax=structurex+structureouterradius;
		clusterboxymin=structurey-structureouterradius;
		clusterboxymax=structurey+structureouterradius;
		clusterboxzmin=structurez-structureouterradius;
		clusterboxzmax=structurez+structureouterradius;
		
		clusterboxmask=(clusterboxxmin<=initialrandomx).*(initialrandomx<=clusterboxxmax).*(clusterboxymin<=initialrandomy).*(initialrandomy<=clusterboxymax).*(clusterboxzmin<=initialrandomz).*(initialrandomz<=clusterboxzmax);
		
		clusterboxrandomx=initialrandomx(find(clusterboxmask>0));
		clusterboxrandomy=initialrandomy(find(clusterboxmask>0));
		clusterboxrandomz=initialrandomz(find(clusterboxmask>0));
		
		distance=((clusterboxrandomx-structurex).^2+(clusterboxrandomy-structurey).^2+(clusterboxrandomz-structurez).^2).^(1/2);
		
		clusterspheremask=distance<=structureouterradius;
		
		spherexs=clusterboxrandomx(find(clusterspheremask>0));
		sphereys=clusterboxrandomy(find(clusterspheremask>0));
		spherezs=clusterboxrandomz(find(clusterspheremask>0));
		
		azimuth=structureinfo(8)+pi;
		elevation=structureinfo(9);
		
		planevector=[cos(azimuth)*sin(elevation);sin(azimuth)*sin(elevation);cos(elevation)];
		
		a=planevector(1);
		b=planevector(2);
		c=planevector(3);
		
		d=-dot(planevector,[structurex;structurey;structurez]);
		
		perpendicularplanedistance=(abs(a*spherexs+b*sphereys+c*spherezs+d))./(a^2+b^2+c^2);
		planemask=perpendicularplanedistance<=structurescale;
		
		structurexs=spherexs(find(planemask>0.5));
		structureys=sphereys(find(planemask>0.5));
		structurezs=spherezs(find(planemask>0.5));
	end
	if structuretype==5
		numberofpoints=structureinfo(8);
		
		structurexs=transpose(linspace(structurex,structureinfo(5),numberofpoints));
		structureys=transpose(linspace(structurey,structureinfo(6),numberofpoints));
		structurezs=transpose(linspace(structurez,structureinfo(7),numberofpoints));
	end
	mapx=[mapx;structurexs];
	mapy=[mapy;structureys];
	mapz=[mapz;structurezs];
end
if volumetype==1
	contained=(xmin<=mapx).*(mapx<=xmax).*(ymin<=mapy).*(mapy<=ymax).*(zmin<=mapz).*(mapz<=zmax);
	
	mapx=mapx(find(contained>0));
	mapy=mapy(find(contained>0));
	mapz=mapz(find(contained>0));
end
if volumetype==2
	contained=(((mapx-centrex).^2+(mapy-centrey).^2+(mapz-centrez).^2).^(1/2))<=radius;
	
	mapx=mapx(find(contained>0));
	mapy=mapy(find(contained>0));
	mapz=mapz(find(contained>0));
end
if volumetype==3
	contained=((((mapx-centrex).^2+(mapy-centrey).^2).^(1/2))<=radius)+((zmin<=mapz).*(mapz<=zmax));
	
	mapx=mapx(find(contained>0));
	mapy=mapy(find(contained>0));
	mapz=mapz(find(contained>0));
end
if plotpoints==1
	figure
	hold on
	box on
	plot3(mapx,mapy,mapz,'k.')
end
if savepoints==1
	structuredpoints=[mapx,mapy,mapz];
	save('structuredpoints','structuredpoints','-ascii')
end