conversion calculation and optimal calculation
11 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
tic
p=0.0; %percentage of dumbbells
ps=1-p; %percentage of spheres
intruder_diam=5;
r=(0.5/(2^(0.5)));%dumbell radius
Area_dumb=2*pi*r*r;
r_s=(Area_dumb/pi)^0.5;%Multiply by 1.1 to set to maximum radius
NoOfParticles=180000;
Lx_t=[-150 150]; Ly_t=[-600 600];
Lx=[Lx_t(1)+2*r_s Lx_t(2)-2*r_s];
Ly=[Ly_t(1)+2*r_s Ly_t(2)-2*r_s];
%Predefined
rad_wall=0.5;
trash1=(Lx_t(1)+rad_wall:(2*rad_wall):Lx_t(2)-rad_wall)';
walls1=[trash1 rad_wall+Ly_t(1)+0*trash1 rad_wall+0*trash1];
walls2=[trash1 -rad_wall+Ly_t(2)+0*trash1 rad_wall+0*trash1];
walls=[walls1; walls2];
NoOfDumbs=floor(p*NoOfParticles);
NoOfSpheres=NoOfParticles-NoOfDumbs;
template=[-r 0; r 0];
template2=[0 0];
NoOfTrials=5000000;
%config=[template [r; r]];
%configuration
NoOfWallParts=length(walls(:,1));
config=zeros(NoOfDumbs*2+NoOfSpheres+NoOfWallParts+1,3);
config(2:NoOfWallParts+1,:)=walls;
config(1,:)=[0 0 intruder_diam/2];
pos=randi([0,1000000000000],NoOfTrials,2)/10^12;
pos(:,1)=pos(:,1)*diff(Lx)+Lx(1);
pos(:,2)=pos(:,2)*diff(Ly)+Ly(1);
angle=rand(NoOfTrials,1)*pi;
radsdumb=(0.2*(rand(NoOfTrials,1)-0.5)+1)*r;
rads=(0.2*(rand(NoOfTrials,1)-0.5)+1)*r_s;
counter=0;
label=NoOfWallParts+1;
for i=1:NoOfTrials
if(counter<NoOfDumbs)
coords=[pos(i,1)+radsdumb(i)*cos(angle(i)) pos(i,2)+radsdumb(i)*sin(angle(i)) radsdumb(i); pos(i,1)-radsdumb(i)*cos(angle(i)) pos(i,2)-radsdumb(i)*sin(angle(i)) radsdumb(i)];
arr=[(config(1:label,1)-coords(1,1)).^2+(config(1:label,2)-coords(1,2)).^2-(config(1:label,3)+coords(1,3)).^2; (config(1:label,1)-coords(2,1)).^2+(config(1:label,2)-coords(2,2)).^2-(config(1:label,3)+coords(2,3)).^2];
minim=min(arr);
if(minim>=0)
init=label+1;
label=label+2;
counter=counter+1;
config(init:label,:)=coords;
end
elseif (counter>=NoOfDumbs && counter<NoOfParticles)
coords=[pos(i,1) pos(i,2) rads(i)];
arr=[(config(1:label,1)-coords(1,1)).^2+(config(1:label,2)-coords(1,2)).^2-(config(1:label,3)+coords(1,3)).^2];
minim=min(arr);
if(minim>=0)
label=label+1;
counter=counter+1;
config(label,:)=coords;
end
end
end
checker=[1;2*ones(NoOfWallParts,1);3*ones(2*NoOfDumbs,1);4*ones(NoOfSpheres,1)];
Final=[checker config];
%Final=[ones(NoOfDumbs*2,1) config(NoOfWallParts+1:NoOfWallParts+2*NoOfDumbs,:); 2*ones(NoOfSpheres,1) config(NoOfWallParts+1+2*NoOfDumbs:NoOfWallParts+2*NoOfDumbs+NoOfSpheres,:); [3*ones(NoOfWallParts/2,1); 4*ones(NoOfWallParts/2,1)] walls];
lammps=[(1:length(Final(:,1)))' Final(:,1) Final(:,4)*2 1+Final(:,1)*0 Final(:,2:3) Final(:,1)*0 Final(:,1)*0 Final(:,1)*0 Final(:,1)*0];
molecules0=[(1:(NoOfWallParts+1))' (1:(NoOfWallParts+1))'];
molecules1=[NoOfWallParts+1+(1:2*NoOfDumbs)' NoOfWallParts+1+floor((2:1+2*NoOfDumbs)'/2)];
%molecules2=[(2*NoOfDumbs+1:2*NoOfDumbs+NoOfSpheres+NoOfWallParts)' floor((1+2*NoOfDumbs)'/2)+(1:NoOfSpheres+NoOfWallParts)'];
molecules2=[NoOfWallParts+1+2*NoOfDumbs+(1:NoOfSpheres)' NoOfWallParts+1+NoOfDumbs+(1:NoOfSpheres)'];
molecule=[molecules0;molecules1;molecules2];
toc
1 commentaire
Réponses (0)
Voir également
Catégories
En savoir plus sur Language Fundamentals dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!