conversion calculation and optimal calculation

11 次查看（过去 30 天）
Racharla Vignesh
Racharla Vignesh 2022-11-25
评论： Varun 2023-2-23
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

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Language Fundamentals 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by