Zeracles: clarification

2009-09-17T06:46:58Z

clarification

← Older revision		Revision as of 06:46, 17 September 2009
Line 16:		Line 16:

	% format of input file by column is: type, x, y, z, scale, density, outerradius (types 1,2,3)		% format of input file by column is: type, x, y, z, scale, density, outerradius (types 1,2,3)
	% type, x, y, z, thickness, density, radius, azimuth, elevation (type 4)		% 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		% 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)		% type, x1, y1, z1, x2, y2, z2, number of points (type 5)

Zeracles: structure generation in matlab

2009-09-17T03:48:40Z

structure generation in matlab

New page

See also [http://www.physics.usyd.edu.au/~asmith/files/structures.m 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, 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

User:Zeracles/Background/Code/Structures - Revision history

Zeracles: clarification

Zeracles: structure generation in matlab