- #1
fab13
- 312
- 6
I need to generate initial conditions for modeling galactic spiral arms.
I start with the following polar equation:
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 :
I thought this issue comes from the subsampling of my segments for little theta, so I tried to fix that with this modified source :
Anyone could see how to reproduce this central bar... ?
Thanks
I start with the following polar equation:
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:
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:
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