# Galaxy - random sampling issue

1. Aug 2, 2013

### fab13

I need to generate initial conditions for modeling galactic spiral arms.

rho = a. / (log (b * tanh (theta / (2 * n))) with a, b ​​and n are parameters to choose from.

to give a thickness along the curve for the generated points, I did the following things:

- I generate points in a rectangle of width hx and height hy.
- I divided the polar curve in several segments delta(theta) then I calculate the value of the tangent vector to the polar curve at each point of the segments
- I calculate the angle between the tangent vector and the vector ex (canonical basis vector ex)
- I apply the rotation matrix with this angle to the set of points generated.
- I applies a translation of the new points with the angle theta of the considered point (the point on the curve corresponding to the segment in question)

Note that I managed to reduce the hy thickness gradually (as theta increases).

but now, I can't produce a thickness for small theta, i.e reproduce the bar in the middle that joins the two spiral arms. You can see the result on the following figure :

You can see that there are artifacts on the figure for small theta.

Here's the original matlab source for this figure :

Code (Text):

function test_arms

% Total number of particles
num_particles=25600;
% Number of segments
num_seg=800;
% Number of particules per segment
num_part_seg=num_particles/(2*num_seg);

% Parameters for curbe
a1=10;
b1=0.5;
n1=4;

% Theta max
theta_max=2.25*pi;
% Width along curve
width=2;
% Hx width
hx=theta_max/num_seg;
% Hy width : decreasing with theta
hy=width-width.*(1:num_seg)/num_seg;
% Theta array
t=0:hx:theta_max-hx;
% Function values
f1=a1./(log(b1*tanh(2*pi*t/(theta_max*(2*n1)))));
% Derivates
f11_prim=diff(f1)./diff(t);
f12_prim=diff(-f1)./diff(t);

f11_prim(num_seg)=0;
f12_prim(num_seg)=0;
% Tangent vector coordinates
x_dom=(f11_prim-f1.*sin(t));
y_dom=(f12_prim+f1.*cos(t));
% Angle value with ex
alpha11=-atan(y_dom./x_dom);
alpha12=atan(-y_dom./x_dom);

for i=1:num_seg

% Random box coordinates for f1
x11_first=hx*(rand(num_part_seg,1)-0.5);
y11_first=hy(i)*(rand(num_part_seg,1)-0.5);

% Rotation matrix + translation
x11((i-1)*num_part_seg+1:i*num_part_seg)=(x11_first*cos(alpha11(i))+y11_first*sin(alpha11(i)))+f1(i)*cos(t(i));
y11((i-1)*num_part_seg+1:i*num_part_seg)=(-x11_first*sin(alpha11(i))+y11_first*cos(alpha11(i)))+f1(i)*sin(t(i));

% Random box coordinates for -f1
x12_first=hx*(rand(num_part_seg,1));
y12_first=hy(i)*(rand(num_part_seg,1)-0.5);

% Rotation matrix + translation
x12((i-1)*num_part_seg+1:i*num_part_seg)=(x12_first*cos(alpha12(i))+y12_first*sin(alpha12(i)))-f1(i)*cos(t(i));
y12((i-1)*num_part_seg+1:i*num_part_seg)=(-x12_first*sin(alpha12(i))+y12_first*cos(alpha12(i)))-f1(i)*sin(t(i));

end

figure(1);
polar(t,f1);
hold on;
polar(t,-f1);
hold on;
scatter(x11,y11,1);
hold on;
scatter(x12,y12,1);

end

I thought this issue comes from the subsampling of my segments for little theta, so I tried to fix that with this modified source :

Code (Text):

function test_arms_2

% Total number of particles
num_particles=50000;
% Number of segments
num_seg=5000;
% Number of particules per segment
num_part_seg=num_particles/(2*num_seg);

% Parameters for curbe
a1=10;
b1=0.5;
n1=4;

%num_seg of bar
num_seg_bar=2000;
num_seg_spiral=num_seg-num_seg_bar;

% Theta min
theta_min=0.25*pi;
% Theta max
theta_max=2.25*pi;
% Width along curve
width=2;
% Hx bar width
hx_bar=theta_min/num_seg_bar;
% Hx spiral width
hx_spiral=(theta_max-theta_min)/num_seg_spiral;
% Hy width : decreasing with theta
hy=width-width.*(1:num_seg)/num_seg;
% Theta array
t=[linspace(0,theta_min,num_seg_bar) linspace(theta_min,theta_max,num_seg_spiral)];
hx(1:num_seg_bar)=hx_bar;
hx(num_seg_bar+1:num_seg)=hx_spiral;
% Function values
f1=a1./(log(b1*tanh(2*pi*t/(theta_max*(2*n1)))));
% Derivates
f11_prim=diff(f1)./diff(t);
f12_prim=diff(-f1)./diff(t);

f11_prim(num_seg)=0;
f12_prim(num_seg)=0;
% Tangent vector coordinates
x_dom=(f11_prim-f1.*sin(t));
y_dom=(f12_prim+f1.*cos(t));
% Angle value with ex
alpha11=-atan(y_dom./x_dom);
alpha12=atan(-y_dom./x_dom);

for i=1:num_seg

% Random box coordinates for f1
x11_first=hx(i).*(rand(num_part_seg,1)-0.5);
y11_first=hy(i)*(rand(num_part_seg,1)-0.5);

% Rotation matrix + translation
x11((i-1)*num_part_seg+1:i*num_part_seg)=(x11_first*cos(alpha11(i))+y11_first*sin(alpha11(i)))+f1(i)*cos(t(i));
y11((i-1)*num_part_seg+1:i*num_part_seg)=(-x11_first*sin(alpha11(i))+y11_first*cos(alpha11(i)))+f1(i)*sin(t(i));

% Random box coordinates for -f1
x12_first=hx(i).*(rand(num_part_seg,1));
y12_first=hy(i)*(rand(num_part_seg,1)-0.5);

% Rotation matrix + translation
x12((i-1)*num_part_seg+1:i*num_part_seg)=(x12_first*cos(alpha12(i))+y12_first*sin(alpha12(i)))-f1(i)*cos(t(i));
y12((i-1)*num_part_seg+1:i*num_part_seg)=(-x12_first*sin(alpha12(i))+y12_first*cos(alpha12(i)))-f1(i)*sin(t(i));

end

figure(1);
polar(t,f1);
hold on;
polar(t,-f1);
hold on;
scatter(x11,y11,1);
hold on;
scatter(x12,y12,1);

end

Anyone could see how to reproduce this central bar... ?

Thanks

2. Aug 2, 2013

### Staff: Mentor

Where is the point in those rotated squares (I think their overlap generates the dense regions, and probably some other issues)? You can use coordinates in a circle, then you don't need those rotations. Just make sure that your sampling of the curve is much denser than the width of the curve there. To improve this: the outer regions can use circles with only a few points inside.
That simplifies the code significantly, it gets rid some strange features and it is easier to change things.