giug ha scritto:
in quale riga? (se mi metti tutte le informazioni insieme facciamo prima...)
si incazza subito prima del plot finale, perchè n resta più grande di quello che ha estratto due volte...
A=-1; B=1;
v=A+(B-A)*rand(2000000,3); %cubo uniforme
d = v(:,1).^2 + v(:,2).^2 + v(:,3).^2;
posizioni=find(d<=1 & v(:,3)>=0 & v(:,1)~=0 & v(:,2)~=0 ); %cerco quei punti che stanno nella semisfera di raggio 1 e con z>=0 e x,y ~=0
condizione=v(posizioni,:); %salvo i punti che stanno nella sfera
modulo_quadro=condizione(:,1).^2+condizione(:,2).^2+condizione(:,3).^2;
modulo=sqrt(modulo_quadro);
versore=condizione./repmat(modulo,1,3) ;
z=0.999; %angular spread %%%INPUT
n=10000; %numero diparticelle %%%INPUT
vers=find(versore(:,3)>=z); %Selezioni i punti che stanno in una calotta sferica
vers_select=versore(vers,:); %creo il vettore con i versori da usare per le proiezioni della velocità
r = randperm(size(vers_select,1)); %creates a vector of random number between 1 and size(vers_select)
if n>length(vers_select)
r1=randperm(n-size(vers_select,1));
r=[r,r1];
end
v=vers_select(r(1:n),:); % versori per il calcolo delle componenti della velocità
figure()
plot3(v(:,3),v(:,2),v(:,1), '.', 'MarkerSize',1)
axis equal