Salve ho scritto un codice che mi genera una figura geometrica di base quadrata, dandomi le coordinate dei punti (x,y,z) , i bond tra i vari punti ( come se fossero dei legami ) , sto avendo difficolta a generare un ciclo for che mi permetta di trovare le coordinate degli angoli generati tra questi bonds.
Con un esempio forse chiarisco l'idea.
Immaginiamo un cubo con 8 vertici, ognuno di essi crea un "legame" con il vicino e secondo vicino, creando il "lato" del cubo e la "diagonale" della faccia. adesso dovrei creare un codice che mi permetta di trovare angoli di 90, 45, 60 gradi tra i vari possibili con la condizione che questi angoli vengano formati da "bond con vertice comune" per esempio l'angolo generato tra un lato e diagonale o tra due diagonali. Il vertice comune non e' banale, basti pensare a figure geometriche complesse fatte da centinaia di punti. il cubo e' una semplificazione a titolo di esempio
Potreste aiutarmi a implementare il codice per favore?
clc;
clear;
diary('LAMMPS_bonds.txt')
diary on
dL=1; % distanza tra punti
Nx=3; % numero punti per asse
Ny=3; % numero punti per asse
Nz=3; % numero punti per asse
n = 1;
dr_Bond(1) = dL;
dr_Bond(2) = sqrt(2)*dL;
dummy = size(dr_Bond);
n_bond_types = dummy(2);
disp('Atoms Positions');
for i=1:Nx % ciclo for che genera coordinate per punti della mia figura
for j=1:Ny
for k=1:Nz
x(n) = i*dL;
y(n) = j*dL;
z(n) = k*dL;
if (i==1||i==Nx) % questa parte serve per differenziare i tipi di punti se stanno o no sulla superficie, utile nel secondo ciclo che genera diversi tipi di bond ( lati & diagonali)
sx(n) = 1; %atom on the surface x
else
sx(n) = 0; %atom not on the surface x
end
if (j==1||j==Ny)
sy(n) = 1; %atom on the surface y
else
sy(n) = 0;
end
if (k==1||k==Nz)
sz(n) = 1; %atom on the surface z
else
sz(n) = 0;
end
disp([n, x(n), y(n), z(n)]);
n = n+1;
end
end
end
Natoms = n - 1;
disp(['Number of Atoms= ', num2str(Natoms)]);
disp(' ');
nb = 1;
disp('Bonds');
for i=1:Natoms % combina i punti alla ricerca di lati e diagonali
for j=i+1:Natoms
dx = x(j)- x(i);
dy = y(j)- y(i);
dz = z(j)- z(i);
dr = sqrt(dx^2+dy^2+dz^2);
for k = 1:n_bond_types
if (dr> 0.9999*dr_Bond(k) && dr <1.0000001*dr_Bond(k))
b1(nb) = i;
b2(nb) = j;
btype(nb) = k;
bs(nb) = sx(i)*sx(j)+sy(i)*sy(j)+sz(i)*sz(j);%how many surface the bond is on
nb = nb+1;
end
end
end
end
Nbonds = nb - 1;
for n=1:Nbonds % differenzia lati e diagonali trovando quali esterni e interni alla figura
if (btype(n)==1)
if (bs(n)==0)
btype(n)= 1;
elseif (bs(n)==1)
btype(n)= 2;
elseif (bs(n)==2)
btype(n)= 3;
else
disp('something wrong with bs');
end
else
if (bs(n)==0)
btype(n)= 4;
elseif (bs(n)==1)
btype(n)= 5;
else
disp('something wrong with bs');
end
end
disp([n, btype(n), b1(n), b1(n)]);
end
disp(['Number of bonds= ', num2str(Nbonds)]);
disp('');
disp('Note: ');
disp('Type 1: side spring in the bulk (k=k1x4)');
disp('Type 2: side spring on the surface (k=k1x2)');
disp('Type 3: side spring on the edge (k=k1)');
disp('Type 4: diagonal spring in the bulk (k=k2x2)');
disp('Type 5: diagonal spring on the surface (k=k2)');
diary off