Ho fatto un main di prova:
z=0.95;
n=100;
v_ion_modulo=2.867e4;
v_ion=repmat(v_ion_modulo,n,3);
v0=feval(@God_Speed,n,z,v_ion);
che chiama la funzione che ho sistemato così:
function [v0] = God_Speed(n,z, v_ion)
A=-1; B=1;
v=A+(B-A)*rand(200000,3); %cubo di lato 2 di numeri random uniformemente distribuiti di lato 2
% figure()
% plot3(v(:,1),v(:,2),v(:,3),'.b','MarkerSize',1)
d = v(:,1).^2 + v(:,2).^2 + v(:,3).^2; %equzione della sfera
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
% figure()
% plot3(v(:,1),v(:,2),v(:,3),'.y','MarkerSize',1)
% hold on
% plot3(condizione(:,1),condizione(:,2),condizione(:,3),'.k','MarkerSize', 1)
modulo_quadro=condizione(:,1).^2+condizione(:,2).^2+condizione(:,3).^2; %normalizzo i punti della semisfera in modo da ottenere una superfice
modulo=sqrt(modulo_quadro);
versore=condizione./repmat(modulo,1,3) ; % Calcolo i versori
% repmat: Replicate and tile array. B = repmat(A,m,n) creates a large matrix B
% consisting of an m-by-n tiling of copies of A. The size of B is [size(A,1)*m, (size(A,2)*n].
% Devo dividere le tre colonne di condizione per il vettore modulo,
% quindi faccio una matrice che contenga tre volte il vettore modulo per poter fare la divisione
% elemento per elemento. repmat(modulo,1,3) è uguale a scrivere [modulo modulo modulo].
figure()
plot3(versore(:,3),versore(:,2),versore(:,1), '.', 'MarkerSize',1)
axis equal
grid on
xlabel('beam direction')
ylabel('y axis')
zlabel('x axis')
% z=0.95; % angular spread %%%INPUT
% n=100; % 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à
% figure()
%
% plot3(vers_select(:,3),vers_select(:,2),vers_select(:,1), '.', 'MarkerSize',1)
% axis equal
% grid on
% xlabel('beam direction')
% ylabel('y axis')
% zlabel('x axis')
r = randperm(size(vers_select,1)); %creates a vector of random number between 1 and size(vers_select)
vers=vers_select(r(1:n),:); % versori per il calcolo delle componenti della velocità
figure()
plot3(vers(:,3),vers(:,2),vers(:,1), '.', 'MarkerSize',1)
axis equal
% v_ion_modulo=2.8e4; %Modulo delle velocità %%%INPUT
% v_ion=repmat(v_ion_modulo,n,3);
v0=vers.*v_ion; %Matrice delle velocità iniziali %%%OUTPUT
figure()
plot3(v0(:,3),v0(:,2),v0(:,1), '.', 'MarkerSize',1)
axis equal
e funziona, a meno che tu non vedi errori...
Mi rimane da capire come sistemare la domanda di prima.