%%% GENERATES UNITARY VECTORS - BT 050314

function  crds = PrtclGenMulti(ind, crds, Rsurf, sp, vp)

num = length(ind);


%% restart clock
crds(ind,9) = 0;


%% redefine latitude as angle from north pole
lat_polar = pi/2 - sp.lat(ind);


%% generate pseudorandom position and velocity
% weight theta distribution to compensate for polar coordinate distortion
if sp.SrfSpd <= pi  % all in one hemisphere
    theta = acos(1-(1-cos(sp.SrfSpd/2))*rand(num,1)); % here the sine dependence is integrated, normalized, and inversed
elseif sp.SrfSpd > pi  % split the integration since inverse cosine cannot yield half spreads greater than pi/2
    angle_section = (sp.SrfSpd-pi)/2;
    integral_section = 1 - cos(angle_section);
    rand_num = rand(num,1);
    if rand_num <= 1/(1+integral_section)
        theta = acos((1+integral_section)*rand_num);
    else
        theta = pi - acos(2-(1+integral_section)*rand_num);
    end
end
% weight theta2 distribution to compensate for polar coordinate distortion
if vp.AngSpd <= pi  % all in one hemisphere
    theta2 = acos(1-(1-cos(vp.AngSpd/2))*rand(num,1)); % here the sine dependence is integrated, normalized, and inversed
elseif vp.AngSpd > pi  % split the integration since inverse cosine cannot yield half spreads greater than pi/2
    angle_section = (sp.SrfSpd-pi)/2;
    integral_section = 1 - cos(angle_section);
    rand_num = rand(num,1);
    if rand_num <= 1/(1+integral_section)
        theta2 = acos((1+integral_section)*rand_num);
    else
        theta2 = pi - acos(2-(1+integral_section)*rand_num);
    end
end
phi = 2*pi*rand(num,1);
phi2 = 2*pi*rand(num,1);
%v_mag = vp.v0(ind) + vp.v_sd(ind).*sqrt(2)*erfinv(2*rand(num,1)-1);  %creates gaussian distribution using inverse error function


%% initialize an upward pointing jet at top of sphere
crds(ind,1) = Rsurf*sin(theta).*cos(phi); %x
crds(ind,2) = Rsurf*sin(theta).*sin(phi); %y
crds(ind,3) = Rsurf*cos(theta); %z
crds(ind,4) = sin(theta2).*cos(phi2);
crds(ind,5) = sin(theta2).*sin(phi2);
crds(ind,6) = cos(theta2);
% rot about Y
vx = crds(ind,4).*cos(theta) + crds(ind,6).*sin(theta);
vy = crds(ind,5);
vz = -crds(ind,4).*sin(theta) + crds(ind,6).*cos(theta);
% rot about Z
vx2 = vx.*cos(phi) - vy.*sin(phi);
vy2 =  vx.*sin(phi) + vy.*cos(phi);
vz2 = vz;
crds(ind,4:6) = [vx2 vy2 vz2];


%% tilt jet (about y-axis)
vx = crds(ind,4).*cos(vp.tlt(ind)) - crds(ind,6).*sin(vp.tlt(ind));
vy = crds(ind,5);
vz = crds(ind,4).*sin(vp.tlt(ind)) + crds(ind,6).*cos(vp.tlt(ind));
crds(ind,4:6) = [vx vy vz];


%% rotate jet to specified azimuth (i.e. about z-axis)
vx = crds(ind,4).*cos(vp.az(ind)) + crds(ind,5).*sin(vp.az(ind));
vy = -crds(ind,4).*sin(vp.az(ind)) + crds(ind,5).*cos(vp.az(ind));
vz = crds(ind,6);
crds(ind,4:6) = [vx vy vz];


%% rotate source to specified latitude (about y-axis)
x = crds(ind,1).*cos(lat_polar) - crds(ind,3).*sin(lat_polar);
y = crds(ind,2);
z = crds(ind,1).*sin(lat_polar) + crds(ind,3).*cos(lat_polar);
vx = crds(ind,4).*cos(lat_polar) - crds(ind,6).*sin(lat_polar);
vy = crds(ind,5);
vz = crds(ind,4).*sin(lat_polar) + crds(ind,6).*cos(lat_polar);
crds(ind,1:6) = [x y z vx vy vz];


%% rotate source to specified longitude (i.e. about z-axis)
x = crds(ind,1).*cos(-sp.lon(ind)) + crds(ind,2).*sin(-sp.lon(ind));
y = -crds(ind,1).*sin(-sp.lon(ind)) + crds(ind,2).*cos(-sp.lon(ind));
z = crds(ind,3);
vx = crds(ind,4).*cos(-sp.lon(ind)) + crds(ind,5).*sin(-sp.lon(ind));
vy = -crds(ind,4).*sin(-sp.lon(ind)) + crds(ind,5).*cos(-sp.lon(ind));
vz = crds(ind,6);
crds(ind,1:6) = [x y z vx vy vz];



end