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
  1 件のコメント
Varun
Varun 2023 年 2 月 23 日
Hello! Can you please clarify what your doubt is?

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeLanguage Fundamentals についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by