2010.1/12
前回で一次処理が終わったので、今回からはその画像を使ってサイエンスをしていきたい。とりあえずは等級を測ってカラーを見ることを目指す。
単に等級を測るだけならSExtractorでAUTO_MAGを見れば良いが、カラーを測ることを考えるとApertureを固定したAperture magnitudeの方が良い。ここで問題になるのはApertureをどのくらいにとるかということだが、S/N比を高く維持するためには、点源の拡がり=PSF程度が良いと思われる。PSFは大体、星のFWHMの二倍と考えて、まず星のFWHMを測りApertureを決めてやる。
それが済んだらいよいよサイエンスに入るが、今日はとりあえずcl0024画像(perfect1222.fits)の等級毎の天体数(単位面積あたり)を測る。Luminosity Functionもどき。
☆其の一 一次処理し終わった画像(perfect1222.fits)の測光からPSF(星の拡がり具合)の決定
◎Apertureを適当に設定して星のFWHM調べる
・とりあえずsextractorかけて天体のFWHMを調べる(この段階でのApertureは適当でいい)
以下のようなMAG_AUTO,MAG_APER,FWHMを吐き出すようなパラメータファイル作る。
→112mag1.sex
112mag1.param
実行
$ sex perfect1222.fits -c 112mag1.sex
↓
結果の確認(MAG_AUTO=$4とFWHM=$7の関係を可視化)
$ gnuplot
gnuplot> plot "112mag1.cat" u 4:7 , 5 w l →こんなの(epsファイル)
→FWHM=5より少し下に水平分布してるのが星！
・適当な明るさがあり(18~16等級)、FWHM<５の星のみでリスト化
$ awk '$4>16 && $4<18 && $7<5 {print}' 112mag1.cat > 112strmag1.cat
・星だけを狙ってFWHMを丁寧に調べてやる
$ more 112strmag1.cat
15 16.4011 17.3004 16.4150 1385.599 63.655 4.68 0
94 16.9298 17.9647 16.9444 2372.658 217.584 4.89 0
112 16.0986 16.9096 16.1200 2449.045 233.658 4.41 0
129 17.1669 17.9828 17.1802 3715.242 271.549 4.43 0
208 17.6537 18.4964 17.6622 2788.159 357.005 4.62 3
216 17.8746 18.5937 17.8751 819.300 344.976 4.28 0
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・（略）・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（左から順にNo,MAG_ISO,MAG_APER,MAG_AUTO,x,y,FWHM.FLAGS）
→[2372:217]に星があるはずなのでds9上で探しそこでimexamして詳細調べる
cl> !ds9 &
cl> reset stdimage=imt8192
cl> display perfect1222.fits 1 zs- zr- z1=-20 z2=200
imexamというds9画面上で色々な情報引き出せるタスクを使う
cl> imexam
(カーソルを目標天体に合わせて"a")
# COL LINE COORDINATES
# R MAG FLUX SKY PEAK E PA BETA ENCLOSED GAUSSIAN DIRECT
2372.83 217.40 2372.83 217.40
13.38 12.15 137875. -1.493 4162. 0.08 0 2.25 4.43 4.49 4.46
（右端３つ、ENCLOSED,MOFFAT,DIRECTがFWHMを示してる）
ちなみにimexamについて
(カーソルを目標天体に合わせて"e")→contour
(カーソルを目標天体に合わせて"r")→fitting具合が見れる
("q")→imexamモード抜けられる
更にそうしたのを見た結果・・・imexamについてイジリタイ場合は
cl> epar rimexam
I R A F
Image Reduction and Analysis Facility
PACKAGE = tv
TASK = rimexam
(banner = yes) Standard banner
(title = ) Title
(xlabel = Radius) X-axis label
(ylabel = Pixel Value) Y-axis label
(fitplot= yes) Overplot profile fit?
(fittype= moffat) Profile type to fit #fittingの仕方。gausianまたはより複雑なmoffat#
(center = yes) Center object in aperture?
(backgro= yes) Fit and subtract background?
(radius = 5.) Object radius #中心から5pixelでfitting行う#
(buffer = 5.) Background buffer width
(width = 5.) Background width
(iterati= 3) Number of radius adjustment iterations
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・（略）・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
でいろいろいじれる。今回はfittype=moffatでいく。
こうした調子で112strmag1.catから一つ一つ星の座標を探してimexamしていくのもいいが、非常に面倒なので、ds9の機能を使ってもう少し楽をする。
星リストから更に天体の(x,y)座標だけのリストをつくる
awk ' {print $5,$6}' 112strmag1.cat > 112str.reg
ds9の上部で
Region→file format→XY
もう一回 Region→ load region → 112str.regを選択
で星可視化されるはず
そこ狙ってimexam,aの連発(ただし目で見て明らかにおかしいものはとばす)
cl> imexam
# COL LINE COORDINATES
# R MAG FLUX SKY PEAK E PA BETA ENCLOSED GAUSSIAN DIRECT
849.60 2485.25 849.60 2485.25
11.17 12.88 70612. -0.0055 3424. 0.07 -21 3.20 3.82 3.76 3.72
865.69 1916.24 865.69 1916.24
11.63 13.22 51670. 13.8 2132. 0.09 20 2.60 3.92 3.92 3.88
169.74 1401.37 169.74 1401.37
12.11 12.73 80845. 5.598 3202. 0.07 27 3.09 4.04 4.08 4.04
1441.04 1476.70 1441.04 1476.70
11.27 12.02 156271. 1.779 6908. 0.07 -11 3.91 3.83 3.92 3.76
1258.45 826.19 1258.45 826.19
11.47 12.49 101122. 13.46 4432. 0.06 -23 3.21 3.87 3.90 3.82
323.62 715.46 323.62 715.46
13.50 13.17 54035. 9.964 1714. 0.18 -36 2.76 4.47 4.51 4.50
819.39 344.96 819.39 344.96
11.97 13.04 60598. 2.216 2509. 0.05 -26 4.20 3.99 4.07 3.99
1922.05 516.44 1922.05 516.44
11.84 12.06 150314. 4.35 6269. 0.02 -18 3.89 3.96 4.04 3.95
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
右３列がFWHM
→星のFWHMは大体４くらい
◎Apertureを星のFWHMの2倍程度にして測光
上で星のFWHMが４くらいと分かったので、112mag1.sexのphoto_aperを4*2=8に設定してまわす。Aperture=8に固定してSN比の良い測光を行う。
→112mag+.sex
$ sex perfect1222.fits -c 112mag+.sex
→結果リスト112mag1.catの３列目にはS/N比の良いAperture Magnitudeが入ってる。
$ more 112mag+.cat
# 1 NUMBER Running object number
# 2 MAG_ISO Isophotal magnitude [mag]
# 3 MAG_APER Fixed aperture magnitude vector [mag]
# 4 MAG_AUTO Kron-like elliptical aperture magnitude [mag]
# 5 X_IMAGE Object position along x [pixel]
# 6 Y_IMAGE Object position along y [pixel]
# 7 FWHM_IMAGE FWHM assuming a gaussian core [pixel]
# 8 FLAGS Extraction flags
1 12.2754 13.0052 12.3065 2108.532 213.747 8.54 4
2 15.5282 16.7430 15.5235 3398.547 65.518 9.63 0
3 17.9392 19.1567 17.8459 2511.405 78.741 11.36 0
4 16.3707 17.5633 16.2672 436.725 35.946 8.96 0
5 17.0392 18.2219 16.9716 1745.814 35.953 15.43 16
6 17.6960 18.3718 17.6288 1843.614 22.517 9.21 16
7 20.3052 20.7145 18.7174 3434.373 20.454 11.65 18
8 18.5439 18.9015 18.2890 2427.395 19.200 5.62 0
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・（略）・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
てな感じ
☆其の二 等級毎の天体数を調べる。ただし単位面積あたりの個数
銀河団銀河以外の銀河や銀河系内の星なども入っているが一応Luminosity Functionのようなものを作りたい。
◎exposure mapをみながら一定以上の深い観測出来てる領域をきりぬき
積分時間が少ない所では当然暗い天体が検出されにくいので、真中の領域だけ切り出す
cl> display exp1222.fits 2
・exp1222.fitsを目で見るかぎり十分深く観測出来てる領域は
→[600:3250,500:1850]！
・112mag+.catから[600:3250,500:1850]内の天体データだけ取り出してリスト化
$ awk '$5>600 && $5<3250 && $6>500 && $6<1850 {print}' 112mag+.cat > 112mag2.cat
・更にmag_Autoだけをリスト化
$ awk '{print $4}' 112mag2.cat > 112mag3.cat
◎等級bin毎の天体数を数えてグラフ化
用意は整ったので（十分深い観測が出来ている領域の天体の等級だけが入ったリスト112mag3.cat）、等級bin毎の天体数を数えていく。iraf上のタスクphistogramを使えば簡単。
・等級bin毎の天体数を弾き出すタスクphist使う
cl> epar phist
I R A F
Image Reduction and Analysis Facility
PACKAGE = plot
TASK = phistogram
input = 112mag3.cat The input data #先に作った考えてる領域の天体のAuto magだけのファイルをinputとする#
(z1 = 14.) Minimum histogram intensity #14〜24等級だけで処理する#
(z2 = 24.) Maximum histogram intensity
(binwidt= 0.5) Resolution of histogram in intensity units #bin一つあたりの幅。例えば14から14.5等の天体数を数える#
(nbins = 512) Number of bins in histogram #上のbinwidtがINDEFになってる場合nbinの数だけのbinを作ってくれる#
(autosca= yes) Adjust nbins and z2 for integer data?
(top_clo= no) Include z2 in the top bin?
(hist_ty= normal) Type of histogram
(listout= yes) List instead of plot histogram?
(title = imtitle) Title for the plot
(xlabel = Data values) X-axis label
(ylabel = Counts) Y-axis label
(wx1 = INDEF) Left user x_coord if not autoscaling
(wx2 = INDEF) Right user x_coord if not autoscaling
(wy1 = 0.) Lower user y_coord if not autoscaling
(wy2 = INDEF) Right user y_coord if not autoscaling
(logx = no) Log scale x-axis?
(logy = yes) Log scale y-axis?
(round = no) Round axes to nice values?
(plot_ty= line) Type of vectors to plot
(box = yes) Draw a box around periphery of window?
(ticklab= yes) Label the tick marks?
(majrx = 5) Number of major divisions along the X axis
(minrx = 5) Number of minor divisions along the X axis
(majry = 5) Number of major divisions along the Y axis
(minry = 5) Number of minor divisions along the Y axis
(fill = yes) Fill viewport vs enforce unity aspect ratio
(vx1 = 0.) Left limit of device viewport (0.0:1.0)
(vx2 = 1.) Right limit of device viewport (0.0:1.0)
(vy1 = 0.) Lower limit of device viewport (0.0:1.0)
(vy2 = 1.) Upper limit of device viewport (0.0:1.0)
(append = no) Append to an existing plot?
(pattern= solid) Line pattern
(device = stdgraph) Output graphics device
(mode = ql)
(:go)
14.25 1
14.75 7
15.25 8
15.75 8
16.25 20
16.75 32
17.25 43
17.75 37
18.25 59
18.75 60
19.25 51
19.75 73
20.25 98
20.75 105
21.25 126
21.75 145
22.25 87
22.75 46
23.25 3
23.75 0
結果を保存してみる
cl> phist 112mag3.cat > 112hist3.cat
とりあえず可視化
$ gnuplot
gnuplot> plot "112hist3.cat" w histep →こんなの(epsファイル)
・今上で求めた数は全視野の数になってるので、視野面積で割って単位面積あたりの個数にする
今画像の視野面積を計算する。
1pixelを見込む角度はMOIRCSの場合0.117秒。また今pixel数で[600:3250,500:1850]の領域見てるのでこの領域の面積を「分^2」で表すと
S=(3250-600)*(1850-500)*0.117*0.117/3600=13.60 #3600で割ることで秒＾２→分^2に変換#
よって視野の大きさは13.60（平方分）
横軸に等級、縦軸に単位面積(1arcmin^2)あたりの個数をとってプロットさせる
gnuplot> plot "112hist3.cat" u 1:($2/13.60) w histep →こんなの(epsファイル)
暗い天体ほど数が多く、23等あたりが検出限界といったことが一目瞭然
◎銀河団のない領域と比較！
112mag3.catだけ見てても銀河団cl0024の性質が見えてこないので別の適当な領域のデータと比較してみる
・鍛冶澤先生のページから一般領域(blank field)のK-band等級カタログをとってくる。一応ここにも置いとく→Kmag_MODSwide.txt(視野:103.333arcmin^2)
（一列目から順にNo,K_AUTO,S/Nという並びになってる）
こちらも等級だけのデータのリストつくっとく(phistにかける用)
$ awk '{print $2}' Kmag_MODSwide.txt > 112Kmag_MODS.cat
先と同様にphistにかけて等級bin毎の天体数を求める
cl> phist 112Kmag_MODS.cat > 112KMODShist.cat
・先の112hist3.catと比べる
$ gnuplot
gnuplot> set logscale y #縦軸はlogスケールで#
gnuplot> plot "112hist3.cat" u 1:($2/13.60) w histep , "112KMODShist.cat" u 1:($2/103.333) w histep #あくまで1arcmin^2あたりの数を比べないといけない#
→こんなの(epsファイル)
・・・銀河団の影響で明るい天体の数に違いが出ているのがよく分かる。