Skip to content

Commit f7edd6a

Browse files
ixakalabadieixakalabadie
authored
Add files via upload
1 parent 48cdda4 commit f7edd6a

4 files changed

+1058
-0
lines changed

‎Means and Figures.ipynb

+309
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,309 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import matplotlib.pyplot as plt\n",
10+
"\n",
11+
"%matplotlib inline"
12+
]
13+
},
14+
{
15+
"cell_type": "code",
16+
"execution_count": null,
17+
"metadata": {},
18+
"outputs": [],
19+
"source": [
20+
"import healpy as hp\n",
21+
"import numpy as np\n",
22+
"from scipy import optimize\n",
23+
"import pandas as pd"
24+
]
25+
},
26+
{
27+
"cell_type": "code",
28+
"execution_count": null,
29+
"metadata": {},
30+
"outputs": [],
31+
"source": [
32+
"beta=hp.ang2vec(264.021,48.253,lonlat=True) #direction of CMB dipole"
33+
]
34+
},
35+
{
36+
"cell_type": "markdown",
37+
"metadata": {},
38+
"source": [
39+
"The newt cell reads data from the optimization, for all fluxes, and all samples. The first four arrays of each list correspond to the lowest flux, the second four to the second lowest and so on."
40+
]
41+
},
42+
{
43+
"cell_type": "code",
44+
"execution_count": null,
45+
"metadata": {},
46+
"outputs": [],
47+
"source": [
48+
"full=['full','wo0_5','noMo'] \n",
49+
"mask=['maskfullz','maskwo0_5','masknoMo'] \n",
50+
"flux=['1','5','10','']\n",
51+
"#from full-sky simulations\n",
52+
"amp_f=[] #amplitudes\n",
53+
"dire_f=[] #direction vector\n",
54+
"dire_f_lb=[] #direction in galactic coordinates\n",
55+
"#from masked-sky simulations\n",
56+
"amp_m=[]\n",
57+
"dire_m=[]\n",
58+
"dire_m_lb\n",
59+
"for f in flux:\n",
60+
" for n in full:\n",
61+
" amp_f.append(np.array(pd.read_csv('path/simul%s/%s/amp.dat'%(f,n),header=None,sep=' ')))\n",
62+
" dire_f.append(np.array(pd.read_csv('path/simul%s/%s/dir_vec.dat'%(f,n),header=None,sep=' ')))\n",
63+
" dire_f_lb.append(np.array(pd.read_csv('path/simul%s/%s/dir_lonlat.dat'%(f,n),header=None,sep=' ')))\n",
64+
" if n=='noMo':\n",
65+
" amp_f.append(np.array(pd.read_csv('path/simul%s/structure/amp_deg_full.dat'%f,header=None,sep=' ')))\n",
66+
" pix=np.array(pd.read_csv('path/simul%s/structure/dir_deg_full.dat'%f,header=None,sep=' '))\n",
67+
" vec=[]\n",
68+
" for p in pix:\n",
69+
" vec.append(hp.pix2vec(32,p,nest=True))\n",
70+
" dire_f.append(np.array(vec))\n",
71+
" lon,lat=hp.pix2ang(32,pix,nest=True,lonlat=True)\n",
72+
" dire_f_lb.append(np.array(np.hstack([lon,lat])))\n",
73+
" for n in mask:\n",
74+
" amp_m.append(np.array(pd.read_csv('path/simul%s/%s/amp12.dat'%(f,n),header=None,sep=' ')))\n",
75+
" dire_m.append(np.array(pd.read_csv('path/simul%s/%s/dir22.dat'%(f,n),header=None,sep=' ')))\n",
76+
" dire_m_lb.append(np.array(pd.read_csv('path/simul%s/%s/dir22_lonlat.dat'%(f,n),header=None,sep=' '))) \n",
77+
" if mas=='maskwo0_5noMo':\n",
78+
" amp_m.append(np.array(pd.read_csv('path/simul%s/structure/amp_deg_mask.dat'%f,header=None,sep=' ')))\n",
79+
" pix=np.array(pd.read_csv('path/simul%s/structure/dir_deg_mask.dat'%f,header=None,sep=' '))\n",
80+
" vec=[]\n",
81+
" for p in pix:\n",
82+
" vec.append(hp.pix2vec(32,p,nest=True))\n",
83+
" dire_m.append(np.array(vec))\n",
84+
" lon,lat=hp.pix2ang(32,pix,nest=True,lonlat=True)\n",
85+
" dire_m_lb.append(np.array(np.hstack([lon,lat])))"
86+
]
87+
},
88+
{
89+
"cell_type": "code",
90+
"execution_count": null,
91+
"metadata": {},
92+
"outputs": [],
93+
"source": [
94+
"#separate the directional parameters in galactic longitud and latitude\n",
95+
"l_f,b_f,l_m,b_m=[],[],[],[]\n",
96+
"for i in range(len(dire_f_lb)):\n",
97+
" l_f.append(dire_f_lb[i].T[:][0])\n",
98+
" b_f.append(dire_f_lb[i].T[:][1])\n",
99+
" l_m.append(dire_m_lb[i].T[:][0])\n",
100+
" b_m.append(dire_m_lb[i].T[:][1])"
101+
]
102+
},
103+
{
104+
"cell_type": "code",
105+
"execution_count": null,
106+
"metadata": {},
107+
"outputs": [],
108+
"source": [
109+
"#normalize direction vectors\n",
110+
"dirnorm_f,dirnorm_m=[],[]\n",
111+
"for d in range(len(dire_f)):\n",
112+
" vec_f,vec_m=[],[]\n",
113+
" for v in dire_f[d]:\n",
114+
" vec_f.append(v/np.linalg.norm(v))\n",
115+
" dirnorm_f.append(vec_f)\n",
116+
" for a in dire_m[d]:\n",
117+
" vec_m.append(a/np.linalg.norm(a))\n",
118+
" dirnorm_m.append(vec_m)"
119+
]
120+
},
121+
{
122+
"cell_type": "code",
123+
"execution_count": null,
124+
"metadata": {},
125+
"outputs": [],
126+
"source": [
127+
"#calculate angle between estimated direction and CMB direction\n",
128+
"ang_f,ang_m=[],[]\n",
129+
"for m in range(len(dirnorm_f)):\n",
130+
" ang=[]\n",
131+
" for d in dirnorm_f[m]:\n",
132+
" cos=np.dot(beta,d)\n",
133+
" ang.append(np.arccos(cos))\n",
134+
" ang_f.append(np.array(ang))\n",
135+
" ang=[]\n",
136+
" for d in dirnorm_m[m]:\n",
137+
" cos=np.dot(beta,d)\n",
138+
" ang.append(np.arccos(cos))\n",
139+
" ang_m.append(np.array(ang))"
140+
]
141+
},
142+
{
143+
"cell_type": "markdown",
144+
"metadata": {},
145+
"source": [
146+
"The following cell prints the mean of the amplitude (times 1e3) for each flux threshold and sample.\n",
147+
"The parameter of which the mean is wanted to be obtained can be changed."
148+
]
149+
},
150+
{
151+
"cell_type": "code",
152+
"execution_count": null,
153+
"metadata": {},
154+
"outputs": [],
155+
"source": [
156+
"pm=u\"\\u00B1\" #define plus-minus sign\n",
157+
"for i in amp_m:\n",
158+
" print(str(round(np.mean(i*1000),2))+pm+str(round(np.std(i*1000),2)))"
159+
]
160+
},
161+
{
162+
"cell_type": "markdown",
163+
"metadata": {},
164+
"source": [
165+
"Below different models of how the plots have been obtained are presented."
166+
]
167+
},
168+
{
169+
"cell_type": "markdown",
170+
"metadata": {},
171+
"source": [
172+
"The following cell creates the amplitude histograms."
173+
]
174+
},
175+
{
176+
"cell_type": "code",
177+
"execution_count": null,
178+
"metadata": {},
179+
"outputs": [],
180+
"source": [
181+
"masks=['S>1 '+uJy,'S>5 '+uJy,'S>10 '+uJy,'S>22.8 '+uJy]\n",
182+
"fig, ax=plt.subplots(2,2,figsize=(14,10))\n",
183+
"fig.subplots_adjust(wspace=0.2,hspace=0.35)\n",
184+
"c=0\n",
185+
"for j in range(2):\n",
186+
" for i in range(2):\n",
187+
" kin=amp_m[0+4*c]\n",
188+
" wol=amp_m[1+4*c]\n",
189+
" nomo=amp_m[2+4*c]\n",
190+
" struc=amp_m[3+4*c]\n",
191+
" ax[j][i].set_title(masks[c],fontsize=24) \n",
192+
" ax[j][i].hist(struc,bins=np.arange(0,0.007,0.0003),edgecolor='black',fc=(0,1,0,0.5),lw=0.5)\n",
193+
" ax[j][i].hist(nomo,bins=np.arange(0,0.007,0.0003),edgecolor='black',fc=(1,1,0,0.5),lw=0.5)\n",
194+
" ax[j][i].hist(kin,bins=np.arange(0,0.007,0.0003),edgecolor='black',fc=(0,0,1,0.5),lw=0.5)\n",
195+
" ax[j][i].hist(wol,bins=np.arange(0,0.007,0.0003),edgecolor='black',fc=(1,0,0,0.5),lw=0.5)\n",
196+
" ax[j][i].axvline(0.00462,c='black',lw=4)\n",
197+
" ax[j][i].tick_params(axis='both', which='major', labelsize=14)\n",
198+
" ax[j][i].set_ylabel('counts',fontsize=18)\n",
199+
" ax[j][i].set_xlabel('amplitude',fontsize=18)\n",
200+
" c+=1\n",
201+
"fig.savefig('amplitude_hist.pdf',bbox_inches='tight')"
202+
]
203+
},
204+
{
205+
"cell_type": "markdown",
206+
"metadata": {},
207+
"source": [
208+
"The next model has been used for direction scatter plots of all samples.\n",
209+
"In the for loop a transformation is made in the parameters to fit this modules coordinate intervals:\n",
210+
" $l=[ -\\pi , \\pi ]$ and $b=[ -\\pi/2 , \\pi/2 ]$\n",
211+
"for the structure part another the same transformation has to be done but it is more complex because of the\n",
212+
"randomness of the directions."
213+
]
214+
},
215+
{
216+
"cell_type": "code",
217+
"execution_count": null,
218+
"metadata": {},
219+
"outputs": [],
220+
"source": [
221+
"fig=plt.figure(figsize=(8,6))\n",
222+
"ax=fig.add_subplot(111, projection=\"mollweide\")\n",
223+
"la,ba=[],[]\n",
224+
"for i in range(2):\n",
225+
" la.append(np.deg2rad(360-l_m[i]))\n",
226+
" ba.append(np.deg2rad(b_m[i]))\n",
227+
"ax.scatter(l_m[3],b_m[3],c='g',marker='.',label='structure')\n",
228+
"ax.scatter(la[0],ba[0],c='b',marker='.',label='kinematic')\n",
229+
"ax.scatter(la[1],ba[1],c='r',marker='.',label='kinematic w/o local structure')\n",
230+
"ax.scatter(np.deg2rad(360-264.021),np.deg2rad(48.253),c='black',marker='o',label='CMB')\n",
231+
"ax.grid()\n",
232+
"ax.set_xlabel(r\"$l$\",fontsize=20)\n",
233+
"ax.legend()\n",
234+
"h=ax.set_ylabel(r'$b$',fontsize=20)\n",
235+
"h.set_rotation(0)\n",
236+
"ax.set_title('S>1 '+uJy,fontsize=28)\n",
237+
"ax.set_xticklabels([150,120,90,60,30,0,330,300,270,240,210,180])\n",
238+
"ax.tick_params(axis='both', which='major', labelsize=14)\n",
239+
"fig.savefig('mollweide_dirs.pdf',bbox_inches='tight')"
240+
]
241+
},
242+
{
243+
"cell_type": "markdown",
244+
"metadata": {},
245+
"source": [
246+
"Longitude, latitude and theta histograms have been obtained following the next model."
247+
]
248+
},
249+
{
250+
"cell_type": "code",
251+
"execution_count": null,
252+
"metadata": {},
253+
"outputs": [],
254+
"source": [
255+
"colors=[(0,0,1,0.5),(1,0,0,0.5),(1,1,0,0.5),(0,1,0,0.5)]\n",
256+
"fig, ax=plt.subplots(3,4,figsize=(24,20))\n",
257+
"ax[0][0].set_title(r'$S\\!>1\\,\\mu$Jy',fontsize=26)\n",
258+
"ax[0][1].set_title(r'$S\\!>5\\,\\mu$Jy',fontsize=26)\n",
259+
"ax[0][2].set_title(r'$S\\!>10\\,\\mu$Jy',fontsize=26)\n",
260+
"ax[0][3].set_title(r'$S\\!>22.8\\,\\mu$Jy',fontsize=26)\n",
261+
"for i in range(ax.shape[0]):\n",
262+
" for j in range(ax.shape[1]):\n",
263+
" li=l[i+4*j]\n",
264+
" if i!=3:\n",
265+
" ax[i][j].hist(li,bins=np.arange(220,300,10),edgecolor='black',fc=cols[i],lw=0.5)\n",
266+
" ax[i][j].axvline(np.mean(li),c='brown',lw=3)\n",
267+
" ax[i][j].axvline(264.02,c='black',lw=3)\n",
268+
" else:\n",
269+
" ax[i][j].hist(li,bins=np.arange(0,360,45),edgecolor='black',fc=cols[i],lw=0.5)\n",
270+
" ax[i][j].axvline(np.mean(li),c='brown',lw=3,label=str(round(np.mean(li),1)))\n",
271+
" ax[i][j].legend(fontsize=15)\n",
272+
" ax[i][j].tick_params(axis='both', which='major', labelsize=18)\n",
273+
" if j==0:\n",
274+
" ax[i][j].set_ylabel('counts',fontsize=22)\n",
275+
" if i==3:\n",
276+
" ax[i][j].set_xlabel(r'$\\l$ / deg',fontsize=22)\n",
277+
"fig.savefig('longitude.pdf',bbox_inches='tight')"
278+
]
279+
},
280+
{
281+
"cell_type": "code",
282+
"execution_count": null,
283+
"metadata": {},
284+
"outputs": [],
285+
"source": []
286+
}
287+
],
288+
"metadata": {
289+
"kernelspec": {
290+
"display_name": "Python 3",
291+
"language": "python",
292+
"name": "python3"
293+
},
294+
"language_info": {
295+
"codemirror_mode": {
296+
"name": "ipython",
297+
"version": 3
298+
},
299+
"file_extension": ".py",
300+
"mimetype": "text/x-python",
301+
"name": "python",
302+
"nbconvert_exporter": "python",
303+
"pygments_lexer": "ipython3",
304+
"version": "3.7.4"
305+
}
306+
},
307+
"nbformat": 4,
308+
"nbformat_minor": 2
309+
}

0 commit comments

Comments
 (0)