CLASS是gildas中的一个软件,之前对CLASS有所接触,但仅限于相关的文件格式转换(http://blog.sciencenet.cn/home.php?mod=space&uid=117333&do=blog&id=303542)。最近课题组需要用CLASS处理数据,所以我们在阅读CLASS手册(http://www.iram.fr/IRAMFR/GILDAS/doc/pdf/class.pdf)以及一个网上简易教程(http://www.mpifr-bonn.mpg.de/staff/bparise/gildas/2010-pety-class-demo-max-planck.pdf)、张智昱博士的网站(http://159.226.71.202/wiki/FAST_2013)一些网页(http://www.astro.umass.edu/~fcrao/library/manuals/map.html)以及德令哈观测站编写的CLASS使用手册(http://www.radioast.csdb.cn/tools.php)的同时做一些笔记,形成此CLASS极简教程。有一些命令是新版本的用fortran 90重写的CLASS里才有的,在这种情况下,命令前面的提示符用LAS90>表示。
本手册不定期更新,牛人无须看。
1. FITS文件(.fits)转换为gildas文件(.gdf, .30m) 天文中一些不同领域定义了不同的文件格式,不同的数据处理软件有时候也有自己的数据格式,比如说gildas相关软件用的就是自己的文件格式(后缀为.gdf)。不同的领域交流时通常使用标准的FITS文件,所以将FITS文件转换为gildas文件是使用CLASS的一个基础操作。.30m文件是CLASS文件默认的后缀名。
最简单的转换方法是使用vectorfits工具(vector包里的fits工具)。在shell里敲class,回车,打开class后,使用如下命令
LAS> vectorfits spectra.fits to spectra.gdf
就可以把FITS文件spectra.fits转换为gildas文件spectra.gdf。不过要注意,FITS文件需要满足一定要求,否则转换出来的gildas文件也不能为CLASS所用。比如,如果spectra.fits是数据块(data cube),那么转换出来的文件CLASS是无法直接使用的。
CLASS能接受的文件是光谱文件,如果有FITS格式的光谱文件(spectra),那么可以如上述将文件转换为gildas格式,也可以进行如下操作将fits文件转换为CLASS能使用的文件(.30m),其中spectra.fits是一个FITS格式的光谱文件。
LAS> set angle sec
LAS> lasfits read spectra.fits
LAS> sic delete spectra.30m
LAS> file out spectra.30m new
LAS> write 1
LAS> file in spectra.30m
LAS> get f
LAS> plot
注意到,读入fits文件用的是lasfits,和vectorfits是不同的。上述操作第一步是设置角度的单位,就文件转换而言是不必要的。第二步就是读入FITS格式的光谱文件。第三步是删除已有的CLASS文件,第四步是写入新的CLASS文件。之后就是读入新生成的CLASS文件,画图。
对于数据块fits文件(比如叫做qhz_cube.fits),可以使用lmv读入、转换。具体如下
LAS90> file out qhz.14m m
LAS90> lmv qhz_cube.fits
这样,数据块fits文件qhz_cube.fits就转换为calss文件qhz.14m了。
2. gildas文件(.30m)转换为FITS文件(.fits)
这个操作就是将前面一个操作倒过来,不过,在这个情况下,vectorfits就不顶用了,得用lasfits,将上面一个操作中的后一个方法倒过来。关于此操作,还可以参看(http://iram.fr/IRAMFR/GILDAS/doc/html/class-html/node15.html )。
LAS> set angle sec
LAS> file in spectra.30m
LAS> find
LAS> get f
LAS> set fits mode spectrum
LAS> lasfits write spectra.fits
其中get f是读出第一条光谱。
在CLASS90里可以将一系列光谱文件转换为FITS二进制表文件(生成FITS数据块文件还不知道如何实现)。
LAS90> file in spectra.30m
LAS90> find /range -200 200 -200 200
LAS90> fits write spectra.fits /mode index
3. gildas文件(.30m)转换为文本文件(.dat)
最简单的光谱文件就是在文本文件里写两列数,第一列是波长或者频率,第二列是强度。将CLASS文件转换为文本文件可以方便其他一些程序进行处理。具体操作如下,其中spectra.30m是输入的CLASS文件
LAS>file in spectra.30m
LAS>find
LAS>get f
LAS>sic output spectra.dat
LAS>for i 1 to channels
LAS: say ’rx[i]’ ’ry[i]’ /format g12.4 g12.4
LAS: next
LAS>sic output
4.读入文本文件(.dat)中的光谱并画图假设spectra.dat是含有光谱的文本文件。LAS> dev xl w
LAS> greg1column x 1 y 2 /file ’’spectra.dat’’
LAS> model y x /regular
LAS> plot
5. 减基线
减基线是谱线数据处理中经常需要进行的操作,不过光看手册的话会比较迷惑,参考对Effelsberg望远镜数据处理的简介(http://www.astro.uni-bonn.de/~rcbruens/tutorials/effberg/HI.html),我们总结出以下操作。
LAS> file in spectra.30m
LAS> find
LAS> get f
LAS> plot
LAS> set cursor on
LAS> set window
这个时候,当把鼠标放到画图窗口的时候会出现一条横线和一条竖线,交点是鼠标的位置。不看说明的话,一般都不知道怎么退出这种模式(答案是按“e”键)。这个模式是让你到画图窗口标明一下强发射线的位置,鼠标到某条发射线(吸收线)的左边,按空格键,到右边,按空格键,以此类推,对所有发射线(吸收线)进行此操作。这个过程你会感到心虚,因为你看不到任何变化。不要急,完成上面步骤之后,按“e”键退出。(也可以直接指定窗口的位置,比如LAS> set window -120 -100 0 20。其中的数就是各窗口左侧和右侧的坐标。)
LAS> draw window
此时就可以看到在下方出现一个矩形,标明发射线的位置(反射性所在的频率窗口(波长窗口))。然后就可以拟合基线了(比如用三个自由参数的,也就是二次多项式拟合,并画出基线)。
LAS> base 3 /plot
这时就看到拟合的基线了。然后
LAS> plot
就可以看到减完基线的谱线了。
结合下面提到的用CLASS90将光谱排好画出来
LAS90> find
LAS90> load
LAS90> plot /index
可以在二维上定义减基线的窗口
LAS90> set window /poly 1 ! Define 1 polygon
6. 构建数据块
成图观测最终能方便使用的数据时数据块,即三维数组,三维分别是赤经、赤纬和速度(频率)。但是成图观测的原始数据是一系列光谱,把这些光谱按坐标排好才能形成数据块。具体构建数据块的步骤如下(其中demo.30m(http://www.mpifr-bonn.mpg.de/staff/bparise/gildas/demo.30m)是相关数据)。
LAS90> file in demo.30m
LAS90> find
LAS90> list
这一步是看看文件中有哪些观测。
LAS90> find /tel 30M-V01-A100
寻找一批标号为30M-V01-A100的观测,然后检查一致性
LAS90> consistency
如果数据满足一致性,就可以进行以下操作
LAS90> let name thecube
LAS90> let type lmv
LAS90> table ’name’ new
LAS90> xy_map ’name’
至此就生成了数据块文件thecube.lmv,然后可以查看
LAS90> go view
LAS90> vectorfits outfile.fits from infile.lmv
7. 对多条谱线进行减基线操作
通常一个数据文件中有很多条谱线,
如果每条谱线都手动进行减基线操作将会非常费时费力。通常一次观测的谱线窗口相差不大,所以可以进行批处理。下面是从读入数据、设定窗口、到循环减基线的脚本。最后画出了所有谱线。
file in spectra.30m
set window 10 20
find /number 59 100
for i 59 to found
get i
base 3
next
stamp
将上面的语言写到一个叫baseline的文件中,在class程序里
LAS> @baseline
就可以执行了。
8. 谱线拟合
减完基线就可以进行谱线拟合了。读入数据画图的步骤和前面相同(qhz.14m是青海站观测的数据http://www.radioast.csdb.cn/tools.php)
LAS> file in qhz.14m
LAS> find
LAS> get f
LAS> plot
下面就是谱线拟合了,假定有四条谱线
LAS> line 4
此时就会出现和set window之后类似的情况,不过要简单一点。四条谱线共需要八个位置来定义,鼠标点八次(每条谱线的左右),然后自然就回到命令窗口。
LAS> min
就开始拟合,完了可以查看
LAS> vis
9. 存ps、eps文件画完图,如果不希望退出程序图就消失,那么就需要把图存在文件中。一般把图存成ps或eps格式,发表文章时也可以用。使用下面的命令就可以完成这个操作
LAS> hardcopy spectra.ps /dev ps fast
用来存ps文件。
LAS> hardcopy spectra.eps /dev eps color
用来存eps文件。
10. 画p-v图画p-v图需要取数据块中沿某个空间方向的一片数据,这可以通过strip命令实现(不过首先要通过find命令找到这一片数据),具体如下。
LAS> file in qhz.14m
读入数据,然后设定角度的单位为角分
LAS> set angle min
结下来用find命令寻找数据片
LAS> find /offset * 5
找y方向偏离5角分的一片数据(x-v二维数组),然后用strip命令将此数据存成一个文件
LAS> strip qhz.gdf
随后就可以将这片数据读入,画p-v图了。
LAS> define image a qhz.gdf read
设定画图的范围
LAS> limit 0 50 * * /rgdata a
然后画出数据,对数据进行重新的线性标度
LAS> greg2plot /scal linear 0 20
然后画出等值线
LAS> level 1 to 10
LAS> rgmap /abs 2
然后就是画轴和做标记
LAS> lab "Velocity (km/s)" /x
LAS> lab "gDDEC. (arcmin)" /y
LAS> axis xl /unit r
LAS> axis xu
LAS> axis yl /unit m
LAS> axis yr
将光谱排好画出来CLASS90还可以这么做
LAS90> find offset * 5
LAS90> load
LAS90> plot /index
11. 定义数组在CLASS中可以定义数组,并对数组进行运算,这对于数据处理是非常必要的。定义10个元素整数数组(数组的指标从1开始)可以用如下命令
LAS> define integer a[10]
如果要定义双精度数数组(我不知道单精度数怎么弄,囧rz……)可以用如下命令
LAS> define double b[10]
定义二维数组可以用如下命令
LAS> define double c[10,10]
多维数组也可以类似定义。经本人研究,最多能定义4维,比IDL的8维的上限要弱小一些。不过亲~,4维一般够用了吧。
12. 数组赋值、查看等操作 数组赋值可以用let命令
LAS> let a[1] 1
查看数组值可以用exam命令
LAS> exam a[1]
13.去除RFI(坏通道)
基本想法就是用上面提到的并排画光谱的操作将光谱排好,如果有RFI(坏通道),那么某几个通道的值就会很大,超过采样最大值。这个之后需要把这些通道用高斯噪声填满,这就是去RFI(坏通道)的基本方法。
LAS90> file out a100-fill single /over
LAS90> find
LAS90> for ient 1 to found
LAS90> get next
LAS90> fill -120 -100 60 75 /noise ! File contaminated channels with Gaussian noise
LAS90> write
LAS90> next ient
14. 思考题 如何画任意方向的p-v图?
spectra.30m
15. 画积分强度图
本段是向张智昱、潘之辰学习而来的。首先要对一个数据文件作积分(在CLASS里):
LAS90> file in co65_final.apex
LAS90> find
LAS90> get first
LAS90> set unit v
LAS90> plot ! 查看一下谱线的速度范围
LAS90> print moment 0 500 /out moment.dat
LAS90> exit ! 退出CLASS
打开GREG,读入之前得到的数据文件
GREG> column x 2 y 3 z 4 /file moment.dat
GREG> set box m
GREG> lim /rev x
GREG> random 100 /bl -1 !将数据转换为二维图像
GREG> plot !画积分强度图
GREG> box /abs !画坐标框
GREG> wed ! 画颜色棒
张智昱指出,在CLASS里是可以直接调用GREG里的函数的,在函数名前面加上greg就可以了。
16. 快速查看谱线常用命令
LAS90> file in xxx.cso ! 后缀名一般是观测用的望远镜的名称,格式是一样的。
LAS90> find /source xxx /telescope "yyy" ! 寻找用yyy仪器观测的xxx源的数据
LAS90> list
LAS90> get nnn ! 读入第nnn条光谱
LAS90> set nomatch ! 不匹配坐标
LAS90> set u v ! 将横坐标的单位设为速度(而不是频率)
LAS90> average ! 把上面找到的所有谱线平均(对找弱谱线尤其有用)
LAS90> set mode x a b ! 取一段x坐标a到b的数据
LAS90> set cursor on ! 设置鼠标可以在窗口读数
LAS90> set window ! 设置窗口,用鼠标读数(原来是靠左键点,现在靠按n,退出读数读取还是按e)
LAS90> plot
附录
一、脚本
1. fitscube2bur.class
FITS数据块文件转换为bur文件,用法(假设FITS文件为data.fits):
LAS90> @fitscube2bur data
=======================================
sic delete &1.bur
file out &1.bur m
lmv &1.fits
=======================================
2. bur2fitscube.class
将bur文件转换为FITS数据块(先转换为lmv文件,然后再转换为FITS),用法(假设bur文件为data.bur):
LAS90> @bur2fitscube data
=======================================
file in &1.bur
find
consistency
let name thecube
let type lmv
table ’name’ new
xy_map ’name’
vectorfits thecube.fits from thecube.lmv
=======================================
转载本文请联系原作者获取授权,同时请注明本文来自钱磊科学网博客。
链接地址:http://blog.sciencenet.cn/blog-117333-567008.html
上一篇:思考题(十三)赤道附近的望远镜有什么优势?
下一篇:科学家、科研技术员和科研工人