In [303]:
%pylab inline
In [304]:
from IPython.display import display
from IPython.display import Image
In [305]:
Image(url="http://upload.wikimedia.org/wikipedia/commons/thumb/0/0b/Gravitationell-lins-4.jpg/290px-Gravitationell-lins-4.jpg")
Out[305]:
In [306]:
Image(url='http://www.space2099theseries.com/blogs/attachment46d1326840084-chandra-observatory-finds-largest-galaxy-cluster-universe-clusters_xray.jpg')
Out[306]:
In [435]:
data= loadtxt('../data/MDR1_BDMW_snap_85_z_200.0_240.0.csv', comments='#', skiprows=37, delimiter=',')
In [442]:
axis_ratio_2_1 = data[:,21]
good_signal_id = where(axis_ratio_2_1[:]<0.7)
data = data[good_signal_id[0],:]
x_dir = data[:,23]
y_dir = data[:,24]
z_dir = data[:,25]
x_pos = data[:,4]
y_pos = data[:,5]
z_pos = data[:,6]
print shape(data)
In [437]:
def get_alignment_stats(x_pos, y_pos, z_pos, x_dir, y_dir, z_dir):
n_points = size(x_pos)
dist_pair_2D = empty((0))
dist_pair_3D = empty((0))
cos_theta_pair_all_2D = empty((0))
cos_theta_pair_all_3D = empty((0))
for i in range(n_points-1):
my_pos_x = x_pos[i]
my_pos_y = y_pos[i]
my_pos_z = y_pos[i]
delta_x = 500.0 - my_pos_x
delta_y = 500.0 - my_pos_y
delta_z = 500.0 - my_pos_z
my_pos_x = my_pos_x + delta_x
my_pos_y = my_pos_y + delta_y
my_pos_z = my_pos_z + delta_z
tmp_pos_x = (x_pos[i+1:] + delta_x)%1000.0
tmp_pos_y = (y_pos[i+1:] + delta_y)%1000.0
tmp_pos_z = (z_pos[i+1:] + delta_z)%1000.0
tmp_dir_x = x_dir[i+1:]
tmp_dir_y = y_dir[i+1:]
tmp_dir_z = z_dir[i+1:]
my_dir_x = x_dir[i]
my_dir_y = y_dir[i]
my_dir_z = z_dir[i]
norm_all = sqrt(tmp_dir_x**2 + tmp_dir_y**2)
my_norm = sqrt(my_dir_x**2 + my_dir_y**2)
cos_theta_pair_2D = (tmp_dir_x*my_dir_x + tmp_dir_y*my_dir_y)/(norm_all * my_norm)
cos_theta_pair_3D = (tmp_dir_x*my_dir_x + tmp_dir_y*my_dir_y + tmp_dir_z*my_dir_z)
radius_2D = sqrt((my_pos_x - tmp_pos_x)**2 + (my_pos_y-tmp_pos_y)**2)
radius_3D = sqrt((my_pos_x - tmp_pos_x)**2 + (my_pos_y-tmp_pos_y)**2 + (my_pos_z-tmp_pos_z)**2)
dist_pair_2D = append(dist_pair_2D, radius_2D)
cos_theta_pair_all_2D=append(cos_theta_pair_all_2D, cos_theta_pair_2D**2)
dist_pair_3D = append(dist_pair_3D, radius_3D)
cos_theta_pair_all_3D=append(cos_theta_pair_all_3D, cos_theta_pair_3D**2)
return [dist_pair_2D, cos_theta_pair_all_2D, dist_pair_3D, cos_theta_pair_all_3D]
In [438]:
def get_pointing_stats(x_pos, y_pos, z_pos, x_dir, y_dir, z_dir):
n_points = size(x_pos)
dist_pair_2D = empty((0))
dist_pair_3D = empty((0))
cos_theta_pair_all_2D = empty((0))
cos_theta_pair_all_3D = empty((0))
for i in range(n_points-1):
my_pos_x = x_pos[i]
my_pos_y = y_pos[i]
my_pos_z = y_pos[i]
delta_x = 500.0 - my_pos_x
delta_y = 500.0 - my_pos_y
delta_z = 500.0 - my_pos_z
my_pos_x = my_pos_x + delta_x
my_pos_y = my_pos_y + delta_y
my_pos_z = my_pos_z + delta_z
tmp_pos_x = (x_pos[i+1:] + delta_x)%1000.0
tmp_pos_y = (y_pos[i+1:] + delta_y)%1000.0
tmp_pos_z = (z_pos[i+1:] + delta_z)%1000.0
tmp_dir_x = x_dir[i+1:]
tmp_dir_y = y_dir[i+1:]
tmp_dir_z = z_dir[i+1:]
my_dir_x = x_dir[i]
my_dir_y = y_dir[i]
my_dir_z = z_dir[i]
pointer_x = tmp_pos_x - my_pos_x
pointer_y = tmp_pos_y - my_pos_y
pointer_z = tmp_pos_z - my_pos_z
norm = sqrt(pointer_x**2 + pointer_y**2 + pointer_z**2)
pointer_x = pointer_x/norm
pointer_y = pointer_y/norm
pointer_z = pointer_z/norm
norm_pointer = sqrt(pointer_x**2 + pointer_y**2)
norm_all = sqrt(tmp_dir_x**2 + tmp_dir_y**2)
my_norm = sqrt(my_dir_x**2 + my_dir_y**2)
#first compute the angles from the pointing vector with all the clusters besides pivot
cos_theta_pair_2D = (tmp_dir_x*pointer_x + tmp_dir_y*pointer_y)/(norm_all * norm_pointer)
cos_theta_pair_3D = (tmp_dir_x*pointer_x + tmp_dir_y*pointer_y + tmp_dir_z*pointer_z)
radius_2D = sqrt((my_pos_x - tmp_pos_x)**2 + (my_pos_y-tmp_pos_y)**2)
radius_3D = sqrt((my_pos_x - tmp_pos_x)**2 + (my_pos_y-tmp_pos_y)**2 + (my_pos_z-tmp_pos_z)**2)
dist_pair_2D = append(dist_pair_2D, radius_2D)
cos_theta_pair_all_2D=append(cos_theta_pair_all_2D, cos_theta_pair_2D**2)
dist_pair_3D = append(dist_pair_3D, radius_3D)
cos_theta_pair_all_3D=append(cos_theta_pair_all_3D, cos_theta_pair_3D**2)
#compute the angles from the pointing vectors with the pivot
cos_theta_pair_2D = (my_dir_x*pointer_x + my_dir_y*pointer_y)/(my_norm * norm_pointer)
cos_theta_pair_3D = (my_dir_x*pointer_x + my_dir_y*pointer_y + my_dir_z*pointer_z)
radius_2D = sqrt((my_pos_x - tmp_pos_x)**2 + (my_pos_y-tmp_pos_y)**2)
radius_3D = sqrt((my_pos_x - tmp_pos_x)**2 + (my_pos_y-tmp_pos_y)**2 + (my_pos_z-tmp_pos_z)**2)
dist_pair_2D = append(dist_pair_2D, radius_2D)
cos_theta_pair_all_2D=append(cos_theta_pair_all_2D, cos_theta_pair_2D**2)
dist_pair_3D = append(dist_pair_3D, radius_3D)
cos_theta_pair_all_3D=append(cos_theta_pair_all_3D, cos_theta_pair_3D**2)
# print size(cos_theta_pair_all_3D)
return dist_pair_2D, cos_theta_pair_all_2D, dist_pair_3D, cos_theta_pair_all_3D
In [438]:
In [439]:
def final_stats(cos_theta2, r_pair):
min_log_r_array = linspace(0.0,2.0,11)
max_log_r_array = min_log_r_array + 0.2
center_log_r_array = empty((0))
mean_theta2_array = empty((0))
for log_min_r, log_max_r in zip(min_log_r_array, max_log_r_array):
min_r = 10**log_min_r
max_r = 10**log_max_r
index = where((r_pair<max_r) & (r_pair>min_r))
# print min_r, max_r, size(index)
cos_theta_pair_tmp = cos_theta2[index]
mean_theta2 = mean(cos_theta_pair_tmp)
center_log_r_array = append(center_log_r_array, 0.5*(log_min_r + log_max_r))
mean_theta2_array = append(mean_theta2_array, mean_theta2)
return center_log_r_array, mean_theta2_array
In [443]:
r_pair_2D, cos_theta_2_2D, r_pair_3D, cos_theta_2_3D = get_pointing_stats(x_pos, y_pos, z_pos, x_dir, y_dir, z_dir)
In [444]:
print shape(r_pair_2D), shape(r_pair_3D), shape(cos_theta_2_2D), shape(cos_theta_2_3D)
log_r, mean_theta2 = final_stats(cos_theta_2_2D, r_pair_2D)
plot(log_r, mean_theta2)
Out[444]:
In [348]:
min_r = 3.0
max_r = 6.0
index = where((dist_pair_2D<max_r) & (dist_pair_2D>min_r))
print size(index)
cos_theta_pair_tmp = cos_theta_pair_all_2D[index]
print mean(cos_theta_pair_tmp), std(cos_theta_pair_tmp)
n, bins, patches = hist(cos_theta_pair_tmp, 20, normed=1, color = 'silver')
In [286]:
min_r = 0.0
max_r = 2.0
index = where((dist_pair_2D<max_r) & (dist_pair_2D>min_r))
print size(index)
cos_theta_pair_tmp = cos_theta_pair_all_3D[index]
print mean(cos_theta_pair_tmp), std(cos_theta_pair_tmp)
n, bins, patches = hist(cos_theta_pair_tmp, 20, normed=1, color = 'silver')
In [288]:
n, bins, patches = hist(dist_pair_2D, 20, normed=1, color = 'silver')
In [295]:
n_points = 30000
phi = 2.0*pi*random.random(n_points)
theta = arccos((random.random(n_points)-0.5)*2.0)
rand_x = sin(theta)*cos(phi)
rand_y = sin(theta)*sin(phi)
rand_z = cos(theta)
In [296]:
cos_angle_3D = (rand_x[1:]*rand_x[0] + rand_y[1:]*rand_y[0] + rand_z[1:]*rand_z[0])**2
In [298]:
n, bins, patches = hist(cos_angle_3D, 30, normed=1, color = 'silver')
In [299]:
print mean(cos_angle)
In [300]:
norm_all = sqrt(rand_x[1:]**2 + rand_y[1:]**2)
norm_one = sqrt(rand_x[0]**2 + rand_y[0]**2)
cos_angle_2D = ((rand_x[1:]*rand_x[0] + rand_y[1:]*rand_y[0])/(norm_all * norm_one))**2
In [302]:
n, bins, patches = hist(cos_angle_2D, 30, normed=1, color = 'silver')
In [313]:
n, bins, patches = hist(axis_ratio_2_1, 30, normed=1, color = 'silver')
In [375]:
def test(n_points):
a =arange(n_points)
b = arange(n_points)
c = arange(n_points)
return a,b,c
In [378]:
l, n,m = test(3)
print l, m, n
In [ ]: